HD 74438: A tight spectroscopic quadruple as possible progenitor of sub-Chandrasekhar Type Ia supernovae


 Stars often form in multiple systems and may follow a complex evolution involving mass transfer and collisions, leading to mergers that are possible progenitors of Type Ia supernovae (SNe) [1, 2]. The progenitors of such explosions are still highly debated [3]. While binaries have received much attention so far, higher-order stellar systems show a wide variety of interactions especially in tight systems, like long-term gravitational effects playing a key role in triple (where they are called von Zeipel-Lidov-Kozai , [4, 5], hereafter ZLK, oscillations) and quadruple systems. Here we report on the properties of the first spectroscopic quadruple (SB4) found within a star cluster: the 2+2 hierarchical system HD 74438 [6]. Its membership in the open cluster IC 2391 makes it the youngest (43 My) SB4 discovered so far. The eccentricity of the 6 y outer period is 0.46 and the two inner orbits, with periods of 20.5 d and 4.4 d, and eccentricities of 0.36 and 0.15, are not coplanar. Using an innovative combination of ground-based high resolution spectroscopy [7, 8, 9, 10] and Gaia/Hipparcos astrometry [11, 12, 13, 14], we show that this system is undergoing secular interaction that likely pumped the eccentricity of one of the inner orbit higher than expected for the spectral types of its components. We compute the future evolution of HD 74438 by considering gravitational dynamics, stellar evolution, and binary interactions [15], and show that this system is an excellent candidate progenitor of sub-Chandrasekhar Type Ia supernova through white dwarf (WD) mergers. This specific type of SNIa better accounts for the chemical evolution of iron-peak elements in the Galaxy [16].

The integrated spectral type of the HD 74438 system is A5 [17]; it belongs to one of the closest young clusters, IC 2391 containing 254 stars [18], located at 146 +8 −7 pc in the Vela constellation, aged 43 +15 −7 My [19,20] and of solar metallicity [21]. It was expected [22] that this system should at least be a triple star because it lies 0.9 mag above the main sequence (MS) of the Hertzsprung-Russell (HR) diagram of the parent cluster. It was detected as a four-component spectroscopic system (SB4) only recently [6] within the ground-based, large spectroscopic Gaia-ESO Survey (GES) [7,8]. Follow-up observations with HRS/SALT [9] and HERCULES/UCMJO [10] allowed us to derive the orbital solution independently for the brighter (components A and B) and the fainter (components C and D) pairs (Table 1 and Methods). The systemic velocities of each pair (v 0 in Table 1) bracket the velocity of the parent cluster [23] (14.98 ± 0.17 km s −1 ), a first hint towards the pairs AB and CD being gravitationally bound. The component mass ratios imply similar brightnesses in each of the inner pairs, well in line with their detection as SB2s.
We adjusted synthetic spectra (see Methods) on two HRS/SALT spectra displaying four well-separated components. One of these observed spectra is shown in Fig. 1, together with the synthetic spectra of the four components and the resulting combined synthetic spectrum. The derived effective temperatures and spectral types are listed in Table 1.
Masses, luminosities and radii are then inferred from a stellar isochrone [24] of 43 My at solar metallicity, since the system components necessarily lie (provided no mass transfer has occurred) on the cluster main sequence (see Fig. 2 and Methods). An important validation of our method is that the spectroscopic mass ratios thus derived are in good agreement, within the uncertainties, with the dynamical ones (Methods and Table 1). The dynamical mass ratio q = (M C +M D )/(M A +M B ) = 0.692 ± 0.003 could thus be derived, and appears to be among the lowest q of SB4s known so far (see Supplementary Information). Adding together the four component luminosities, a total luminosity of 15.7 ± 1.8 L ⊙ is obtained, which is in excellent agreement with the independent determination provided by the Gaia DR2 luminosity [11] (15.3 -0.1 +0.2 L ⊙ , Fig. 2).
We can deduce, following Eqs. 1, 5 and 6 of Methods, the orbital inclinations with respect to the sky plane : i AB = 52.5 ± 1.5°for the massive pair and i CD = 84.0 ± 0.9°for the low-mass pair ( Table 1). The two SB2 orbits are thus far from being coplanar. The CD pair should moreover show eclipses that could however not be detected from the inspection of the TESS [25] photometric data. Thanks to archival ESO observations and the recent RV follow-up, it became possible to constrain the mutual orbit of the two pairs. The period and eccentricity of the wide pair are found to be 5.7 y and 0.46 (Methods and Table 1). Moreover, thanks to the knowledge of the masses, the inclination of the wide pair on the sky is found to be 73.2°± 2.7°. However, the current mutual inclinations (i AB / AB-CD and i CD / AB-CD ) are not known because they require the knowledge of the longitudes of the ascending nodes (Eq. (10) of Methods), which can only be determined by astrometric or interferometric measurements, should they be able to turn HD 74438 in a visual system. In the Methods section, we show that Hipparcos [13,14] and Gaia [11,12] astrometric data succeed in doing so for the outer orbit, giving access to the corresponding longitude of the ascending node. However, a similar application of the astrometric method to the inner orbits is impossible because of the closeness of the pairs, and in this case, even an interferometric imaging approach is challenging. This innovative combination of Hipparcos and Gaia astrometric data will be especially useful in the context of the forthcoming Gaia data releases.
Concerning its birth condition, HD 74438's properties may be compared to those of the much younger quadruple system discovered in the star-forming core Barnard 5 and consisting of a young protostar and three gravitationally-bound dense gas condensations [26]. The latter system has a much wider separation, but it is recognised that separations tend to decrease over time as the protostar grows in mass and dynamically interacts with the local gas reservoir [27]. These two systems thus represent different stages in the evolution of multiple-star systems within their birth environment.
HD 74438 is a rare case with a relatively short outer period of 5.7 y when compared with the other known 2+2 quadruples from the Multiple Star Catalogue [28]. Such systems are uncommon in OCs; actually, our system has the shortest outer period known so far in an OC, as shown in the left panel of Fig. 3, μ Ori [29] being the previous record holder (see also Supplementary Information). It is currently unclear whether this results from an observational bias or from the difficulty to form them [30]. Orbital shrinking caused by gas accretion and dynamical friction mechanisms that took place during core fragmentation [27] seems to be required to account for such short outer periods.
The high level of characterisation of this quadruple system allows us to use it as a test bench to get insight into the secular evolution of multiple hierarchical systems [31]. Long-term gravitational effects in which orbit-averaged torques change the orbital angular momenta and eccentricities play a key role in triple and quadruple systems; they can indeed drive the eccentricity to high values in quadruple systems. Despite its importance, the ZLK mechanism was ignored for many years, but was revived some 10-20 years ago for triple systems with the detection of the eccentric planet 16 Cyg B [32] or the close to perpendicular orbits in the Algol system [33,34]. The double astrometric binary μ Ori has also been reported to experience ZLK cycles [29]. Within triple systems, the inner and outer orbits exchange angular momentum, which leads to mutual inclination and eccentricity oscillations on timescales much longer than the orbital periods. To explore the evolution of HD 74438 (up to 10 Gy in the future), we use the Multiple Stellar Evolution code [15] that takes into account a wide range of processes, most importantly gravitational dynamics, stellar evolution, and binary interactions such as mass transfer and common-envelope evolution. We employ a Monte-Carlo approach to sample a set of realisations of HD 74438 taking into account the observational uncertainties (see Methods). The first interesting result is the fact that ZLK oscillations do indeed occur in the CD pair (right panel of Fig. 4), due to angular momentum transfer between the CD and AB-CD pairs, for a very wide range of mutual orbital inclinations between these pairs. Incidentally, the merging of the AB pair is even more probable (left panel of Fig. 4) because the ZLK timescale of AB is shorter than the one of CD (see Methods). Hence for systems like HD74438 one or more mergers occur in almost 50% of the simulated cases. We show in the Methods the relative fractions of all mergers in our simulations (collisions, as well as CE evolution). Collisions between main-sequence stars are common, as well as CE events involving giant stars. Also possible are WD-WD mergers.    . We highlighted blue quadruples belonging to clusters as cross-matched with SIMBAD. HD74438 appears to have the smallest outer period. The solid line represents the condition for 2+2 quadruples to be dynamically stable [5]. Right: Location of the two inner SB2 (AB in blue, CD in red) in the eccentricity-period diagram, compared with A-type primaries (matching the AB pair; grey dots in the top panel) and G-type primaries (matching the CD pair; grey dots in the bottom panel) of the SB9 catalogue [59] of spectroscopic binaries and doubly eclipsing binaries [60] in yellow circles. In the case of the e AB amplitude, there is a large gap near 90°inclinations, since these have merged during their previous evolution. The e CD -plot shows that high eccentricity amplitudes are possible for a large range of initial mutual inclinations i CD / AB-CD . Due to quadruple dynamics, the range of mutual inclination angles is much enhanced, as can be seen by noting the colours: R's close to unity are required to get high amplitudes for small mutual inclinations, whereas if i CD / AB-CD is close to 90°, high amplitudes can also be attained regardless of R values.

Methods
Observations and data reduction. The details of the observations are summarized in Table 2.
The spectra were obtained with three different high-resolution spectrographs: (i) UVES/VLT, (ii) HERCULES/UCMJO, (iii) HRS/SALT, and with one medium-resolution spectrograph (GIRAFFE/VLT). In addition, archival data from ESO with the mid-resolution spectrograph GIRAFFE/VLT [53] allows us to add three epochs in 2004 and help to constrain the longer period but only for the AB pair because C and D components are not resolved with GIRAFFE.
Cross correlation function (CCF) calculations. The CCFs were computed by cross-correlating normalized observed spectra with a unique template built as follows. The initial guess for the template was based on the spectral type attributed to the unresolved multiple system (A2 [17]). We thus adopted a Kurucz model atmosphere [56] with T eff = 9000 K. Nevertheless a lower temperature (T eff = 7000 K) provided more contrasted CCFs. Therefore we adopted a Kurucz model atmosphere of T eff = 7000 K, log g = 4.0 and [Fe/H] = 0. We computed a synthetic spectrum from 3700 to 8900 Å to cover the three spectrographs spectral range. Strong lines were masked (Balmer lines, H&K Ca II, Ca II IR triplet, and Na I D lines). We then used Detection of Extrema (DOE) [6] to extract the positions and depths of 6123 lines, building a comb of rectangular functions 0.01 Å-wide, and with a height corresponding to the line depth. The CCFs were then computed by cross-correlating the normalized spectra and the comb of rectangular functions.
Radial-velocity measurements. The identification of the number of radial velocity (RV) components, their positions and the RV uncertainties was performed with DOE [6]. First, the third derivative of the CCF was used to identify the number of components. Second, a multi-Gaussian fit was performed on the CCF. Examples of fits to multiple CCF peaks are illustrated in Fig. 6 for HRS/SALT and HERCULES/UCMJO. Three components are always visible and a fourth one can also be distinguished in most CCFs; they are labelled A, B, C, D by decreasing CCF peak intensity. While the assignment of the A and B components is straightforward, several trials and errors were necessary for a correct assignment of the C and D components because their peak intensities are similar, especially among the HRS spectra where the sparse sampling does not allow to easily follow each component over time.
Orbital parameters of the inner pairs. The HD 74438 system is a 2+2 hierarchical system. The time interval between UVES and HERCULES/UCMJO + HRS/SALT data is about 5 y, enough for the wider AB-CD pair to imprint a secular trend; therefore only the HERCULES/UCMJO and HRS/SALT data sets were used to derive the AB and CD short-period orbits. We also assumed that the instrumental zero-point RV offset between the two spectrographs is negligible. In Fig. 7 , with temperatures ranging from 4000 to 10000 K (with a step of 250 K), convolved with a Gaussian of standard deviation of 5 km s -1 . No significant broadening was detected beyond the instrumental one. The best composite synthetic spectrum was retained on the basis of the smallest standard deviation of the (observed -calculated) residuals. The adopted T eff of each component is the average of the values obtained from the two spectra quoted above. The T eff uncertainties were computed by adding quadratically the error on the mean (since the T eff determination was performed on two spectra, providing slightly different results for the C and D components arising from the small differences between the synthetic and observed spectra at different orbital phases) and the grid step. As an illustration of the quality of the fit, Fig. 1  Radii are also derived. Finally, the spectroscopic mass ratios for each pair and for the AB-CD system, i.e. M CD /M AB , are derived. The astrophysical parameters of HD 74438's individual components are listed in Table 1. The large temperature range from the Gaia-ESO Survey (7600 ± 750 K) well brackets the temperatures of the brightest-pair components while the Gaia temperature [70] (7423 ± 52 K) agrees well with component B.

Inclinations and separations.
Combining orbital and astrophysical parameters, it is possible to deduce the inclination of each SB2 pair: (1) where Kp, P, Mp, e and q are the RV amplitude of the primary, the orbital period, the mass of the primary of the pair, the eccentricity and the mass ratio of the secondary over the primary. We obtain i AB = 52.5 ± 1.5°and i CD = 84.0 ± 0.9°. Knowing the inclinations, it is then possible to derive the geometrical separation of each pair: (2) where the p index stands for the primary of each pair. We obtain a AB = 0.21 ± 0.25 au and a CD = 0.068 ± 0.004 au, where uncertainties on inclination dominate. . We now compare this luminosity to the one obtained when adding the four component luminosities. Since the temperatures of the four components have been spectroscopically determined, the cluster PARSEC isochrone at 43 My directly provides their luminosities. Summing up these luminosities leads to a total luminosity of 15.7 ± 1.8 L ⊙ for the combined system. It is represented by a horizontal line in Fig. 2, because the temperature of the SB4 considered as a single star (for instance as given by the GES temperature of 7600 ± 750 K; blue shaded area in Fig. 2) has no physical meaning. The combined luminosity is thus in excellent agreement with the Gaia DR2 one. The above value can be compared to the luminosity also derived from the Gaia DR2 catalogue [11, 18], but dereddening the photometry, using the usual relation:  Table 4 from [62]. We finally obtain from Eq. (3) an extinction-corrected Gaia luminosity of L = 21.1 ± 9.5 L ⊙ , as represented in Fig. 2 by the brown dot. The uncertainty on the derived luminosity is obtained by propagating uncertainties from G, , and BC G and is admittedly large, because of the uncertainty on the bolometric correction. The cluster stars whose membership probability is higher than 90% are also represented as grey dots in Fig. 2, with effective temperatures as provided by the GES iDR4 [19].
Parameters of the outer orbit. We retrieved observations with GIRAFFE spectrograph taken in 2004. Combined with our recent follow-up with HRS/SALT and HERCULES/UCMJO, this allows us to cover more than half of the orbital phase. We computed the center of mass velocity of each inner pair and fit a Keplerian orbit. The fitted RV curve is given in Fig. 9. We obtained a period of about 5.7 y. The dynamical stability of hierarchical systems imposes the outer period to be at least approximately five times longer than the longest inner period [30], which is clearly the case here. The high eccentricity is not surprising for such a young system. For the derivation of the inclination of the outer pair, we used Eq. (1) while for the semi-major axis, we directly used the third Kepler law, the values being reported in Table 1.

Dynamical masses.
An independent way to derive the stellar masses (besides the location of the stars in the HR diagram) resorts to the velocity semi-amplitudes through the set of equations: The third relation leads to (6) which allows to derive three masses once one has been fixed. Adopting the lower possible spectroscopic mass value M A (i.e. derived from the location in the HR diagram; Table 1, row M (a)), we obtained dynamical masses consistent with the spectroscopic ones for all the components. The masses derived in this way are listed in Table 1 in row M (b). Because both subsystems are SB2, the dynamical mass ratios can be derived as q = M s /M p = K p /K s , where s and p denote the secondary and primary components, and K is the RV amplitude.

Astrometric constraints. Astrometry from Gaia DR2 [11], Gaia eDR3 [12], Tycho-1 [13]
and Tycho-2 [14] provide constraints on the orientation on the sky of the relative orbit AB-CD, as we now show. As may be seen in the Δ μ α * and Δ μ δ rows of Table 3, the differential proper motion of the photocentre of HD 74438 with respect to the centre of mass of the cluster varies according to both the epoch of the observations and the time span covered by the position observations used to derive the proper motion. That differential motion is compatible with zero for the Tycho-2 data since these long time-span data (about one century-long) average out the AB-CD orbital motion occurring on a much shorter time scale (5.7 y, as indicated in Table 1). On the contrary, the much shorter time spans of the Tycho-1, Gaia DR2 and eDR3 proper motions reveal the orbital motion on the sky. This motion is easily modelled by the Thiele-Ines constants [63] A, B, F, G, which allow to express the position (x, y) of a component on the sky with respect to a reference point (both will be specified below) as follows: with (x', y') the coordinates in the orbital plane: and v is the true anomaly, E the eccentric anomaly, e the orbital eccentricity, r the length of the radius-vector and a the semi-major axis (whose exact definition will be specified below). After derivating the above relations with respect to time, the proper motion components write being the Thiele-Ines constants divided by n a", where n = 2π / P and a" is the orbital semi-major axis expressed in arcsec, i.e., a" = a au ὼ where ὼ is the parallax in arcsec. The angles ω and Ω are the argument of periastron and longitude of the ascending node, respectively. In the following, since we are dealing with astrometric and spectroscopic orbits, these angles refer to the absolute orbit of the considered component (or photocentre) around the centre of mass of the considered system.
To go further, it is thus necessary to specify that the orbit probed by the observed differential proper motion is that of the photocentre of the AB-CD pair around the centre of mass of the system. The semi-major axis of the corresponding photocentric orbit (a phot AB-CD ) is related to that of the relative AB-CD orbit (a rel ) by the relation  (8) may be converted successively in a relation for the apparent semi-major axis on the sky (a" phot AB-CD = a phot AB-CD ὼ ) and in a relation for the relative orbit proper motion (μ rel = v" rel = 2 π a" rel / P = n a" rel ). For HD 74438, a" rel = 36.37 mas and μ rel = 40.37 mas y -1 . Thus μ phot AB-CD = μ rel (κ -β) = 13.03 mas y -1 .
From the above considerations, it results that the proper-motion components observed by the astrometric satellites (and listed in Table 3) correspond to those of the photocentre of the AB-CD pair around the centre of mass of the system, i.e., μ phot AB-CD , so that the value of a" entering Eqs. (7a-b) is actually a" phot AB-CD . Therefore, one may define Ẋ ≡ μ α * / μ phot AB-CD and Ẏ ≡ μ δ / μ phot
The quantities Ẋ and Ẏ are proper-motion observables, obtained at some given epoch, whereas dξ/dt and dη/dt, as well as ẋ' and ẏ', depend on time through the orbital motion. Equations (9a-b) then allow to derive Ω. Since the Gaia DR2 proper motion is the one obtained on the shortest time span among all those listed in Table 3, it most accurately represents the orbital motion and will be used in the astrometric analysis. The orbital ephemeris obtained with the orbital elements listed in Table 1 where i 1 , i 2 are the inclinations of the long-and short-period orbit on the sky, and Ω 1 , Ω 2 their respective longitudes of the ascending nodes. A definite evaluation of Φ is not possible, however, in the absence of the knowledge of the longitude of the ascending node of the short orbit. Nevertheless, given the large eccentricity (e = 0.15) for this short-period system (P = 4.4 d; compare its location in the e -P diagram of similar systems with GV primaries in the bottom-right panel of Fig. 3), it seems likely that the ZLK phenomenon is at work through the gravitational interaction of the CD pair with AB-CD.
Birth scenario discussion. Multiple stars are thought to be produced according to the following scenarios: (i) cluster fragmentation [65], i.e. core fragmentation for the outer (large period) subsystem and disk fragmentation for the inner ones (namely the 2 shorter period pairs) or (ii) dynamical and resonant capture. Several observables of HD 74438 contradict the capture scenarios. Some simulations of dynamical evolution within young clusters in the Solar Neighborhood have initial conditions (N = 182, simulation time = 50 My [66]) that match well with our case (IC 2391, with 254 objects, is 43 My old). These simulations, including primordial single, binary, triple but not quadruple stars, under-produce the latter compared to the observed fractions by a factor of three, suggesting that a primordial population of quadruples must be present. Moreover, some magneto-hydrodynamics simulations coupled with N-body, stellar evolution and binary interaction codes [67] show that the binary fraction drops far below the observed statistics for low-mass stars, presumably because primordial binary formation was not included. These results tend to support a formation of lower-mass multiple systems like HD 74438 predominantly by primordial core and disk fragmentations rather than by capture within the cluster.

ZLK oscillation timescales.
To zeroth-order approximation, the ZLK oscillation timescale depends on the masses, the periods and the eccentricity of the outer period. Using Eq. 1 of [68] we estimated these timescales to be about 240 y for the pair AB around the center of mass of CD and to be 620 y for the pair CD around the center of mass of AB. These correspond to 10⁴ to 10⁵ secular ZLK oscillation timescales since the birth of the cluster, enough to allow interactions (triple common envelope evolution, merger, ejection, etc.) to occur.
Relativistic precession. Depending on the system parameters, general relativity effects may cause pericenter precession, preventing the inner pair from reaching high eccentricities and thereby avoiding close encounters. We estimated the relativistic timescales (~0.1 and~0.01 My for the AB and CD orbit, respectively) to be much longer than the ZLK time scales quoted above. Therefore, relativistic precession does not decrease the maximum eccentricity reachable by the two inner pairs in HD 74438.
Future evolution of HD 74438. The MSE code models the evolution of multiple-star systems of arbitrary configurations [15], taking into account gravitational dynamics, stellar evolution, and binary interactions such as mass transfer and CE evolution. Also included are simple recipes for 'triple' interactions such as triple CE evolution, when a tight binary star enters the envelope of a giant star. We employ Monte Carlo methods to sample a set of realisations of the observed system taking into account the observational uncertainties. Specifically, for each system we sample all parameters (the four masses and five orbital elements for each of the three orbits, i.e., 19 parameters in total) from Gaussian distributions centred at the observed values, and with standard deviations given by the observed error bars. An exception to the latter applies to the longitudes of the ascending nodes, Ω i , which are not all observationally constrained and which were sampled from flat distributions. With this Monte Carlo approach, we sample N MC = 10⁴ systems and evolve them with MSE until a system age of 10 Gy was reached, or until the maximum allowed CPU wall time of 20 hr was reached, which occurred for about 67% of the systems. Given the highly compact nature of HD74438, the system is computationally prohibitively expensive to evolve (even using the secular approach in MSE), since the secular time-scales are short compared to 10 Gy. Most first mergers occur at an early age of ∼10-100 My, significantly younger than the attained age of most of the systems for which the CPU wall time was exceeded (peaking near 10 9 y, see Fig. 10). This justifies the imposed maximum wall time. A summary of our results is given in Table 4, which shows the fractions per occurring physical processes, number of remnants and outcomes of our simulations. We distinguish between different events, i.e., each system can be attributed one or more events, and where we distinguish between non interacting systems, those with one or more mergers (including both CE evolution and direct collisions), CE evolution, direct collisions, dynamical instability, and RLOF of a tertiary star onto an inner binary. In the latter case, mass transfer can proceed stably (not likely), or unstably leading to triple CE evolution. About half of realisations do not lead to interactions such as mergers, whereas the other half involves interactions such as stellar mergers, and triple CE evolution. The high fraction of triple Roche Lobe Overflow (RLOF) and triple CE (close to 40%) is much higher compared to the expectation for the entire population of isolated triple systems [69]. We classify the final outcome of the system in terms of the number of remaining objects (which can be stars or compact objects), as well as in terms of their hierarchical configuration. The non-interacting systems correspond to stable quadruple systems (fraction ∼ 0.5), which are the only final configuration with four remnants (e.g., the 'Triple+Single' outcome does not occur). When three remnants remain (fraction ∼ 0.1), they can be either in a stable triple (likely), or a binary with an unbound single (not likely). If there are two remnants (fraction ∼ 0.1), the fractions of them being bound and unbound to each other are ∼0.03 and ∼0.06, respectively. A single remnant occurs for a fraction of ∼ 0.2 of realisations, whereas there are no cases with no remaining remnants. In almost all instances when the tertiary fills its Roche lobe around an inner binary, the subsequent evolution is unstable and leads to triple CE evolution. In addition, the simulation show that there is a clear preference for interacting systems (in particular, mergers and triple CE systems) to have i AB / AB-CD close to 90°, whereas the non-interacting systems favor initial i AB / AB-CD further from 90°, with peaks near the smallest i AB / AB-CD around 20°, and the largest around 120°, as shown in Fig. 11. This can be easily understood by noting that highly mutually inclined systems will lead to high eccentricities due to strong secular evolution [42]. The overwhelming majority (∼ 98%) of first collisions indeed occur in orbit AB. On the other hand, i CD / AB-CD shows less distinct features among the different outcomes, with the exception that non-interacting systems favor smaller i CD / AB-CD compared to interacting ones (i.e., further separated from 90°). Fig. 12 illustrates one possible future evolution of HD 74438 towards a WD with a sub-Chandrasekhar mass: the AB pair rapidly collides after 170 My because secular evolution drives eccentricity close to unity. The merger remnant, a 3.3 M • MS star, subsequently evolves and fills its Roche lobe around the companion binary at 540 My. The outcome of the triple CE is an unstable system in which a collision quickly occurs between the two components of the inner binary. A WD+MS binary remains. As the MS companion evolves, it fills its Roche lobe at around 2 Gy, ultimately causing a merger. The final remnant is a single 1.3 M • WD.

Supplementary information
The Five septuples have also been reported with most of them being part of moving groups or associations. Due to the strong statistical biases, any conclusion on the frequency of multiples in the field and in OC is still premature. Large on-going and future surveys will allow us to increase the statistics of these interesting high-order multiples. The vertical shift of the spectra is arbitrary. Spectral fitting was performed using the two spectra marked with an asterisk and plotted in bold.         Table 3. Proper motions in right ascension (RA, including the cos De factor) and declination (De).
The differential proper motion Δμ α * and Δμ δ denote the differential stellar motion with respect to the cluster. The differential proper motions from Tycho-1, Gaia DR2, and Gaia eDR3 reflect the orbital motion of the ABCD photocentre around the centre of mass of the cluster.  Table 4. Outcome fraction per physical process (several processes can actually occur for a given simulated system), number of remnants and final outcomes of MSE simulations. Statistical error bars are given (based on Poisson statistics). 'Mergers all' means that any kind of merger occurred at any point, irrespective of whether or not an ejection event happened. Dynamical instability represents instability triggered by stellar evolution (e.g. orbital expansion due to wind mass loss). In some cases, the outcome was indeterminate from the simulation data in relation to the CPU wall time being exceeded.