Young and eccentric: the quadruple system HD 86588

High-resolution spectroscopy and speckle interferometry reveal the young star HD 86588 as a quadruple system with a 3-tier hierarchy. The 0.3"resolved binary A,B with an estimated period around 300 years contains the 8-year pair Aa,Abc (also potentially resolvable), where Ab,Ac is a double-lined binary with equal components, for which we compute the spectroscopic orbit. Despite the short period of 2.4058 day, the orbit of Ab,Ac is eccentric (e=0.086+-0.003). It has a large inclination, but there are no eclipses; only a 4.4 mmag light modulation apparently caused by star spots on the components of this binary is detected with Evryscope. Assuming a moderate extinction of A_V = 0.5 mag and a parallax of 5.2 mas, we find that the stars are on or close to the main sequence (age>10 Myr) and their masses are from 1 to 1.3 solar. We measure the strength of the Lithium line in the visual secondary B which, together with rotation, suggests that the system is younger than 150 Myr. This object is located behind the extension of the Chamaeleon I dark cloud (which explains extinction and interstellar Sodium absorption), but apparently does not belong to it. We propose a scenario where the inner orbit has recently acquired its high eccentricity through dynamical interaction with the outer two components; it is now undergoing rapid tidal circularization on a time scale of ~1 Myr. Alternatively, the eccentricity could be excited quasi-stationary by the outer component Aa.


Introduction
The object of this study is a young, chromospherically active star HD86588 (Table 1). Alcala et al. (1995) classified it as a weak-lined T Tau star because of the strong lithium line; the star is not known to be variable. Several authors (e.g., Frink et al. 1998) attribute it to the Chamaeleon I (Cha I) star-forming region. The Gaia DR2 distance, 118 pc, and the corresponding heliocentric velocity (U, V, W)=(+5.4, −3.6, −3.2) km s −1 (the U axis is directed away from the Galactic center) are indeed similar to those of the ChaI members, although HD86588 is located on the sky outside the boundaries of known molecular clouds. No infrared excess indicative of a debris disk was found by Spangler et al. (2001). There are no emission lines in the spectrum. Lopez Martí et al. (2013) attribute this star to the field by its kinematics. Covino et al. (1997) detected triple lines in the spectrum, but they have not followed to determine the spectroscopic orbit. Later, the same team found this star to be "single", based on three spectra (Guenther et al. 2007). Meanwhile, it was resolved in 1996 by Koehler (2001) into a 0 27 visual pair denoted as KOH86 in the WDS (Mason et al. 2001). This object clearly deserves further study.
We found that this system is quadruple, rather than triple. Its hierarchical structure is shown in Figure 1. The brightest and most massive star Aa is a rapid rotator, so its lines, blended with other components, were not previously detected in the spectra. The star Aa is orbited by the spectroscopic binary Ab,Ac. The visual secondary B is at the upper level of the three-tier hierarchy.

Observational Data
High-resolution optical spectra used here were taken with the 1.5 m telescope sited at the Cerro Tololo Inter-American Observatory (CTIO) in Chile and operated by the SMARTS Consortium, 6 using the CHIRON optical echelle spectrograph (Tokovinin et al. 2013). Observations were conducted from 2018 March to June. Most spectra are taken in the slicer mode with a resolution of R=80000 and a signal-to-noise ratio of 25-40. Thorium-Argon calibrations were recorded for each spectrum.
Radial velocities (RVs) were measured by cross-correlation of the reduced spectrum with the binary mask. Further details of this procedure can be found in Tokovinin (2016). Figure 2 shows the cross-correlation functions (CCFs) recorded on two successive nights. The strong central component B has a constant RV of 2 km s −1 , while the two "satellites" Ab and Ac move around it rapidly, indicating presence of the short-period subsystem. Covino et al. (1997) measured the RVs of three components at +2, +85, and −98 km s −1 , in agreement with these CCFs. The fine details in the CCFs are likely produced by spots on rapidly rotating stars in the close binary. Naturally, these stars are chromospherically active, explaining the X-ray detection.
In most resolved triple-lined multiple systems, the stationary component of the CCF corresponds to one of the resolved stars, while the rapidly moving lines are produced by the close pair identified with another star. However, the fluxes of the components B and Ab+Ac are comparable, while speckle photometry indicates that B is fainter than A by 2 mag. One notes in Figure 2 that the amplitude of the "satellites" is changing, while the central component is stable in both RV and amplitude. When the CCF is fitted by 3 Gaussians, the area of the components Ab and Ac depends on their RV in a systematic way: the satellites become stronger when they approach the center and weaker when they are widely separated. Given the lack of strong photometric variability, evidenced below in Section 3, this effect cannot be caused by variable circumbinary extinction. Instead, it is explained by blending of the satellites with the fourth broad CCF component Aa that can be guessed by looking at the dashed line in Figure 2.
It turns out that the broad central component of the CCF, Aa, has the largest area and corresponds to the brightest star in this quadruple system. The six CCFs with well-separated satellites were modeled by fitting four Gaussians, and average parameters of each component were computed. The remaining CCFs where Aa is heavily blended with the other components were fitted by four Gaussians while the amplitude, width, and RV of the broad component were fixed to the average values. Figure 3 illustrates these two cases. One can see in the right panel that blending with Aa increases the apparent depth of the satellites and shifts their positions when they are close to each other. In the four-component fits, the dependence of the area of components Ab and Ac on their position (i.e., on the RV) vanishes. The rms residuals of the four-component fits to the CCFs are typically about 0.001. Table 2 lists the mean amplitude a and the mean dispersion σ of the Gaussians fitted to the subset of CCFs with wellseparated dips. The product aσ is proportional to the area of each component. The rms scatter of the amplitude is about 0.001 and the rms scatter of σ is about 0.4 km s −1 . Considering the complex nature of the CCF profiles and the small     amplitudes of the dips, the formal errors of the fitted parameters should be less than their real errors, dominated by the remaining systematic distortions of the CCFs. For this reason, we do not provide formal errors of mean parameters in Table 2. Table 2 also lists the average RVs of the two constant dips Aa and B and estimates of the projected rotation velocities V i sin , computed from σ by the approximate formula given in Tokovinin (2016); as this formula is valid for solar-like stars with σ<12 km s −1 , our estimates of V i sin are crude, except for B.
The individual RVs derived by fitting CCF with four Gaussians are listed in Table 3. The RVs of the broad-lined component Aa are determined from the six CCFs where the satellites are widely separated, with large errors (rms scatter 2.5 km s −1 ); their mean value is 7.8 km s −1 . The errors for the other components are on the order of 1 km s −1 . The RVs of the central component B range from 0.9 to 3.3 km s −1 ; the rms scatter of 0.8 km s −1 is caused by blending with other components.

Spectroscopic Orbit
As the two satellites Ab and Ac are equal and move rapidly, it was difficult to tell which is which in any individual observation. One of us (F.W.) took three spectra during one night and two more spectra on the following night. This established the orbital period at about 2 days, helping to tag the components and to derive the orbit. The elements of the 2.4 days spectroscopic orbit are listed in Table 4, while Figure 4 gives the RV curve. Despite the short period, the spectroscopic orbit is significantly non-circular. This result, based on the de-blended RVs derived from the four-component fits, also holds when using the original RVs derived from the three-component fits that neglect Aa. The rms residuals to the orbit are 1.1 and 1.0 km s −1 for Ab and Ac, respectively (the residuals were 2.4 and 1.6 km s −1 before de-blending). The residuals are unusually large for CHIRON spectra. However, they are explained by the complex and variable nature of the CCF profiles, illustrated in Figure 2. If a circular orbit is enforced, the residuals increase to 5 km s −1 .
A star of one solar radius synchronized with the orbit would rotate at 20.9 km s −1 , similar to our crudely estimated V i sin . This fact, together with the large M i sin 3 , implies a highly inclined orbit of Ab,Ac.

The Lithium Line
Knowing the RVs of each component, we can extract their individual spectral features by modeling and subtracting the contributions of other components, as explained in Tokovinin (2016). Briefly, the depth of the lithium line in Ab, Ac, and B (three parameters, as we ignore here the line of Aa, too wide to be measurable) is found by the simultaneous linear fit to all spectra using the RVs and the preliminary values of the line width. Then, for one component, the lines of the two remaining components are subtracted from each spectrum, all spectra are shifted to the zero velocity, and co-added. Figure 5 shows the lithium line in three components derived by this method. It is clearly present in B, with an equivalent width (EW) of 40±5 mÅ (or ∼200 mÅ if corrected for dilution by other components that contribute 0.8 of the light). The wide lithium lines in the satellites Ab and Ac are detectable, but their EW is not well measured. The EW of 200 mÅ reported by Alcala et al. (1995) was derived from the low-resolution spectra where all components were blended.

Sodium Lines
Narrow and deep interstellar lines of sodium are present in the spectra ( Figure 6). The CCF with a binary mask containing only those two lines shows that the RV of these features is constant at +14.5 km s −1 , while the dispersion of the Gaussian fitted to the CCF is only 3.4 km s −1 . The narrow sodium lines are caused by material on the line of sight that is unrelated to the multiple system. The RV of the sodium absorption matches  Figure 4. Spectroscopic orbit of the close pair Ab,Ac (see Table 4). Squares and full line correspond to the component Ab, triangles and dashed line to Ac.
the RV of stars and gas in the ChaI dark cloud (Covino et al. 1997). In addition to the absorption, some (but not all) spectra have a weak emission component in the D-lines (see the right panel of Figure 6); the RV of the emission is close to zero. The narrow sodium emission could come from the active chromosphere of the component B or from the sky light pollution.

Variability
Photometric variability of HD86588 is expected for several reasons: eclipses in the close binary, ellipsoidal variation, star spots, variable circumstellar extinction. Photometry in the SDSS g′ band was provided by the Evryscope instrument sited at Cerro Tololo (Law et al. 2015). A total of 17640 measurements covering the period from JD 2457755 to 2458078 with a typical cadence of 2.2 minutes yield the mean magnitude of 8.503 mag (in uncalibrated instrumental system) and the rms scatter of 19.84 mmag. Therefore, to the first order, the star is not variable on timescales from minutes to one year. The close binary is not eclipsing. Variable circumstellar extinction is also excluded. Figure 7 shows the periodogram in the period range from 0.8 to 5 days (there is no significant variability at longer periods). It is computed by fitting sine and cosine terms at each trial period, subtracting the fit, and estimating the variance of the residuals. The strongest details are located near the frequency of 0.4 cycles per day; the orbital frequency of the close binary is 0.4155 day −1 . Subtraction of the largest sine term with the frequency of 0.403 day −1 reduces the rms from 19.84 to 19.60 mmag, hence the rms variability at this frequency is 3.07 mmag and the corresponding sine-wave amplitude is 4.4 mmag. The micro-variability is apparently caused by star spots in the components Ab and Ac. They are rotating with periods that are close but not exactly equal to the binary period.
Components of slightly eccentric binaries usually rotate pseudo-synchronously with a frequency 1+6e 2 times the orbital frequency (see Equation (43) in Hut 1981), or 0.4338 day −1 for the Ab,Ac binary. The two strongest details in Figure 7 have frequencies of 0.403 and 0.429 day −1 . If they correspond to the rotation periods of the two components, one of them can be almost pseudo-synchronized, while the other rotates a bit slower. The tidal pseudo-synchronization timescale is orders of magnitude faster than the tidal circularization timescale (see Section 6), so it is consistent with a pseudosynchronized binary still having a slightly eccentric orbit. In such a case, the X-ray emission would be driven by tidal spin up as the system is currently experiencing rapid tidal circularization.
The variability at the frequency of 0.6 day −1 , also seen in Figure 7, could be caused by star spots on another rapidly rotating star in this system, either Aa or B.
Assuming the mass sum of 2   for Ab,Ac, the semimajor axis of the close binary is 0.044 au or 9.5 R e . The absence of eclipses restricts the inclination of the spectroscopic binary to i Ab,Ac <78°, or i sin 0.935 3 < . The minimum masses of the stars Ab and Ac are, therefore, 1.06   .

Positional Measurements
Apart from the discovery measure in 1996 by Koehler (2001), observations of the relative motion of the outer binary A,B (where A refers to the photo-center of the unresolved inner system) come from two sources. Vogt et al. (2012) monitored this pair with adaptive optics at VLT from 2007 to 2011. Furthermore, the relative position was measured by speckle interferometry at the Southern Astrophysical Research Telescope (SOAR) from 2011 to 2018 (see Tokovinin et al. 2018, and references therein). Most speckle measures at SOAR were made in the I band and result in the mean magnitude difference ΔI AB =1.80 mag, with the rms scatter of 0.12 mag. Additionally, we measured ΔV AB =1.90 mag. The totality of measures provides a relatively dense coverage from 2007 to 2018.
The pair A,B is in slow retrograde motion. It turned by 36°in 20 years since its discovery in 1996. This motion corresponds to an orbital period of a few hundred years. The period estimated from the projected separation is of the same order of magnitude. However, as shown in Figure 8, the observed motion is "wavy", suggesting presence of unresolved subsystem with a period of ∼8 years. This additional system presumably corresponds to the Aa,Abc pair, as shown in the diagram in Figure 1. Table 5 gives two sets of orbital elements describing the observed motion of A,B. The outer orbit is not constrained by the short observed arc. It is chosen to represent the data and to match the mass sum and parallax from the system model developed in Section 5. The circular inner orbit in Table 5 is sufficient to model the wobble, but the actual inner orbit can be eccentric just as well. These elements fit the data, but are by no means unique; they are one possible solution. The rms residuals to this model are 1.3 mas in both coordinates; they increase to 4 mas if the wobble is ignored. The measurements and residuals to the tentative orbits are given in Table 6. Note that the systems A,B and Aa,Abc have opposite directions of orbital motion, meaning that their orbits are definitely not coplanar.
The period of 8 years, mass sum of 3.4  , and parallax of 5.2 mas (see Section 5) correspond to the semimajor axis of 6 au or 31 mas on the sky. The wobble amplitude of 6.6 mas is only 21% of the semimajor axis, being reduced by the comparable fluxes and masses of Aa and Abc. The subsystem Aa,Abc should be resolvable by speckle interferometry. It was definitely unresolved at SOAR in 2018.0 and 2018.4, while all other observations do not reach the full diffraction-limited resolution of SOAR, 40 mas.
The center-of-mass velocity of Ab,Ac and the velocity of Aa should vary in anti-phase with a period of ∼8 years and an amplitude of the order of ∼10 km s −1 . At present, the RV difference between Aa and Abc (the center of mass of the inner binary) is 7.8-1.3=6.5 km s −1 . The RVs of Ab and Ac measured by Covino et al. (1997) imply that Abc had an RV of (85-98)/2=−6.5 km s −1 , different from its actual value of 1.3 km s −1 . Hence, the 8-yr orbit does produce the expected RV signature. Guenther et al. (2007)    −3.1) mas yr −1 , differs from the long-term Tycho-2 PM reported in Table 1. If the latter corresponds to the center-ofmass motion, the differential PM in 2015.5 (the DR2 epoch) was 10.5 masyr −1 , directed at 133°angle. The orbit of Aa,Abc in Table 5 predicts motion of the photo-center between 2015 and 2016 with a speed of 11.7 masyr −1 directed at 77°. The PM difference between Gaia and Tycho-2 roughly matches the proposed orbit in speed but not in direction. We predict that Gaia will soon detect astrometric acceleration and, possibly, will derive the astrometric orbit of Aa,Abc.

Modeling
We make the first, preliminary attempt at modeling the system by evaluating individual magnitudes of the components and comparing them to the isochrones. The photometry by Vogt et al. (2012) yields a good measurement of Δm AB = 1.50±0.12 mag at 2.18 μm wavelength that can be identified with ΔK AB . They also measured Δm AB =1.55 mag at 1.265 μm. Speckle interferometry at SOAR gives ΔI AB = 1.8 mag and ΔV AB =1.9 mag with errors of ∼0.1 mag. This, together with the total brightness of the system, defines the individual V and K magnitudes of the resolved components A and B.
The areas of the CCF dips, corrected slightly for the temperature dependence by using preliminary assumed effective temperatures, lead to the flux fractions of 0.54, 0.16, 0.16, and 0.14 for Aa, Ab, Ac, and B, respectively, in the V band. The corresponding ΔV AB =1.9 mag matches the speckle photometry. The relative fluxes of the unresolved components Aa, Ab, and Ac in the K band are not measured. They are computed by assuming ΔK ≈0.65ΔV, as derived from standard relations for main-sequence stars of these masses; the measured V-band flux ratios are re-scaled to the K band. This defines the individual magnitudes of all components listed in the first two lines of Table 7.
We begin by comparing these stars to the 1 Gyr solarmetallicity Dartmouth isochrone (Dotter et al. 2008). The Gaia DR2 parallax of 8.50±0.77 mas places the components slightly below the main sequence. The corresponding masses of Ab and Ac are 0.85   , contradicting their minimum mass of 1.06   derived from the orbit. However, the DR2 astrometry of this star has unusually large errors, likely caused by the duplicity and by the acceleration in the 8-yr orbit. Poor quality of the astrometric solution for HD86588 is manifested by the excess noise of 5.2 mas (its typical value in DR2 is 0.2 mas). Therefore, the DR2 parallax can be biased.
The total extinction in the direction of HD86588 is A V =0.48 mag according to Schlafly & Finkbeiner (2011). 7 The agreement between model and observations is improved if we adopt A V =0.5 mag and a smaller parallax of 5.2 mas. This makes the stars brighter and more massive (circles in Figure 9).
The main-sequence masses of all components, also given in Table 7, then roughly match the minimum masses of Ab and Ac. The extinction is confirmed by the presence of interstellar lines in Figure 6. The RV of these lines indicates that the dust is related to the ChaI dark cloud, while the RV of HD86588 is different by ∼10 km s −1 .
We stress that the proposed system model is based on several assumptions and is by no means definitive. The V−K color of Ab and Ac is not measured directly, hence their location on the main sequence is artificial. The measured color of the component B places it above the main sequence, under the 25 Myr isochrone.
The age of the system can be inferred from the strength of the lithium 6708 Å line. The component B contributes about 0.2 fraction of the total flux at this wavelength, so the measured EW of 40 mÅ translates to the intrinsic EW of 0.2 Å. The effective temperature inferred from the isochrone is 5920 K ( T log 3.77 e = , spectral type G0V). According to Figure 6 of Covino et al. (1997), the star is located above the upper envelope for the Pleiades, somewhere among massive members of ChaI. Apparently, the system is younger than the Pleiades (i.e., younger than ∼150 Myr), but we cannot say how much younger.
Stellar rotation provides another diagnostic of age. The projected rotation velocity of the solar-mass component B, 12.2 km s −1 , implies its rotation period P rot,B <4.1 days. The fast rotation of Aa, estimated very crudely from its CCF, corresponds to P rot,Aa <1 day. Adopting the B−V colors of 0.63 and 0.43 mag for B and Aa, respectively, derived from the isochrones, the age-rotation relation of Barnes (2007) leads to the minimum ages of 150 and 100 Myr. If the periodogram detail at 0.6 day −1 corresponds to the true rotation period of B, the same relation gives the age of 25 Myr.

Origin of the Eccentric Inner Binary
Nonzero eccentricity of the inner binary Ab,Ac provides an estimate of its "tidal" age, i.e., the time of tidal orbit circularization. Using the equilibrium-tide model with convective damping and a scale factor of F tid,con =20, as  parameterized in Belczynski et al. (2008), a pair of stars with M M 1 1 2  = =  , R 1 =R 2 =1 R e and a period of 2.4 days evolves from e=0.3 to e=0.1 in 0.8 Myr. At the pre-mainsequence (PMS) stage, when the stellar radii are larger, the tidal evolution is much faster; PMS binaries with periods shorter than ∼6 days are already circularized (Meibom & Mathieu 2005). So, the inner binary in HD86588 has a tidal age of ∼1Myr, much less than the actual system age. Somehow, the binary had to become eccentric after its components reached the main sequence.
As suggested by Moe & Kratter (2018), most eccentric binaries with periods below the circularization period acquired their eccentricity relatively recently through dynamical interaction in hierarchical systems (Kozai-Lidov, or K-L, cycles). Therefore, their "tidal age" is a small fraction of the true age. However, compact hierarchies with outer periods as short as 8 years tend to have nearly coplanar orbits, preventing the K-L cycles. Moreover, even if the orbits were inclined, the period of K-L cycles and the associated timescale of orbit evolution are much shorter than the system age, leaving the paradox of young tidal age unsolved.
The fact that this system is quadruple, rather than triple, comes to the rescue. Hamers et al. (2015) studied dynamics of quadruple systems with 3+1 hierarchy composed of "nested" binaries X, Y, and Z (in our case, the 2.4 days, 8 years, and 300 years systems). They found that the inner binary X can acquire a large eccentricity even if it were originally coplanar with the intermediate binary Y. This happens when the periods of K-L cycles in the inner triple X-Y and outer triple Y-Z are comparable. The K-L period for the inner triple is given by their Equation (11) The K-L period in the outer triple Y-Z is computed in the same way. Using the parameters determined above, we estimate for HD86588 P KL,X−Y ≈26 kyr and P KL,Y−Z ≈47 kyr. So, these periods are comparable and the dynamical evolution of the inner binary to large eccentricity by this mechanism is possible, at least in principle.
An alternative mechanism, frequently invoked in the literature, is constant excitation of the inner binary's eccentricity by the outer companion (in our case Aa), accompanied by tidal damping. The small eccentricity is then defined by the balance of these opposite processes and persists for a long time. The long-lived eccentric inner binary makes its discovery more probable, compared to the chance of catching it during the fast tidal damping episode. If the eccentricity of the Aa,Abc orbit is e Y <0.8, then general relativistic (GR) and tidal precession currently suppresses such continuous excitation of the inner binary's eccentricity (Liu et al. 2015). However, this orbit likely oscillates between small and large eccentricities due to K-L cycles with component B. The K-L timescale P KL,Y−Z =47 kyr is also substantially shorter than the tidal circularization timescale of e e 1 Myṙ . If the 8-yr orbit can reach e Y >0.8 in its K-L cycle, then it can sustain the eccentricity e X =0.09 of the inner binary, counteracting tidal friction as well as GR and tidal precession.

Summary and Discussion
Originally, the interest in HD86588 was driven by its presumed PMS status related to its membership in ChaI. We found that this star is located behind the ChaI molecular cloud, which imprints sodium absorption in its spectrum and causes a mild extinction. The RV of HD86588 differs from the RV of the ChaI group by ∼10 km s −1 . However, this system is definitely juvenile (younger than ∼150 Myr), as evidenced by the presence of strong lithium line and fast rotation. It could be born in a previous episode of star formation related to ChaI.
The complex nature of this three-tier hierarchical system adds uncertainty to its interpretation. We cannot trust the Gaia parallax until the 8-yr wobble is accounted for in its astrometric solution. Our best guess is a parallax of 5.2 mas (distance 190 pc); assuming the extinction of A V =0.5 mag, we propose a model where the stars Aa, Ab, and Ac are located on the main sequence, and their masses match the minimum mass derived from the spectroscopic orbit. The adopted parallax and the system velocity of 3.5 km s −1 (estimated center-of-mass RV of Aa,Abc) lead to the revised heliocentric motion of (U, V, W)= (+8.7, −6.1, −5.3) km s −1 . This velocity does not match any known kinematic group.
The close binary Ab,Ac has a nonzero eccentricity e=0.09, unusual for its short period of 2.4 days. However, this binary is not unique. The multiple star catalog (Tokovinin 2018) contains many spectroscopic binaries belonging to hierarchical systems. Among 184 such pairs with periods less than 4 days and primary mass less than 1.5 solar, 30 have nonzero eccentricity. Eight of them have a complex hierarchy with three or four tiers, reminiscent of HD86588, while the rest are triple or 2+2 quadruple systems with just two tiers. An example of the second group is MSC 16228−2326, a PMS 2+2 quadruple system RX J1622.7−2325 where the pair Ba,Bb composed of two 0.4   stars has P=3.23 days and e=0.30 (Rosero et al. 2011). These authors show the period-eccentricity plot where the PMS and main-sequence binaries occupy the same locus. They argue that eccentric orbits with short periods are found in both groups, independently of age, and suggest that eccentricity can be excited by dynamical interaction with other components in multiple systems.
Although HD86588 is no longer interesting as a calibrator of PMS evolution, its further study can clarify mechanisms of close-binary formation within stellar hierarchies. The next obvious step would be a spatial resolution of the 8-yr subsystem, either by speckle interferometry at 8 m telescopes, or by long baseline interferometers such as VLTI. In parallel, spectroscopic monitoring during several years and the Gaia astrometry will complement interferometry in defining the intermediate orbit. This will provide accurate measurement of stellar masses for all components.