Planet Hunters TESS IV: a massive, compact hierarchical triple star system TIC 470710327

ABSTRACT

1 INTRODUCTION

Despite their intrinsic rarity implied by the initial mass function (IMF; see, e.g. Salpeter 1955; Bastian, Covey & Meyer 2010; Dib, Schmeja & Hony 2017), massive stars (M ≥ 8 M_⊙) provide radiative, dynamical and chemical feedback to their environment, driving evolution on a wide range of scales. The physical processes responsible for the formation, evolution, and death of massive stars, however, are not well understood (e.g. Zinnecker & Yorke 2007; Tan et al. 2014). The study of these processes is further complicated by the fact that such stars often have nearby stellar companions (Sana et al. 2012), which can affect their properties and evolution at all stages of their lives. Large-scale spectroscopic, interferometric, and high contrast imaging surveys of OB stars have demonstrated that most, if not all, massive stars are formed in a binary or higher order multiple system (Kiminki & Kobulnicky 2012; Sana et al. 2013, 2014; Aldoretta et al. 2015; Moe & Di Stefano 2017; Maíz Apellániz et al. 2019; Rainot et al. 2020; Bodensteiner et al. 2021)

The identification and characterization of massive stars in multiple systems is crucial to discriminate between different formation and evolution scenarios, to place constraints on theoretical models and to understand the interactions and effects of multiple stars. Sana et al. (2012) and de Mink et al. (2014) have already shown that the complex interactions between tides, angular momentum exchange, and stellar evolution in binary and higher order multiples fundamentally affect the evolution and thus the end product of nearly 70 per cent of all early-type massive stars. In compact binary orbits, for example, processes including mass transfer, exchange of angular momentum and stellar mergers can open up new evolutionary pathways and end-products such as X-ray binaries, γ-ray bursts, stellar mergers, and gravitational wave events (Sana et al. 2012; de Mink et al. 2013). The addition of a nearby third body further complicates the evolution. Interactions between three stars can induce different evolutionary pathways through von Zeipel–Kozai–Lidov cycles (ZKL; von Zeipel 1910; Kozai 1962; Lidov 1962; Naoz 2016; Ito & Ohtsuka 2019), or result in tertiary-driven mass transfer that can lead to mergers, exotic common-envelope systems, close double or triple degenerate systems, or contribute to the population of walk-away and run-away systems in our galaxy (Antonini, Toonen & Hamers 2017; Renzo et al. 2019; Stephan et al. 2019; Leigh et al. 2020a; Glanz & Perets 2021; Hamers et al. 2021b).

Due to the high intrinsic brightness of massive stars, detecting and characterizing non-eclipsing lower mass close companions using spectroscopy alone is challenging. Detailed characterization of close companions, therefore, often relies on the detection of eclipses or other dynamical effects caused by the presence of a close companion, such as using radial velocity (RV) observations or observations of light travel time effects via eclipse-timing variations (ETVs). Observations of the latter were made possibly with the advent of space-based photometric surveys with long time-base, high-precision, and high-cadence observations such as CoRoT, Kepler, K2, and TESS. These space-based missions have enabled the detection of dozens of new triple and higher order multiple systems through ETVs or multiply eclipsing events (Conroy et al. 2014; Marsh, Armstrong & Carter 2014; Borkovits et al. 2015; Hajdu et al. 2017; Rappaport et al. 2017; Li et al. 2018; Sriram et al. 2018; Borkovits et al. 2021). Detailed modelling of such systems can lead to the determination of the absolute masses of the stellar components in some cases, or minimally the derivation of mass ratios, even in the absence of eclipses (Borkovits et al. 2016).

Further considerations beyond light travel time effects need to be made for triple systems where the inner stellar binary is orbited by a third body with a longer period, known as a hierarchical triple system. Such a configuration can result in measurable dynamical perturbations to the orbit of the inner binary. If the orbital periods and the separations involved are short enough, these effects can be studied with a combination of ETVs (caused by light travel time effects and direct third-body perturbations) and/or RV observations. Measurements of both ETVs and RVs allows for the determination of precise stellar mass ratios and orbital parameters, including the mutual inclination between the inner and the outer orbit of the triple, which is thought to be indicative of the formation history.

In this paper, we present a new compact, hierarchical triple system identified in TESS data consisting of one O- and two B-type stars. The system, which shows large ETV and RV variations, contains an inner ∼1.1 d eclipsing binary and a massive O9.5-B0.5V tertiary orbiting a common centre of mass in a ∼52 d orbit. The discovery of the system and the data are discussed in Sections 2–4. Section 5 outlines the analysis of the photometric and spectroscopic data and Section 6 discusses the system configuration, stability, formation and possible future evolution scenarios. Finally, the conclusions are presented in Section 7.

2 THE TARGET AND ITS SURROUNDINGS

TIC 470710327 (BD + 61 2536, TYC 4285-3758-1, V = 9.6 mag, parallax = 1.06 ± 0.24 mas, distance = 950 ± 220 pc) was initially identified as an early B-type star by Brodskaya (1953) and as a short period eclipsing binary with a period of 1.1047 d using photometric data obtained with the 0.25-m Takahashi Epsilon telescope in Mayhill, New Mexico, USA (Laur et al. 2017). The target is not a known member of a cluster or OB association (Laur et al. 2017). Two epochs of speckle interferometric measurements revealed a close companion at ∼0.5 arcsec, with the position angle and distance of the companion advancing from θ = |$306{_{.}^{\circ}}2$| and ρ = 0.533 arcsec in 1987.7568 (Hartkopf et al. 2000) to θ = |$303{_{.}^{\circ}}7$| and ρ = 0.502 arcsec in 2003.9596 (Hartkopf, Mason & Rafferty 2008).

In order to calculate the magnitude differences between TIC 470710327 and the ∼0.5 arcsec companion, hereafter TIC 470710327′, we performed speckle imaging using the Zorro instrument on the 8.1-m Gemini South telescope on Cerro Pachón, Chile (Matson, Howell & Ciardi (Howell et al. 2011; Matson, Howell & Ciardi 2019). Observations were carried out on 2020 August 15 using the two-colour diffraction-limited optical imager with 60 msec exposures in sets of 1000 frames. The 5σ detection sensitivity and the speckle reconstructed image are shown in Fig. 1. The data confirmed that the companion star is located at an angular separation of 0.529 arcsec with a position angle of |$304{_{.}^{\circ}}5$|⁠, which is in agreement with previous observations (Hartkopf et al. 2000, 2008). The data showed that the companion has a magnitude difference of Δm = 1.17 mag at 562 nm and of Δm = 1.13 mag at 832 nm.

Figure 1.

Contrast curves showing the 5σ detection sensitivity and the speckle reconstructed image for filters centred on 562 nm (blue) and 832 nm (red).

Open in new tabDownload slide

To characterize the target’s surroundings and to quantify the light contribution of nearby stars, we queried the early Gaia Data Release 3 catalog (eDR3; Gaia Collaboration et al. 2021). This search revealed a bright (V = 11.6, ΔT_mag = 2.4) nearby star at a separation of ∼22 arcsec (LS I + 61 72), as well as a further six stars with a ΔT_mag < 5 within 100 arcsec of the target (listed in Table 1). The light contribution of these stars to our photometric observations will be discussed in Section 3.

Overall, the triple system TIC 470710327 presented in this paper consists of a P₁ ∼ 1.1 d eclipsing binary with a P₂ ∼ 52 d non transiting tertiary. Hereafter, the two stars in the inner, short period binary will be referred to as stars A and B, while the tertiary on the wide, outer orbit will be referred to as star C. The 0.5 arcsec companion star will be referred to as TIC 470710327′. Using the Gaia eDR3 list of sources in a 2 arcmin radius to determine the local field density, and the magnitude contrasts given by the speckle observations, we show that the spurious association probability between TIC 470710327′ and TIC 470710327 is ∼ 1 × 10⁻⁵. This means that statistically there is an association between the target and the nearby companion. Assuming that the Gaia distance to this target (∼950 pc; Gaia Collaboration et al. 2018) is accurate, an angular separation between TIC 470710327′ and TIC 470710327 of 0.5 arcsec corresponds to a physical separation of approximately 500 au with an orbital period on the order of 1500-2000 yr, depending on the mass of the companion. Following this, we will for the remainder of this paper assume that their physical separation is too large for TIC 470710327′ to have any dynamical effect on the triple system on the timescales of our observations. Furthermore, we note that both the astrometric excess noise (significant at 4000 σ) and the RUWE parameter (11.9) for TIC 470710327 are large. This is to be expected for unresolved multiple systems, but can contribute to an unreliable astrometric determination (Ziegler et al. 2018; Lindegren et al. 2021). Conversely, the 22 arcsec companion has relatively small values for the RUWE parameter and astrometric excess noise, indicating that it is likely a single object.

3 PHOTOMETRIC DATA AND ETVS

3.1 TESS

TIC 470710327 was identified as a potential multiple system by citizen scientists taking part in the Planet Hunters TESS (PHT) citizen science project (Eisner et al. 2021). PHT, which is hosted by the Zooniverse platform (Lintott et al. 2008, 2011), has engaged nearly 30 000 registered volunteers in the search for planetary transit signals in 2-min cadence light curves obtained by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015). In brief, each light curve is seen by 15 volunteers who identify times of transit-like events. Once a volunteer has classified a target, they are able to discuss the target on a discussion forum, and flag interesting systems to the PHT science team using searchable hashtags. TIC 470710327 was flagged as an interesting system on the discussion forum on 2019 December 6,¹ due to the light curve containing multiple periodic signals: P₁ ∼ 1.10 d, P₃ ∼ 9.97 d, and P₄ ∼ 4.01 d, as shown in the bottom panel of Fig. 2. We note that for the remainder of this paper we assume that the stars associated with the P₃ and P₄ signals are not close enough to the triple to induce observable ETVs and that these signals do not originate from star C, as discussed further in Section 6.1.

Figure 2.

Top panel: full TESS light curve of TIC 470710327 obtained during sectors 17,18, and 24. The dashed lines on the x-axis show a split in the axis. Bottom panel: light curve phase folded on P₁ = 1.1 d, P₃ = 9.97 d, and P₄ = 4.09 d signals. Each panel represents one signal, where the other two signals had been masked out. Combined these three signals make up the light curve seen in the top panel. P₃ and P₄ are assumed to not be dynamically associated with the triple system presented in this paper, as discussed in Section 5.

Open in new tabDownload slide

TIC 470710327 was monitored by TESS in the 2-min cadence data during Sectors 17, 18, and 24 of the nominal mission. The full TESS data set, displayed in the top panel of Fig. 2, spans ∼230 d, with the first subset covering nearly 60 d continuously, followed by a large gap of 140 d, before the second ∼30 d subset of observations. Visual inspection of the full TESS light curve reveals a clear ∼1.10 d eclipsing binary system (Fig. 2). The presence of full eclipses, combined with clear points of ingress and egress implies that the inner binary has a high inclination (i₁), close to 90º.

The pixel aperture used to extract the light curve, as determined by the TESS pipeline at the Science Processing Operations Center (SPOC, Jenkins et al. 2016), is displayed as the red outline in Fig. 3. The orange circles indicate the position of nearby stars that have a TESS magnitude difference ΔT_mag< 5 from the target star, as queried from Gaia eDR3 (Gaia Collaboration et al. 2021). Considering all sources with a ΔT_mag < 5 within 100 arcsec of the target (listed in Table 1), we expect TIC 470710327 (including the 0.5 arcsec companion, TIC 470710327′) to contribute 88–89.5 per cent of the total light in the TESS light curve.

Figure 3.

Average flux per pixel around TIC 470710327 obtained by TESS during Sector 17. The orange dots show all neighbouring stars with V_mag < 5ΔT_mag from the target star, as queried from Gaia eDR3 (Gaia Collaboration et al. 2021). The red outline shows the TESS aperture that was used to extract the flux in Sector 17.

Open in new tabDownload slide

Similarly, using the magnitude contrast derived from the speckle imaging, we estimate that TIC 470710327′ contributes ∼26 per cent of the observed light at 832 nm. Thus, given the ∼10.5–12 per cent light contribution of the neighbour at 22 arcsec, we estimate that TIC 470710327 and TIC 470710327′ contribute ∼65 per cent and ∼23 per cent of the total light observed by TESS, respectively. As we show in Section 6.1, the signals P₃ ∼ 9.97 d and P₄ ∼ 4.01 d cannot be hosted by the tertiary object in TIC 470710327, but likely originate from this composite contaminating light. For the remainder of the paper, we will thus not consider their influence/contribution. Finally, we attempted to use both smaller apertures as well as centroid motions to determine the direction in which the photocentre moved during the P₃ ∼ 9.97 d and P₄ ∼ 4.01 d eclipses. We found no discernible motion in the centroid positions and negligible change in eclipse depths when different apertures were used. This either suggests that all signals originate from a single TESS pixel or that these signals have too small of a light contribution to result in a detectable centroid motion.

3.2 Additional photometric observations

In addition to the TESS data, there are 531 archival photometric observations obtained with the 0.25-m Takahashi Epsilon telescope in Mayhill, New Mexico, USA between 2011 and 2013 (Laur et al. 2017); 108 photometric observations taken over 600 d by the SuperWASP-N camera located on La Palma, Canary Islands (Pollacco et al. 2006); as well as 1784 photometric observations obtained since 2012 by the ASAS-SN network of telescopes (Shappee et al. 2014; Kochanek et al. 2017). We also obtained 456 photometric observations in the Johnson Cousins Bu band between 2020 March 21 and 30, using a 0.41-m RC telescope with an Optical Guidance Systems and an SBIG STL-6303E camera, located at the Sandvreten Observatory, Sweden. Finally, we obtained 1415 photometric observations using the Las Cumbres Observatory (LCO) global network of fully robotic 0.4-m/SBIG and 1.0-m/Sinistro facilities. The LCO data were reduced and calibrated using the standard LCO Banzai pipeline. For both the LCO and the Sandvreten Observatory data, we performed aperture photometry for TIC 470710327 and five comparison stars using the open source sep package (Bertin & Arnouts 1996; Barbary 2016)

These archival and new photometric measurements significantly increased the baseline of the observations, allowing us to refine the period of the inner binary to P₁ = 1.104 686 ± 0.000 004 d. The facilities used to obtain these measurements have significantly smaller pixel scales than TESS, such that the extracted light curves do not include TIC 470710327′. This allowed us to confirm the stability of the dominant periodic signal, P₁ = 1.10 d, and verify that this binary signal does not originate from TIC 470710327′. Due to the sparse sampling of these five data sets, we are unable to use these observations to investigate the stability or the origin of either the P₃ nor the P₄ signals.

3.3 Eclipse-timing variations

Using the TESS data alone, we determined the deviations from this ephemeris following the methodology outlined by Li et al. (2018). In brief, we determined the eclipse regions of the primary and secondary eclipses by extracting the minima of the second derivative around the time of the eclipses in the phase folded light curve. These translate to the phases of ingress and egress. Prior to this, we masked the P₃ = 9.9733 d signal from the light curve. The lower amplitude signal from P₄ was not removed due to the risk of introducing spurious signals.

Next, we generated a model of the primary eclipse by fitting a trapezoid (the shape found to best represent the eclipses) to the smoothed, phase folded and subsequently binned light curve. This model was then fit to each individual primary transit where the only two free model parameters were the time of eclipse and the slope of an underlying linear trend. The latter was to allow for systematic effects that change the slope of the eclipse. The same methodology was independently carried out for the secondary eclipses. The individual fits to all eclipses, including both the primary and the secondary eclipses of the P₁ = 1.1047 d signal, were optimized using a Markov chain Monte Carlo (MCMC) approach, using the open source software exoplanet (Foreman-Mackey et al. 2021). The observed minus calculated times of eclipse (O – C), which show variations on the range of ∼± 6 min, are listed in Tables A2 and A3 for the primary and secondary eclipses, respectively. The periodicity in the O – C curves, of P₂ ∼52 d, is also observed in the radial velocities (Section 4), and is thus attributed to the tertiary star of the triple system, as discussed further in Section 5.

4 SPECTROSCOPIC OBSERVATIONS AND RADIAL VELOCITY EXTRACTION

We obtained 24 spectra between 2020 January 31 and October 5 using the hermes spectrograph (R ∼ 85 000, Raskin et al. 2011) on the 1.2-m Mercator telescope at Observatorio del Roque de los Muchachos at Santa Cruz de la Palma, Canary Islands, Spain. The spectra were reduced (extracted, order-merged, wavelength-calibrated) using the local hermes pipeline (Raskin et al. 2011), and subsequently normalized using a spline fit (Abdul-Masih et al. 2021). The spectra show strong He ii 4686 lines and no sign of He ii 4542 lines, indicating that the dominant signal originates from an early B star. The relative strengths of the He i, He ii, and Mg lines indicates that the dominant spectral contribution is consistent with that from an O9.5-B0.5V star (Sota et al. 2011). The RV shifted and median combined spectrum is shown in Fig. 4, where the bottom panels show the relative depths of pairs of diagnostic spectral lines from Sota et al. (2011) used to determine the spectral classification of the star. This spectral classification agrees with the previous estimates of the spectral class by Brodskaya (1953) and Laur et al. (2017). A spectral class of O9.5-B0.5V nominally corresponds to a 14–17-M_⊙ star (Harmanec 1988; Martins, Schaerer & Hillier 2005; Silaj et al. 2014).

Figure 4.

Top panel: Shifted and median combined hermes spectra of TIC 470710327 A (S/N > 80). Solid vertical blue lines denote He i lines, dashed vertical blue lines denote He ii lines, dotted red lines denote metal lines, and the shaded grey region denotes a diffuse interstellar band. Bottom panel: each panel displays a set of two lines used for diagnostic determination of spectral classification by Sota et al. (2011).

Open in new tabDownload slide

We searched for additional signals in the spectra using Least Squares Deconvolution (LSD, Donati et al. 1997; Tkachenko et al. 2013). In brief, this technique assumes that all lines in a spectrum have a common underlying profile with varying depths depending on the particular line. This common profile is recovered by deconvolving the lines within a particular wavelength region with a line-list template with associated line depths. This method was generalised by Tkachenko et al. (2013) to allow for multiple components in a spectra and to allow for each component to draw from a different line list. Here, we calculated LSD profiles of the Hermes spectra over 4300–5200 Å, using helium, carbon, nitrogen, and silicon lines whose rest wavelengths were computed from the Vienna Atomic Line Database (VALD-II, Kupka et al. 1999). We restricted ourselves to this range to minimize the potential noise contribution from weak lines. The resulting profiles did not reveal the presence of any additional components. The radial velocities of each spectrum were computed as the centre of gravity of the LSD profiles and are listed in Table A1. Inspection of the RVs reveal the same periodic signal, of P₂ = 52 d, as seen in the ETVs (Section 3.3). As the P₂ = 52 d signal is seen in both the ETVs and the RVs, we confirm that the system is comprised of a P₁ = 1.1047 d eclipsing binary (stars A and B) and a P₂ = 52 d outer O9.5-B0.5V tertiary (star C) orbiting around a common centre of mass.

5 JOINT ETV AND RV MODELLING

In this section, we describe the joint modelling of the ETV and the RV signals assuming that they are physically associated as a triple system. There are two main effects responsible for deviations from strict periodicity in eclipse timings: the light travel time effect (LTTE; geometrical contribution) and the dynamical effect. The former is a result of a change in projected distance from the centre of mass of the binary to that of the triple. The dynamical effect, on the other hand, results from physical changes in the orbit of the binary system due to the gravitational influence of the third body (Borkovits et al. 2003).

In order to best exploit the complementary and overlapping information in the RV and ETV data, the two data sets were modelled jointly. This allowed us to simultaneously fit for, and better constrain, parameters that appear in both models (see Table 2). The joint analysis made use of the open source software packages exoplanet and pymc3(Salvatier, Wiecki & Fonnesbeck 2016; Foreman-Mackey et al. 2021). The optimal parameter values and their uncertainties were calculated as the median and 67.8 per cent highest posterior density of the marginalized posterior distributions. The priors and extracted values for all sampled parameters are given in Table 2. The best-fitting model, as constructed from the values in Table 2, for the RVs and ETVs are shown in the left- and right-hand panels of Fig. 5. We note that the residuals in the ETVs and RVs show no evidence for additional periodicities.

Figure 5.

Joint MCMC model of the hermes RV data (left-hand panel) and the extracted ETVs (right-hand panel). The model parameters are presented in Table 2. The bottom panels show the residuals of the best fit. The overall ETV fit is a linear addition of the dynamical effect (small dashed grey line) and the light travel time effect (large dashed red line).

Open in new tabDownload slide

6 SYSTEM CONFIGURATION, STABILITY, AND EVOLUTION

6.1 Configuration

The combination of RV and photometric data revealed a dynamically interacting triple system, comprised of a close 1.1047 d eclipsing binary (stars A and B) with a massive companion on a wide, non-eclipsing, ∼52 d orbit (star C). From the dynamical modelling (presented in Section 5), we derived a mass ratio of the inner binary to the tertiary of q = m_AB/m_C = a_Csin i₂/a_ABsin i₂ = 0.7776. Given the spectral classification of star C, of O9.5V-B0.5V, the tertiary star has a mass in the range of 14–17 M_⊙, meaning that the combined mass of stars A and B is in the range of 10.9–13.2 M_⊙. By considering both the spectroscopic mass range of the tertiary and the estimated mutual inclination i_m, we can derive limits on the inclination of both the orbit of inner binary i₁ and the orbit of the tertiary i₂ on the sky. From these considerations, we find that i₁ lies in the range: |$77{_{.}^{\circ}}9$|≤ i₁ ≤ 90° and i₂ lies in the range: |$62{_{.}^{\circ}}5\,$|≤ i₂ ≤ |$71{_{.}^{\circ}}1\,$|⁠. The high inclination of i₁ agrees with the observed flat bottom, i.e. full, eclipses observed for the short 1.1 d orbit. Conversely, from the presence of flat-bottomed eclipses in the TESS light curve, we can infer that the inner binary has a large inclination: i₁ > =80°. By combining this inferred range with the mutual inclination derived from the ETV modelling, we arrive at an estimate for the inclination of the tertiary on the sky: i₂ ∈ 59°−74°. Using this, we can infer physical ranges for |$m_{\rm AB}\in 8.8{-}15.6\,{\rm M}_{\odot }$| and |$m_{\rm c}\in 10.1 - 20.3\, {\rm M}_{\odot }$|⁠. This wide range is in agreement with the mass ranges inferred from the spectral type of the tertiary star. An overview of the triple system is presented in Fig. 6. We note that while the triple and the 0.5 arcsec companion star, TIC 470710327′, are statistically associated, the companion is not assumed to have a detectable effect on the observed dynamics of the triple given its separation and short period covered by the observations.

Figure 6.

Schematic of the triple system. Relative sizes of the orbits are not to scale and for simplicity the orbits are depicted as circular.

Open in new tabDownload slide

In addition to the eclipsing P₁ = 1.1047 d binary signal, the TESS light curve contains two further periodic eclipsing signals with P₃ = 9.9733 d and P₄ = 4.092 d. The bottom panel of Fig. 2 displays the TESS light curve phase folded on P₃ (middle panel) and P₄ (right-hand panel). We note that due to the lack of visible ‘secondary’ eclipses, we cannot distinguish between periods of P₃ and P₄ or twice those values. We can, however, determine that all three sets of eclipses are of different objects, as the morphology of all of the eclipses are constant in time, and points of overlap between the 1.10 and 9.97 d signal are reproduced as linear additions of the different eclipse signals. Furthermore, we see no dynamical evidence in the RVs or the ETVs of either the 9.97 or the 4.09 d signals being part of the same system. Using Equation 3, we can show that a dynamically interacting body with P₃ = 9.97 d would produce an amplitude in the ETVs that is at least (9.97/1.10)² = 82 times larger than the dynamical amplitude induced by the P₂ = 52 d tertiary. We see no evidence of this signal in the data. Using the same argument we also find no evidence of the P₄ = 4.09 d signal being associated with the tertiary star of the triple. We note that P₁ = 1.10 d and P₃ = 9.97 d is close to a 1:9 resonance (P₃/P₁ = 9.028). While this could indicate that the signals are related, we do not have sufficient data at present to further investigate their association.

As noted above, the presence of two additional bright sources in the TESS aperture, at separations of 0.5 and 22 arcsec, overall contributing over ∼10 per cent of the light in the aperture, make it difficult to determine the exact source of all of the periodic signals seen in the TESS light curve. Through the RVs and ETVs, however, we can associate the triple system with the brightest O9.5V-B0.5V star in the data. This is further corroborated by the data obtained with the Takahashi Epsilon telescope, which has a pixel scale of 1.64 arcsec, which allowed us to confirm that the P₁ = 1.10 d signal lies on the O9.5V-B0.5V target star and not the 22 arcsec companion star that lies within the same TESS aperture.

6.2 SED Modelling

We investigate the properties of the components of this system via grid-fitting analysis of the composite spectral energy distribution (SED) of TIC 470710327 and TIC 470710327′, using photometric data from Vizier (Table A4). In order to reduce the degeneracies present in composit SED modelling, we use the synthetic SEDs from a grid of MIST isochrone models (Choi et al. 2016). This enforces that all components of TIC 4707410327 and TIC 470710327′ are the same age and located at the same distance, and ensures realistic physical parameters for all components. The models are reddened according to |$F_{\text{red}} = 10^{-A_{\lambda }/2.5}F_{\text{mod}}$| and A_λ = E(B − V)R_λ, where E(B − V) is a free parameter, and R_λ is calculated according to Cardelli, Clayton & Mathis (1989). Each model is evaluated against the data using a χ² metric.

We fit a composite four-component SED to the observations, assuming one component for each member of the triple system, and one for the nearby 0.5 arcsec companion, TIC 470710327′. The light contributions derived in Section 3 are used to inform the model, such that the total light contribution of the triple system is fixed to 77 per cent in the V band, with the remaining 23 per cent originating from TIC 470710327′. Additionally, we enforce that the luminosity of the primary component of the close binary must be higher than that of the secondary, i.e. L_A > L_B. Furthermore, we adopt the spectroscopically derived mass range for the tertiary: |$M_{\rm C}\in 14{-}17\, {\rm M}_{\odot }$| and only accept solutions that satisfy the triple system mass-ratio of q_c = 0.776. This results in 1019 possible solutions, which we rank according to their χ².

The solutions favour masses of the OB star of |$M_{\rm C}\in 14.5{-}16\, {\rm M}_{\odot }$|⁠, whereas the solutions for the components of the close binary orbit favour a mass ratio q₁ ∈ 0.9−1.0 with masses for the primary components ranging from |$M_A\in 6{-}7\, {\rm M}_{\odot }$| and |$M_B\in 5.5{-}6.3\, {\rm M}_{\odot }$| We find solutions for TIC 470710327′ in the range of |$11.5{-}14\, {\rm M}_{\odot }$|⁠. We note that the values of the mass of TIC 470710327′ are correlated with the mass of TIC 470710327 C due to the imposed light contributions of these objects. The model and observed SED are shown in Fig. 7. We note that this modelling returns a distance in the range of 4–4.5 kpc, whereas the Gaia parallax returns a distance closer to 1 kpc. In order to arrive at a distance close to the Gaia estimate, the reddening would need to be around 2, as opposed to the value in the range of 0.40–0.44 that we find. Although these two distances do not agree, the complexity of this system could contribute to an erroneous astrometric solution, which will likely be improved upon in future Gaia data releases.

Figure 7.

Observed SED (black points) with composite and individual model contributions for the best-fitting composite model. The individual contributions are unreddened, whereas we plot both the reddened (maroon) and unreddened (blue) composite models.

Open in new tabDownload slide

6.3 Dynamical stability

Even though both the inner and outer orbits of the triple system are relatively compact, the separation between the two orbits (P₂/P₁ ≈ 50) implies that the current configuration is dynamically stable. The stability criterion of Mardling & Aarseth (1999) implies long-term dynamical stability for TIC 470710327 for P₂ ≳ 18 d.

Dynamically stable systems can remain intact for many dynamical time-scales, giving rise to the possibility of observing three-body dynamics. ZKL cycles (Kozai 1962; Lidov 1962), for example, would manifest themselves as cyclic changes in the eccentricity of the inner orbit and in the mutual inclination between the inner and outer orbit. However, classical ZKL resonance can only occur in triples with mutual inclinations between |$39{_{.}^{\circ}}2\,{\rm and}\,140{_{.}^{\circ}}8$|⁠, meaning that with the mutually inclined of TIC 470710327, of |$\approx 16{_{.}^{\circ}}8 _ { -1{_{.}^{\circ}}38} ^ { +4{_{.}^{\circ}}18 }$|⁠, this effect is unlikely to be significant for the future evolution of the system.

The eZKL mechanism is expected to be important for the evolution of the system for values of ϵ_oct ≳ 0.001−0.01, and under the condition that the mass ratio of the inner binary is less than unity. Given the values in Table 2, we find ϵ_oct = 0.001 for q ≈ 0.9 and ϵ_oct = 0.01 for q ≈ 0.3. Given the morphology of the eclipses in the TESS data (and considering the substantial diluting third light), initial modelling suggests that the mass ratio of this inner binary is close to unity. Thus, the octupole term is within the range of relevance for the dynamical evolution of this system. The expected time-scale of the eZLK cycles are thought to be of order τ ≈ 180 yr/|$\sqrt{\epsilon _{\rm oct}}$| (Antognini 2015), i.e. τ ≈1800–5700 yr. We further investigated these time-scales with simulations using the TrES triple star evolution code (Toonen et al. 2016), and found the time-scales to be consistent with our analytical estimates. However, when including tides and gravitational wave radiation, the maximum amplitude of the eccentricity of the inner binary remains on the order of e₁ ∼ 0.001. Detailed constraints on the time-scale and amplitudes of these cycles would require detailed modelling of the inner binary, which remains impossible without full radial-velocity characterization of the inner components.

6.4 Possible formation scenarios

Several theories exist pertaining to the origin of higher order multiple systems. The formation of multiple systems is dependent on fragmentation of the natal material during the formation process. Hierarchical collapse within molecular clouds eventually leads to the formation of dense stellar cores and clumps, whose collapse results in the formation of stars and clusters. Throughout this, equatorial discs are formed through the conservation of angular momentum, which in turn can become involved in the accretion process and the formation of secondary cores. However, there is still debate as to what scale of fragmentation is the main cause of observed massive multiple systems – the fragmentation of the prestellar core or fragmentation of circumstellar discs. In both cases, opacity has a large effect on the initial separations of the systems, which cannot be less than around 10 au due to the opacity limit of fragmentation (Bate 1998; Boss 1998). Therefore, it is assumed that systems at closer separations than 10 au must have migrated to their observed positions (Bate, Bonnell & Bromm 2002). Some studies present evidence of disc fragmentation creating higher order massive multiple systems during their embedded phases (e.g. Megeath, Wilson & Corbin 2005). Other effects such as dynamical interactions between other companions and discs (Eggleton & Kisseleva-Eggleton 2006) could also cause inner binaries in multiple systems to harden into close orbits. Recent modelling by Oliva & Kuiper (2020), for example, has shown how disc fragments in the discs of massive protostars form through hierarchical fragmentation along spiral arms and migrate to spectroscopic orbits.

TIC 470710327 presents an interesting puzzle in terms of its formation given our derived geometry of a close binary with a total mass lower than the tertiary star. More massive stars have shorter Kelvin–Helmholtz time-scales than those of lower masses, so one would expect that the more massive tertiary was the first star to form. However, if this was the case, it is likely that when this star reached the main sequence (before the inner binary) it would have disrupted the remaining natal material and therefore discontinued the central binary’s formation. If the central binary did form first, a more consistent interpretation could be that the inner binary formed through disc fragmentation and the dynamical effects of this binary on the disc could have created a large overdensity at large radii. Such overdensities have been shown to occur in circumbinary discs in works such as Price et al. (2018). This overdensity could have accreted significant mass, perhaps accelerated by the continuing dynamical effects of the inner binary. In order for this scenario to hold the mass of the disc must have been very large, as the disc fragmentation process would not convert all the disc material into the eventual stars, and the combined mass of the stars in the tertiary system is at least 29 solar masses. However, the largest protostellar discs detected around massive young stellar objects are of order |$\sim \! 10 \, {\rm M}_{\odot }$| (e.g. Johnston et al. 2020, Frost et al. 2021).

An alternative explanation for the formation of TIC 470710327 is that this system is a result of the fragmentation of the prestellar core as opposed to fragmentation of a disc (e.g. Krumholz & Bonnell 2009). This is supported by the fact that additional sources surround this triple system in the local region, in the form of 0.5- and ∼22-arcsec distant companions. With a spurious association probability of ∼ 1 × 10⁻⁵ between TIC 470710327 and the 5 arcsec companion, this closer source is expected to have formed from the same core collapse. While the more distant source at ∼22 arcsec may not be currently bound to the tertiary system, it could still have come from the same prestellar core. The core collapse scenario circumvents the mass problem described above for disc fragmentation. Should this have occurred, we can assume that the close-binary formed from one collapse event, the tertiary star from another and the distant source from yet another. The close-binary system could have been formed by disc fragmentation as described above, and its hardening into a close-orbit could have been facilitated by dynamical interactions between the 14–17 M _⊙ star and the distant sources from others. Dynamical effects within the collapsing core could also have led to the capture of the  14–17 M _⊙ star by the close-binary system, forming the tertiary we see today. Dedicated radiative hydrodynamical modelling could help disentangle whether core fragmentation, disc fragmentation or a combination of both resulted in the formation of TIC 470710327, whilst repeat observations of all the sources in the region could help distinguish orbits and determine which stars are bound. Furthermore, as there are currently no resolved Gaia parallaxes for TIC 470710327 and TIC 470710327′, we cannot concretely rule out the possibility that these two targets are just unassociated nearby objects on the sky.

6.5 Future evolution

The observed mass ratio, whereby the tertiary is more massive than the combined mass of the inner binary, implies that the remainder of this system’s evolution will be driven by the evolution of the tertiary as this star will be the first to evolve off the main sequence. Given the expected mass of the tertiary of 14–17 M_⊙, star C is expected to fill its Roche lobe at an age of ∼ 13 Myr, at which point it will start transferring mass towards the inner binary. This type of mass transfer, from an outer star to an inner binary, is expected to occur in |$\sim \! 1{{\ \rm per\ cent}}$| of all triple systems in the Massive Star Catalogue (de Vries, Portegies Zwart & Figueira 2014; Hamers, Glanz & Neunteufel 2021a).

The outcome of such a mass transfer phase is, inherently, a hydrodynamical problem (de Vries et al. 2014; Portegies Zwart & Leigh 2019). If the inner binary is compact enough such that the mass transfer stream intersects itself, a circumbinary disc may form. Leigh et al. (2020b) argue that such a scenario leads to preferential accretion to the lowest mass star of the inner binary and therefore favours evolution towards equal mass inner binary stars. The circularization radius of the mass transfer material (Frank, King & Raine 2002; Toonen et al. 2016) is around |$20{R}_{\odot }$|⁠. As the inner binary has an orbital separation of |$\sim \! 10{R}_{\odot }$|⁠, a disc may form, but it is not clear from simple analytical calculation whether that disc would be stable enough to allow for secular accretion. If the mass transfer stream intersects the trajectory of the inner two stars in their orbit around each other, friction may reduce the orbit (de Vries et al. 2014) to lead to a contact system and/or a merger. Such a merger remnant may be considered a blue straggler for two reasons: it would be rejuvenated due to the merger and due to the accretion from the tertiary star.

Assuming a merger does take place, the triple would reduce to a binary system (stage 3 in Fig. 8). After the mass transfer phase ends the merger remnant would have a mass of 12–18 M_⊙ depending on how efficiently the binary was able to accrete matter. Given the mass of the former donor star it will evolve to become a neutron star. With typical post-mass transfer periods of several hundreds of days (i.e. orbital velocities of around 100 km s^-1), typical natal kicks from the supernova explosion, with magnitudes of several hundreds of km s^-1, would unbind the newly formed binary into two single stars (Hobbs et al. 2005; Verbunt, Igoshev & Cator 2017; Igoshev 2020), reducing the multiplicity of the system once more (5a in Fig. 8). If the orientation of the supernova kick is such that the binary remains intact the binary would likely undergo an additional mass transfer phase when the merger remnant evolves off the main-sequence. Given the large mass ratio of such a binary, the mass transfer would lead to a common-envelope phase (Ivanova et al. 2013). Consequently, this could either lead to a merger and the formation of a Thorne–Żytkow object (whereby the neutron star is enclosed by the red giant star, 6a in Fig. 8; Thorne & Zytkow 1975; Podsiadlowski, Cannon & Rees 1995; Levesque et al. 2014; Tabernero et al. 2021) or, if the binary survives, experience an ultra-stripped supernova (Tauris et al. 2015) and end up as a double neutron star (6b in Fig. 8).

Figure 8.

Possible future evolution of TIC 470710327. Relative sizes of the stars and orbits are not to scale. For simplicity, all orbits are depicted as circular.

Open in new tabDownload slide

7 CONCLUSIONS

TIC 470710327 is a compact, hierarchical triple system consisting of a 1.10 d binary containing two B-type stars and an OB-type tertiary on a wide 52 d orbit. The system was initially identified in TESS data by citizen scientists taking part in the PHT project. Using publicly available TESS data and newly obtained Hermes data, we report on the dynamical modelling of the system to reveal a rare configuration wherein the tertiary object in the wide orbit is more massive than the combined mass of the inner binary (m_AB = 10.9−13.2 M_⊙, m_C = 14−17 M_⊙). This configuration poses several challenges to explain its formation. Given the current mass of the tertiary component, star C is expected to fill its Roche lobe at an age of ∼13 Myr, meaning that the system must necessarily be younger than this age. Considering that the main-sequence lifetime of a 5.5–7 M_⊙ star is between 40–80 Myr, all of the stars in this system are currently in the main-sequence phase of their evolution. Furthermore, given the values from the dynamical modelling and constrained SED modelling, the primary binary component (star A) will fill between 70 and 83 per cent of its Roche lobe before the tertiary evolves off the main sequence, while the secondary binary component (star B) will fill between 68 and 75 per cent of its Roche lobe during this time. This means that the binary will not undergo mass transfer before experiencing Roche lobe overflow from the tertiary.

Given the compact orbits and the unusually high mass of the tertiary object, we speculate that the future evolution of this system will minimally involve one episode of mass transfer as the massive tertiary evolves across the Hertzsprung gap. Based on its initial mass, the tertiary will likely end its life as a neutron star. Alternatively, depending on the rate and efficiency of the mass transfer to the inner binary, the tertiary could evolve into an intermediate mass stripped star (Götberg et al. 2020). Should the binary system remain bound after the expected supernova kick (or supernovae kicks), this system could result in a close double neutron star gravitational wave progenitor, or an exotic Thorne–Żytkow object. Detection of more systems similar to TIC 470710327 would provide constraints on potential progenitor systems to gravitational wave events.

In addition to the triple system we report on two nearby stars that significantly contribute to the TESS aperture, located at angular separations of 22 arcsec (LS I + 61 72) and 0.5 arcsec (TIC 470710327′). Given the field density of stars around TIC 470710327, determined using Gaia eDR3, and the magnitude contrast between the target and the 0.5 arcsec companion, we show that the spurious association probability between TIC 470710327 and TIC 470710327′ is ∼ 1 × 10⁻⁵. These two nearby stars may be the source of the two additional periodic signals (P₃ = 9.97 d and P₄ = 4.01 d) seen in the TESS light curve. Additional photometric observations of the nearby stars are needed in order to probe the origins of the P₃ and P₄ and to determine the true multiplicity of this complex system.

With further observational characterization, particularly aimed at characterizing the nature of the inner binary, this system stands to become an excellent target to scrutinize simulations of massive star formation and evolution. Future spectroscopic observations that are specifically aimed at detecting the RV variations of both components of the inner binary would allow us to place tighter constraints on the dynamics of the system. In particular, these observations would allow us to probe the eZKL mechanism. Finally, RV characterization and detailed eclipse modelling of the inner binary would precisely constrain the light contributions of all components of the triple system. With such constraints, derivations of the atmospheric properties, such as T_eff, log g, v_rot, and L_bol would allow us to further test whether the three stars are coeval.

ACKNOWLEDGEMENTS

We thank the referee for their helpful comments that improved the manuscript. We thank all of the citizen scientists who take part in the Planet Hunters TESS project and who enable the discovery of exciting planet and stellar systems in TESS data. NE also thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant number 1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. Furthermore, NE and SA acknowledge support from the UK Science and Technology Facilities Council (STFC) under grant codes ST/R505006/1 and consolidated grant no. ST/S000488. This work also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 865624). CJ has received funding from NOVA, the European Research Council under the European Union’s Horizon 2020 research and innovation programme (N°670519:MAMSIE), and from the Research Foundation Flanders under grant agreement G0A2917N (BlackGEM). ST acknowledge support from the Netherlands Research Council NWO (VENI 639.041.645 grants). PGB was supported by NAWI Graz. SJ acknowledges support from the FWO PhD fellowship under project 11E1721N. KZAC acknowledges support under grant P/308614 of the IAC, which is financed by funds transferred from the Spanish Ministry of Science, Innovation and Universities (MCIU).

Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Based on observations made with the Mercator Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Based on observations obtained with the HERMES spectrograph on the Mercator telescope, which is supported by the Research Foundation – Flanders (FWO), Belgium, the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genéve, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany.

Finally, NE and CJ wish to thank the Harry Potter franchise for providing us with the in-house nickname for this system of Fluffy, inspired by Hagrid's three-headed dog in Harry Potter and the Philosopher's Stone.

DATA AVAILABILITY

The TESS data used within this paper are hosted and made publicly available by the Mikulski Archive for Space Telescopes (MAST, http://archive.stsci.edu/tess/). The TESS data described here may be obtained from https://dx.doi.org/10.17909/t9-5z05-k040. Similarly, the PHT classifications made by the citizen scientists can be found on the Planet Hunters Analysis Database (PHAD, https://mast.stsci.edu/phad/), which is also hosted by MAST.

Original Hermes spectra, and the newly obtained photometric data are available upon request.

This work made use of Astropy, a community-developed core python package for Astronomy (Astropy Collaboration et al.2013, 2018), matplotlib (Hunter 2007), pandas (McKinney et al. 2010), NumPy (Walt, Colbert & Varoquaux 2011), astroquery (Ginsburg et al. 2019), sklearn (Pedregosa et al. 2011), exoplanet (Foreman-Mackey et al. 2021) and pymc3(Salvatier, Wiecki & Fonnesbeck 2016).

Footnotes

REFERENCES

APPENDIX A: ADDITIONAL TABLES

Summary of the RV observations obtained with the hermes spectrograph (Section 4), ETVs determined from the TESS data (Section 3.3 and inputs used for the SED analysis (Section 6.2).