A Quintuple Star System Containing Two Eclipsing Binaries
We present a quintuple star system that contains two eclipsing binaries. The unusual architecture includes two stellar images separated by on the sky: EPIC 212651213 and EPIC 212651234. The more easterly image (212651213) actually hosts both eclipsing binaries which are resolved within that image at , while the westerly image (212651234) appears to be single in adaptive optics (AO), speckle imaging, and radial velocity (RV) studies. The ‘A’ binary is circular with a 5.1-day period, while the ‘B’ binary is eccentric with a 13.1-day period. The velocities of the A and B binaries are different by 10 km s. That, coupled with their resolved projected separation of , indicates that the orbital period and separation of the ‘C’ binary (consisting of A orbiting B) are years and AU, respectively, under the simplifying assumption of a circular orbit. Motion within the C orbit should be discernible via future RV, AO, and speckle imaging studies within a couple of years. The C system (i.e., 212651213) has a radial velocity and proper motion that differ from that of 212651234 by only 1.4 km s and 3 mas yr. This set of similar space velocities in 3 dimensions strongly implies that these two objects are also physically bound, making this at least a quintuple star system.
Subject headings:stars: binaries (including multiple): close—stars: binaries: eclipsing—stars: binaries: general—stars: binaries: spectroscopic—stars: binaries: visual
The Kepler mission, with its exquisite photometric precision (Borucki et al. 2010; Batalha et al. 2011) has revolutionized the study of stars in general, and binary stars, in particular. There were some 3000 binaries discovered in the Kepler main field (Slawson et al. 2011; Matijevič et al. 2012; Kirk et al. 2016), and there is a growing collection of binaries that have been found to date in the 2-wheel extension of the Kepler mission, called ‘K2’ (Howell et al. 2014; LaCourse et al. 2015). Among this impressive collection of mostly eclipsing binaries, some 220 triple stars have been discovered (Rappaport et al. 2013; Conroy et al. 2014; Borkovits et al. 2015; Borkovits et al. 2016), mostly through eclipse timing variations, but some via 3rd-body eclipses. In addition to the large sample of triple-star systems, a number of higher-order multiple star systems have also been discovered using Kepler data plus follow-up ground-based observations (quadruple 4247791: Lehmann et al. 2012; quadruple 7177553: Lehmann et al. 2016; and quintuple KIC 4150611/HD 181469: Shibahashi & Kurtz 2012, and references therein; Prsa et al. 2016). Other quintuple systems, not found with Kepler, include: 1SWASP J093010.78+533859.5 (Lohr et al. 2015); the young B-star quintuple HD 27638 (Torres 2006); HD 155448 (Schütz et al. 2011); 14 Aurigae (Barstow et al. 2001); Coronae Borealis (Raghavan et al. 2009); GG Tau (Di Folco et al. 2014), HIP 28790/28764 and HIP 64478 (Tokovinin 2016), and V994 Her (Zasche & Uhlař 2013). Aside from being intrinsically fascinating systems, we can learn the most about the formation and evolution of multiple star systems if the constituent stars are close enough to interact dynamically on timescales of less than a few years.
In this work we report the first quintuple star system found in the K2 fields, and one of a relatively few that contain two eclipsing and spectroscopic binaries. This work is organized as follows. In Sect. 2 we present the K2 Field 6 data in which the two eclipsing binaries, ‘A’ and ‘B’, were discovered to be in the same stellar image. Our optical radial velocity studies are discussed in Sect. 3, and Keplerian orbits are fit to the RV data in Sect. 4. We evaluate the full set of parameters for the ‘A’ and ‘B’ binaries in Sect. 5 using a Monte Carlo evaluation process. The spectra are analyzed in Sect. 6 to further confirm our parameter determinations for binaries ‘A’ and ‘B’. We consider the tidal status of binaries ‘A’ and ‘B’ in Sect. 7. In an effort to further probe the structure of the system we obtained adaptive-optics and speckle-interferometric images of the system; these are presented in Sect. 8. The photometric distance to the target is derived in Sect. 9. We discuss in Sect. 10 the constraints that can be set on ‘binary C’ which is composed of binaries ‘A’ and ‘B’ that are gravitationally bound to each other. Binary ‘D’ comprised of binary ‘C’ in orbit about the apparently single star EPIC 212651234, is discussed in Sect. 11. We investigate the single star EPIC 212651234 in Sect. 12. Finally, we summarize our results in Sect. 13.
2. K2 Observations
As part of our ongoing search for eclipsing binaries, we downloaded all available K2 Extracted Lightcurves common to Campaign Field 6 from the MAST111http://archive.stsci.edu/k2/data_search/search.php. The flux data from all 28,000 targets were high-pass filtered with a cut-on frequency of 1 day. The data files were then Fourier transformed to facilitate a search for objects with periodic signatures. The folded lightcurves of targets with significant peaks in their FFTs were then examined by eye to look for unusual objects with periodic features. Unfolded and normalized light curves for targets lacking significant peaks in their FFTs were also examined by eye to look for aperiodic or quasi-periodic signatures.
One object that stood out was EPIC 212651213 (hereafter ‘E1213’) which exhibited eclipses with a 5-day period, but where the folded light curve also showed the presence of eclipses that did not fit this period. The raw light curve is shown in Fig. 1 and the two nearly equal eclipses of the 5-day orbital period are readily apparent. However, there is also another set of unequally spaced, more shallow eclipses whose locations are marked in the figure with short vertical red lines to guide the eye. These correspond to a different eclipsing eccentric binary with a 13-day period.
We also noticed that the target EPIC 212651234 (hereafter ‘E1234’) exhibited the same two sets of eclipses, but with reduced amplitudes. The SDSS image of the sky in the vicinity of these two targets (see Fig. 2) shows two comparably bright stellar images that are separated by 11. The properties of these two stars are summarized in Table LABEL:tbl:mags. Investigation of the Target Pixel Files for E1213 and E1234 with Guest Observer software PyKE222http://keplerscience.arc.nasa.gov/software.html#pyke (Still & Barclay 2012)333http://adsabs.harvard.edu/abs/2012ascl.soft08004S confirmed that the location of the photometric aperture for each of these two targets slightly overlaps the other star. By narrowing these apertures and performing a series of independent extractions on each target we were able to tentatively conclude that the source of both eclipsing binary signals is, in fact, E1213. (We have subsequently confirmed this with narrow-slit and fiber-fed spectroscopic observations – see Sect. 3.)
We present the folded lightcurves for each of the eclipsing binaries in Fig. 3. We used the Ames K2 PDC data, which have been corrected for the dilution by the other star. In each case, we subtracted the folded light curve of the other eclipsing binary in such a way as to effectively ‘clean’ the data of the modulations due to the other binary. To accomplish this, the folded data used for the subtraction was produced while eliminating the phases around the other binary’s eclipses.
The eclipse timing analysis, for determining orbital periods, phase zero of the eclipses, and relative eclipse separations, was carried out as described in Borkovits et al. (2015) and Borkovits et al. (2016). This includes allowance for local trends in the out-of-eclipse regions near the eclipses.
The 5.077-day binary (top panel of Fig. 3), hereafter designated as ‘Binary A’, has two nearly equal depth eclipses that are separated by very close to half an orbital cycle. The fitted fractional separation between eclipses differs by only ppm from half an orbital cycle, and we can then utilize the approximate expression (good to 2nd order in eccentricity ):
to say that , which is indicative of a rather circular orbit unless the longitude of periastron, , for the binary is within of or . In this expression is the fractional time difference between the two closely spaced eclipses. We can also utilize information from the relative widths of the two eclipses, and , to find a measure of . For small and arbitrary :
Finally, we note the out-of-eclipse modulation with a 1% amplitude that is in phase with the binary orbit, but is not due to ellipsoidal light variations (which would have minima at the eclipses), and we ascribe this modulation to star-spots on one or both stars of the “A” binary. Figure 3 also shows that the apparently irregular variations seen in the out-of-eclipse regions in Fig. 1 could be well modeled in this way.
In the bottom panel of Fig. 3 the 13.1947-day eccentric binary is more clearly revealed. For this set of eclipses we find the closely spaced pair to be separated by 0.4019 () orbital cycles. From Eqn. 1 we find . In Sect. 4 we are able to compare this to the results we find by directly measuring and from RV studies, and find excellent agreement.
Notes. (a) Taken from the SDSS image (Ahn et al. 2012). (b) From APASS (https://www.aavso.org/apass). (c) Carlsberg Meridian Catalog (http://svo2.cab.inta-csic.es/vocats/cmc15/). (d) 2MASS catalog (Skrutskie et al. 2006). (e) WISE point source catalog (Cutri et al. 2013). (f) Based on photometric parallax only. (g) The mass-weighted mean of the velocities for the ‘A’ and ‘B’ binaries (see text). (h) From UCAC4 (Zacharias et al. 2013), and Huber et al. (2015).
3. Radial Velocity Observations
We have carried out 27 radial velocity measurements of E1213 on 25 different nights. Eight of these were carried out with the TRES spectrograph at the Fred Lawrence Whipple Observatory, 16 were obtained at Thüringer Landessternwarte Tautenburg (TLS), and 3 were acquired at the Konkoly Observatory of the Hungarian Academy of Sciences. Here we briefly describe the spectroscopic observations at each facility. In addition, we have obtained 11 spectra of E1234.
3.1. RV observations with the TLS Spectrograph
Spectra of E1213 and E1234 were taken using the Coudé-Echelle Spectrograph on the 2-m Alfred Jensch Telescope at TLS. The spectrograph was used with a projected slit width of 2 arcsec providing a resolving power of R = 30,000. We obtained 16 spectra of E1213 and three spectra of E1234 with 40 minutes of exposure time each.
Spectrum reduction was done using standard ESO-MIDAS packages and included bias and stray-light subtraction, filtering of cosmic-ray events, flat fielding using a halogen lamp, optimum order subtraction, wavelength calibration using a ThAr lamp, normalization to the local continuum, and merging of spectral orders. Small shifts in the instrumental zero point were removed from each spectrum using a large number of telluric O lines as reference. The long-time RV accuracy reached with this spectrograph and reduction procedure depends on spectral type and of the observed star and is about 100 ms for sharp-lined, solar-type SB1 stars.
3.2. RV observations with the TRES spectrograph
Spectra of E1213 and E1234 were taken using the Tillinghast Reflector Echelle Spectrograph (TRES; Szentgyorgyi & Furész 2007), on the 1.5 m telescope at the Fred Lawrence Whipple Observatory (FLWO) on Mt. Hopkins, Arizona. TRES is a fiber-fed optical spectrograph with a resolving power of R = 44,000. We obtained 8 spectra of E1213 between UT 2016-02-21 and UT 2016-03-23 with a typical exposure time between 660-1200 seconds, and an average signal-to-noise ratio (S/N) per resolution element of 38. E1234 was observed 6 times between UT 2016-02-21 and UT 2016-05-03 with the exposure time ranging between 360-1500 seconds, and an average S/N of 28. The spectra were extracted following the procedures reported by Buchhave et al. (2010).
3.3. RV observations at the Konkoly Observatory
E1213 and E1234 were observed during three nights in March 2016 with the fiber-fed ACE spectrograph attached to the 1-m RCC telescope of the Konkoly Observatory at Piszkés-tető, Hungary. The instrument covers the 0.415 m – 0.915 m wavelength range with a resolution of R = 20,000. The fiber entrance has a sky-projected diameter of , which assured complete separation of the light from the two stellar images of E1213 and E1234. The light of a Thorium–Argon (ThAr) lamp can be projected on the entrance of the same fiber that is used for stellar observations for precise wavelength calibration.
The individual 1800-s exposures had S/N 25 per pixel around 0.55 m. The RV measurements were averaged for 3 – 4 consecutive exposures, with the exception of the night 2 457 466, when only one exposure of E1213 was obtained.
The spectra were reduced using IRAF444IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. standard tasks including bias, dark, and flat-field corrections, aperture extraction, wavelength calibration (using ThAr exposures) and barycentric correction. The normalisation, cosmic-ray filtering, order merging and cross-correlation were performed by in-house developed programs of ÁS.
The systematic errors introduced by the data processing and the stability of the wavelength calibration system was found to be better than 0.36 km s, based on RV standard observations (Derekas et al. 2016).
|BJD-2457000||Star 3||Star 1||Star 2||Observatory|
Notes. (a) Units are km s and the typical empirical uncertainties on the RV values are 1.7 km s for star 3 and 7 km s for stars 1 and 2. (b) ‘TLS’ stands for the Thüringer Landessternwarte Tautenburg, ‘KO’, the Konkoly Observatory, and ‘TRES” is the Tillinghast Reflector Echelle Spectrograph of the Whipple Observatory.
Notes. (a) Units are km s and the typical uncertainties on the RV values are 0.1 km s for any one observatory. (b) ‘TLS’ stands for Thüringer Landessternwarte Tautenburg, ‘KO’, the Konkoly Observatory, and ‘TRES” is the Tillinghast Reflector Echelle Spectrograph of the Whipple Observatory.
3.4. RV Analysis
For each observation we cross correlated the acquired spectrum with a single reference template stellar spectrum. We utilized the template spectrum that yields the highest peak correlation value after testing a grid of and values, for an assumed solar metallicity. An illustrative cross-correlation function (‘CCF’) is shown in Fig. 4. The reference spectrum was for K and , close to values ultimately derived for star 3. The three clearly discernible peaks correspond to three stars in the two binaries that are detectable in these RV studies. The RVs were measured by modeling the CCFs using two profiles with common centers, one Gaussian plus one Cauchy, per stellar component. When the peaks are clearly separable, as in the example in Fig. 4, the RVs of these three stars can be well measured. When the peaks are closer, or even merged, we may be able to determine only the RV of star 3 which dominates the CCF in all observations. Lines of star 4 could not be detected in any of the spectra. In particular, we checked when stars 1 and 2 have very small RVs, and simultaneously the RV of star 3 is large—star 4 is still not detected in the CCFs, even with a lower temperature reference spectrum. The measured RVs for the three detectable stars in E1213 are summarized in Table 2.
We also obtained 11 spectra for E1234. Only a single peak in the CCF is measured in any of these. The 11 RVs for this image are reported in Table 3. We interpret this in terms of either a single star within E1234 or perhaps a binary with an orbit that is sufficiently wide that no RV changes are detected. We can set a limit of km s for RV variations over the 72-day interval of the observations (for more stringent constraints, see Sect. 12).
|Parameter||Star 1||Star 2|
Notes. (a) All values are cited with 1- uncertainties; times are BJD-2457000. (b) Based on the K2 photometry. (c) From fits to the radial velocity data given in Table 2, with eccentricity fixed at . (d) Mass function for each star based on their respective velocities. (e) From an analysis using the Phoebe binary light curve emulator (Prša & Zwitter 2005). (f) From a Monte Carlo analysis of the system parameters and error propagation (see Sect. 5). (g) The extra constraints needed to yield the full set of system parameters were obtained from observational restrictions on the light ratio between the ‘A’ and ‘B’ binaries, and the Yonsei-Yale evolution tracks (see Sect. 5).
4. Radial Velocity Orbits
4.1. Fitting the ‘A’ binary
Given the very small value of for the ‘A’ binary that we were able to deduce from the K2 light curve (see Sect. 2), we restricted our fits to circular orbits. We fit the RV points shown in the upper panel of Fig. 5 simultaneously for the orbits of stars 1 and 2 using common values of the orbital phase and velocity, as well as and velocities. We fixed the orbital period at d, as determined from the K2 data. This is more accurate than the period we would be able to determine from the RV data alone. We utilized a simple Levenberg-Marquardt fitting routine to determine the four unknown parameters (, , , ) and their uncertainties.
The best-fit model RV curves for the ‘A’ binary are shown in the top panel in Fig. 5, and the fitted orbital parameters are summarized in Table 4. The values of and are 73.5 km s and 78.1 km s, respectively; these are similar but do differ at the 1.6- level. The orbital phase agrees with the phase obtained from the K2 photometry and projected forward by nearly half a year. The agreement is good to within 0.047 days with an uncertainty of 0.025 days (i.e., hr).
We also note that star 1, with the somewhat lower velocity (compared to star 2), and hence the higher mass, is the one whose descending node on the RV curve is an integer number of orbital cycles from the slightly deeper eclipse in the K2 light curve (see the top panel of Fig. 3). This is the time when the more massive star is eclipsed. And, at least for stars on the main sequence, the more massive star should have the higher surface brightness, and hence the deeper eclipse when it is transited—just as we observe.
We use the values to derive the constituent masses in Sect. 5.
Notes. (a) All values are cited with 1- uncertainties; times are BJD-2457000. (b) Based on the K2 photometry. (c) From fits to the radial velocity data given in Table 2. (d) Mass function based on the velocity and the orbital eccentricity. (e) From an analysis using the Phoebe binary light curve emulator (Prša & Zwitter 2005). (f) The extra constraints needed to yield the full set of system parameters were obtained from observational restrictions on the light ratio between the ‘A’ and ‘B’ binaries, and the Yonsei-Yale evolution tracks (see Sect. 5).
4.2. Fitting the ‘B’ binary
We next fit the RV points plotted in the bottom panel of Fig. 5 for the obviously eccentric orbit. We held the orbital period fixed at 13.1947 days, as determined from the K2 photometry, and fit for 5 free parameters: , where and are the longitude and passage time of periastron. For this we utilized a Markov chain Monte Carlo (‘MCMC’) code to map out the uncertainties in the parameters. The results are reported in Table 5. In spite of the fact that this is a single-line spectroscopic binary, we were still able to make useful mass estimates of both star 3 and star 4 via their combined flux compared to the flux from binary “A” (see Sect. 5).
5. Evaluation of Basic Stellar and Atmospheric Parameters
We now have two independent mass functions for star 1 and star 2 in the ‘A’ binary, and a mass function for star 3 in the ‘B’ binary. For both binaries, the orbital inclination angles are determined with the Phoebe binary light curve emulator (Prša & Zwitter 2005), and are very close to 86. To determine the masses and other stellar properties we utilize two ingredients: (i) the three mass functions and their uncertainties (Tables 4 and 5), and (ii) the relative fluxes of the “A” and “B” binaries (see, e.g., Sect. 8.2). For the latter constraint we adopt the following specific equality involving the stellar luminosities:
We then use a Monte Carlo error propagation technique to evaluate all four masses simultaneously. First we choose realizations for stars 1 and 2, using the cited values for the orbital velocities and random Gaussian draws of their cited uncertainties (Tables 4 and 5). This then determines and for that particular realization of and . Next we similarly chose a random realization of the value of using its error bar. We then choose a random mass ratio for between 0 and 1; that, plus the mass function for star 3, yields a specific realization for and .
In order to further constrain the system parameters, we require at least one more substantive constraint on the system. We take this supplemental information to be the measured ratio of the brightness of binary ‘B’ to binary ‘A’ (see Sect. 8). However, in order to incorporate the luminosities of the star, we need to make use of basic evolution tracks for low-mass stars. In order to do that, we have devised non-physical fitting functions to the Yonsei-Yale evolution tracks (see Fig. 6; Yi et al. 2001). The following represent the radius and luminosity evolution of a star of mass and evolution age, .
where is in solar masses, and the evolution age, , is expressed in Gyr. These should be sufficiently accurate for stars of mass .
Thus, our Monte Carlo realizations of the system also require an evolutionary age, which we take to be uniformly distributed between 1 and 10 Gyr. Once we have a set of 4 masses and a randomly chosen evolution age, we can use Eqns. (4) and (5) to evaluate all the system parameters, including masses, radii, luminosities, and values of g. The ratio of is then computed for the given realization. This ratio is then checked against the measured ratio given in Eqn. (3) with a randomly chosen value with respect to the uncertainty given in that equation. If the computed ratio lies outside that of the ‘measurement’, the entire system is discarded.
This process is repeated times, and distributions for the four masses, radii, , and g’s are built up. The results are shown in Fig. 7. A perusal of these mass distributions reveals a few interesting things. For the ‘A’ binary, the masses of stars 1 and 2 are: and , respectively. The most probable values of g are 4.54 and 4.59, with an uncertainty of 0.1 dex, respectively. The locations of stars 1 and 2 in the plane are shown in Fig. 6, and they are clearly in good accord with the corresponding mass of the tracks they are near.
As for the ‘B’ binary, the mass distributions for star 3 and star 4 (see Fig. 7) indicate and . These are rather well defined given that this is a single-line spectroscopic binary, and result from the luminosity comparison of binary “B” with binary “A”. The most probable values for g are and , respectively. The locations of these stars comprising binary “B” are shown in the plane in Fig. 6. Star 3 appears to be slightly evolved.
6. Spectral analysis
The TLS spectra were taken during grey time and some of them are slightly contaminated with moonlight. Though this contamination could be separated from the stellar features in most of the CCFs used for the RV determination, we had to be sure that the contamination from the moonlight does not enter the spectrum analysis. We therefore initially used only nine of the 16 TLS spectra of E1213 listed in Table 2 for the analysis. This number was too small for a secure decomposition of the observed composite spectra into the spectra of the individual stellar components. Nonetheless, we tried the FFT-based KOREL program (Hadrava 1995; 2004) that delivers the orbital parameters together with the RVs and the decomposed spectra, and then the singular value decomposition of Simon & Sturm (1994) with the measured RVs as input values using the CRES code (Ilijic 2004). However, in both cases the subsequent analysis of the decomposed spectra gave inconsistent results. In the end, we decided to analyze only the one composite spectrum that shows the largest separation of the components’ peaks in the CCF (Fig. 4).
Spectrum analysis was performed using the GSSP program (see Lehmann et al. 2011 and Tkachenko et al. 2012 for a description of the method). The program is based on the spectrum synthesis approach. We extended the method to a simultaneous fit for three components in a composite spectrum. The synthetic spectra were computed on the TLS cluster computer with the parallelized version of the SynthV program (Tsymbal 1996), using a library of precomputed line-by-line model atmospheres (Shulyak et al. 2004). The atomic data were taken from the VALD data base (Kupka et al. 1999). We used three grids of atmospheric parameters, one for each stellar component. Free parameters were [M/H], , , microturbulent velocity , and . For the surface metalicities [M/H] we used scaled solar abundances in accordance with Asplund et al. (2009). Moreover, we did a fine-adjustment for the RVs of the contributions of the three components, using the values obtained from fitting the CCF of the spectrum (Sect. 3.4) as the starting values. Normally, when we deal with an SB3 star, the flux ratios of the three stars compared to the total continuum flux are coupled via and can be determined for each grid point from least squares minimisation of two of the flux ratios. In our case, we got a much better fit to the observed composite spectrum (see Fig. 8 with the fit of the Mg I triplet) when adding an additional continuum contribution (veiling), originating from star 4. In that case, is determined from and to are no longer coupled via their flux ratios and were treated as free parameters of the fit.
|Parameter||Star 1||Star 2||Star 3||Star 5|
Notes. (a) Star 5 refers to the image of E1234. (b) Quantities in curly brackets were held fixed at plausible representative values. (c) Fraction of the light in stars 1–4.
We used all five of the above mentioned quantities as free parameters for star 3. The relatively low S/N of 32 for the single spectrum, together with the low line strengths of stars 1 and 2 compared to star 3, did not allow us to do the same for stars 1 and 2. For these, we fixed at 4.5 and 4.6, respectively, as the most probable values obtained in Sect. 5, and [M/H] at the value obtained for star 3 (the latter was done iteratively). After the final optimal solution for all free parameters was found, the parameter errors were determined from statistics using the full grid of atmospheric parameters. In this way, the errors include all interdependencies among the parameters.
A comparison of and for star 3 derived from this spectrum analysis (Table 6) with those inferred from the Monte Carlo evaluation of the system parameters (see Sect. 5) is shown in Fig. 6. The values of differ by 0.28 dex between that found from the spectrum analysis vs. that from the MC analysis. While this difference is only about 1-, it is nonetheless suggestive that star 3 may be somewhat more evolved than the MC analysis indicates. The inferred values of the flux ratio of binary ‘B’ to binary ‘A’ from the spectrum analysis is which is consistent with the value of derived from the speckle imaging (see Sect. 8.2). Overall, we take the spectrum analysis to be a satisfactory check on our more detailed parameter study in Sect. 5.
Finally, we obtained the stellar parameters for E1234 using the Spectral Parameter Classification (SPC) tool developed by Buchhave et al. (2012). SPC cross correlates an observed spectrum against a grid of synthetic spectra based on Kurucz atmospheric models (Kurucz 1992). The weighted average for the six TRES spectra of E1234 is K, , [m/H] = , km s. The values were calculated by taking an average of the stellar parameters that were calculated for each spectrum individually. The values are then weighed according to the cross correlation function peak height.
7. Tidal Synchronization Status of ‘A’ and ‘B’
The spectroscopically obtained projected rotational velocities of stars 1–3 (see Table 6) carry indirect information on the state of the tidal synchronization processes in the two eclipsing binaries. The projected synchronized (or, in the case of an eccentric system: pseudo-synchronized) rotational velocity of the binary members can be calculated simply as
where the orbital eccentricity-dependent ratio of the orbital to the (pseudo-)synchronized rotational period which, in the equilibrium tide model of Hut (1981) can be calculated as
Furthermore, we also assumed that for (pseudo-)synchronized rotation the orbital angular momentum vectors of the star and orbit are aligned. Substituting stellar radii, orbital inclinations, periods, and eccentricities from Tables 4 and 5 we arrive at , , and for stars 1, 2 and 3, respectively. Therefore, we can conclude that the stars in the 5-day circularized binary ‘A’, are rotationally synchronized with a significant likelihood. (Although, theoretically, we cannot exclude the possibility of faster stellar rotations around misaligned rotational axes.) For the primary star of the eccentric binary ‘B’, our result strongly supports a modestly super-synchronous rotation.
According to the formulae of Moreno et al. (2011) and Zahn (2008), for a convective star in an eccentric orbit with a fractional radius of , the expected synchronization time-scale even for an initially ten-times faster orbital rotation should be 1 Gyr for physically realistic values of the viscosity parameter. This timescale is shorter than the likely age we infer for the system. Therefore, we might expect some secondary source working against synchronized rotation, which might be indirect evidence for dynamical perturbations of binary ‘A’ on binary ‘B’.
8. High Resolution Imaging
8.1. Keck Adaptive Optics
To obtain better constraints on this quintuple system, we imaged both E1213 and E1234 on 2016 April 13 UT with the NIRC2 instrument (PI: Keith Matthews) on Keck II using natural guide star imaging with the narrow camera setting (10 milliarcseconds, ‘mas’, pixel). We used a three-point dither pattern to avoid NIRC2’s noisier lower left quadrant and we calibrated the images and removed artifacts using dome flat fields and dark images.
We obtained 12 images of E1213 in the band (central wavelength 2.145 m) for a total on-sky integration time of 300 seconds. Figure 9 (top panel) shows a stacked image of this target, in which we see two peaks. For each calibrated frame, we fit a two-peak PSF to measure the flux ratio and on-sky separation. We model the PSF as a combined Moffat and Gaussian shape. The best-fit PSF was found over a circular area with a radius of 10 pixels around each star (the full width at half-maximum of the PSF was about 5 pixels). More details of the method can be found in Ngo et al. (2015). When computing the separation and position angle, we applied the new astrometric corrections from Service et al. (2016) to account for the NIRC2 array’s distortion and rotation. We find the for binary ‘A’ binary ‘B’ to be , the separation to be mas, and the position angle of ‘A’ relative to ‘B’ to be deg E of N. We also took 6 images in J band (central wavelength 1.248 m). We find the for binary ‘A’ binary ‘B’ to be . The magnitude differences between the ‘A’ and ‘B’ binaries are summarized in Table 7.
We also imaged E1234 with the same setup parameters and show the stacked image in the bottom panel of Figure 9. We do not see more than one light source, so we compute a 5- contrast curve to put an upper limit on an unseen companion. To determine our limiting contrast over a range of distances, we divide the stacked image into a series of annuli centered on the star with a width of 5 pixels. Then, for every pixel in each annuli, we compute the sum of all neighboring pixels in a box. The standard deviation of these values determines the contrast for that separation. The limiting magnitude is determined by dividing each limiting contrast by the flux in a box centered on the primary star. At separations of 0.05, 0.10, 0.15, 0.20, and 0.5 from E1234, we rule out additional companions brighter than , 1.4, 2.8, 3.9, and 6, magnitudes, respectively.
8.2. Speckle Interferometry
We performed speckle interferometric observations at the 3.5-m WIYN telescope located on Kitt Peak using the Differential Speckle Survey Instrument (see Horch et al. 2009; 2011). This speckle camera provides simultaneous observations in two filters by employing a dichroic beam splitter and two identical electron-multiplying CCDs as the imagers. We observed E1213 and E1234 on 17 April 2016 UT obtaining 5 sets of 1000, 40-msec images for each star. The observations were performed simultaneously in narrow “R” and “I” bandpasses, where “R” and “I” have central wavelengths of 0.692 m and 0.880 m, respectively, and the corresponding FWHMs are 0.04 m and 0.05 m. Our final reconstructed images and results of the observations were obtained using the full combined data set. The details of how we obtain, reduce, and analyze the speckle interferometric results are described in Howell et al. (2011).
The speckle imaging also resolves the image of E1213 into two images that we associate with binaries ‘A’ and ‘B’ just as we found from the AO imaging (see Sect. 8.2). The separation on the sky in the speckle image confirms the angular separation at found in the Keck AO image. The position angle of agrees with that determined with the Keck AO to within the speckle uncertainty of 1.3. We use the speckle images to determine the ratio of the fluxes from binary ‘B’ to binary ‘A’, and find at 0.69 m, and at 0.88 m. These ratios are summarized in Table 7. When combined with the three mass functions that we measure for binaries ‘A’ and ‘B’, these flux ratios play an important role in determining the properties of the four stars (see Sect. 5).
9. Photometric Distance Estimate
We can estimate the distance to E1213 and E1234 from photometric parallax. For E1213, comprised of binaries ‘A’ and ‘B’, we have good estimates of the masses, and evolutionary states of the four constituent stars, and hence their combined luminosities. We estimate that stars 1 through 4 have a combined bolometric luminosity of . Since the three stars that dominate the luminosity budget are not far from solar-type stars, we take this ensemble of stars to be visual magnitudes brighter than the Sun (i.e., no significant bolometric correction), and thus it has . The visual magnitude of E1213 is , and therefore the distance modulus to the source is 7.1 magnitudes. Finally, this leads to a distance estimate of 260 pc. When done within our Monte Carlo parameter estimation code (see Sect. 5), we find pc.
For the stellar image E1234, we estimate that this apparently single star has K and (see Table 6). Its location is superposed on evolution tracks in the plane in Fig. 6, and one can readily see that it is substantially evolved. From the location of E1234 in the HR diagram, and the use of the Yonsei-Yale evolution tracks (Yi et al. 2001), we estimate that its bolometric luminosity is , for stars between 1.0 and 1.3 . The apparent magnitude for this star is , which roughly translates to . The distance modulus then comes out to be , corresponding to a distance of pc. This distance estimate is therefore quite consistent with that found above from the brightness of binaries ‘A’ and ‘B’.
|Parameter||Binary A||Binary B|
Notes. (a) Taken from Tables 4 and 5. (b) We apportioned the difference in velocities inversely proportional to the masses of binary A (1.82 ) and binary B (1.73 ). (c) Based on the AO imaging. (d) Based on a Monte Carlo selection of binary inclination angles and phases (see Sect. 10).
10. Binary ‘C’: Comprised of Binaries A & B
Having established the basic properties of binaries ‘A’ and ‘B’, we can now consider them as likely bound members of a higher-order binary (i.e., a quadruple star system) which we call ‘C’. We know three quantitative things about the ‘C’ binary: (1) The relative radial velocity of ‘A’ relative to ‘B’ is 10.6 km s as inferred from the difference in and , hereafter ‘’. (2) The projected angular separation between these two components is , based on the imaging study (see Sect. 8), which corresponds to a physical projected separation of AU at the estimated distance of pc. (3) The relative radial acceleration of the two binaries with respect to each other is constrained to be cm s based on limits to for the ‘A’ and ‘B’ binaries (see Table 8). We can use these three facts to set some interesting constraints on the orbital period and inclination of the ‘C’ binary following the approach used in Lehmann et al. (2016) for the case of the quadruple system KIC 7177553.
Because we have only (i) the relative radial velocity between the two binaries, , (ii) a value for the projected physical separation, , and (iii) an upper limit on the relative radial acceleration of the two binaries, we consider only the simple case of a circular orbit for the ‘C’ binary as being illustrative. These three quantities are related to the unknown orbital radius, , inclination, , and phase, by:
where is the total system mass (i.e., the combined masses of binaries ‘A’ and ‘B’). Note that is taken here to be zero at the time of superior conjunction, which is 90 from the definition used by Lehmann et al. (2015). Following Lehmann et al. (2016), this unknown orbital phase can be eliminated from Eqns. (8) and (9) to find a cubic expression for the orbital separation:
In spite of the fact that we do not know the orbital inclination angle, , we can still produce a probability distribution for (and hence ) via a Monte Carlo approach as follows. For each realization of the system we choose a random inclination angle with respect to an isotropic set of orientations of the orbital angular momentum vector of the ‘C’ binary. In addition, because there is an uncertainty in the determination of (see Tables 4 and 5), and a significant uncertainty in the projected orbital separation, (due largely to the uncertainty in distance to the target), we also choose specific realizations for those two quantities using a Gaussian random draw for both and . We then solve Eqn. (11) for , and we also store the corresponding value of . Finally, we restrict the inferred value of so as to yield a value for the radial acceleration (see Eqn. 10) that is lower than the observed upper limit cited in Table 8. If any inclination angle leads to a non-physical solution of Eqn. (11), or a radial acceleration that is too large, we discard that system.
We repeat this process times to produce our distributions of , and . Probability distributions for , , and are shown in Fig. 10. The orbital period of binary ‘C’ is years (1- uncertainties). The corresponding orbital separation is AU. According to the orbital inclination angle distribution, values of are strongly favored (see bottom panel of Fig. 10). Such orbital periods are sufficiently short that the prospects for near-term detections of the accelerated motions of binaries ‘A’ and ‘B’ are highly promising.
Given that the projected separation of binaries ‘A’ and ‘B’ on the sky is very similar (within the uncertainties) to the peak value of in the distribution of orbital separation (see Fig. 10), i.e., in Eqn. (8), this implies that or 270.
Repeating this exercise while allowing for eccentric orbits tends to broaden the distributions by a factor of about 2. We examined this question extensively in Lehmann et al. (2016). This would not change the conclusions reached here in a material way except to increase the uncertainties.
A sketch of the ‘C’ binary is shown in Fig. 11. The outer orbit of binary ‘A’ around binary ‘B’ is drawn to scale, and the orbit is assumed circular. The orbits of stars 1 and 2 in the circular binary ‘A’ and stars 3 and 4 in eccentric binary ‘B’ are drawn to the correct shape and size relative to each other; however, their absolute size has been artificially scaled up by a factor of 20 so that the orbits can be resolved.
Since we are likely viewing binary ‘C’ with large inclination angles and , we can reduce Eqns. (8) and (9) to:
where is the distance to the target, is the change in projected physical separation (in AU), and is the advance in the orbital phase (in radians) over the next few years. Future RV and high-resolution imaging observations of this object should readily achieve accuracies of 1.5 mas in terms of the angular separation of binaries ‘A’ and ‘B’, and 1 km s for the RV measurements. Solving for the orbital phase changes required to make a detection of orbital motion in binary ‘C’ we find and 0.07 orbital cycles for the high-resolution imaging and the RV observations, respectively. These orbital phase advances correspond to 2 and 4.5 years, respectively. Thus, we advocate further RV and high-resolution imaging of this object starting a couple of years from now.
11. E1213 and E1234 as a Binary Pair
The three facts we have to work with concerning the relationship between objects E1213 (comprised of a quadruple star system) and E1234 (which appears single both in imaging and spectroscopy), are the following. First, the projected separation of the two stars is on the sky. For a distance of about 260 pc, this corresponds to a projected physical separation of about 2800 AU or 0.013 pc. This is well within the allowed separation for stars to remain bound pairs without being tidally disrupted in their passage around the Galaxy (Jiang & Tremaine 2010). Second, the radial velocity of the center of mass of the ‘C’ binary (actually a quadruple system) is km s while the radial velocity of E1234 is -15.0 km s, both with an uncertainty of 1 km s. These are suggestively close to having the same value. Third, the magnitude of the difference in the two proper motion vectors is only mas yr.
Combining these three facts leads to the quite plausible, even highly likely, notion that E1213 and E1234 are themselves a gravitationally bound binary pair. We refer to this as the ‘D’ binary. The angular, velocity, and time scales are too long to hope to see any changes in this binary over the next decades.
12. E1234 as a Single Star
The rms scatter of the RVs for E1234 within the combined TRES and KO sets is 100 m s, while for the TLS observations by themselves, the rms scatter is similar (see Table 3). However, there is a systematic offset between these two sets of 900 m s. If we subtract off that constant offset from the three TLS RV measurements of E1234, the overall rms spread in the RV measurements remains at 100 m s. We are unable to account for the systematic shift between observatories, which has minimal effect on the evaluation of the orbital parameters in binaries ‘A’ and ‘B’. However, for estimating an upper limit to the projected radial acceleration of a possible binary in E1234 this discrepancy matters, and we adopt an rms radial velocity of 100 m s for this purpose. In turn, this yields an upper limit to the projected radial acceleration of 0.003 cm sec.
In addition, the angular separation of any internal stellar pairs is based on the AO and speckle imaging (see Sect. 8).
The limits on angular separation and projected radial accelerations suggest that any binary stellar pair within E1234 has a physical separation, constrained by:
where and are the orbital inclination angle and phase, respectively, of any hypothetical binary in E1234 (as in Eqns. 8 and 10). Here we have assumed a total binary mass of .
Interestingly, unless the orbital phase or inclination angle are somewhat finely tuned, the above constraints already suggest that any putative binary in E1234 must have an orbital separation not too far from 16 AU. If it is much wider it would have been resolved with AO or speckle imaging, and if much closer the acceleration (i.e., change in velocity) would likely have been detected over the 72-day interval of the RV measurements.
To push this argument a bit further, we took advantage of TRES’s 15 m s instrumental stability. Relative RVs were derived by cross-correlating each observed spectrum, order-by-order, against the strongest observed spectrum over the wavelength range 0.435-0.628 m. These yield significantly tighter constraints on changes in velocity over the course of the observations than the absolute velocities given in Table 3. With these multi-order velocities, we can set a limit of 30 m s on the constancy of the radial velocity for E1234. In turn, this lower limit increases the coefficient on the left side of Eqn. (14) to a value of 35. This constraint effectively eliminates nearly all viable solutions for the semi-major axis of any putative circular binary. All this is in keeping with our assertion that image E1234 is that of a single star.
Finally, we note that if the location of star E1234 in the is as shown in Fig. 6, it is significantly evolved, and would have a mass of and an age Gyr. If we take E1213 and E1234 to be coeval, this sets the same lower limit to the age of the two binaries in E1213.
13. Summary and Conclusions
We have presented an unusual quintuple star system consisting of 5 stars arranged as a hierarchical quadruple orbited by a single star. The K2 data show that both binaries comprising the quadruple are eclipsing (Fig. 3). We have measured the spectroscopic orbits of 3 of the 4 stars in the two binaries (Fig. 5). The 5-day binary is highly circular, while the 13-day binary is eccentric with (Fig. 5 and Tables 4 and 5). The separation of the two binaries is resolved by both Keck AO imaging and WIYN speckle interferometry at 90 mas (Fig. 9). The flux ratio of the two binaries was thereby measured quite accurately at four different wavelengths with the AO and speckle imaging (Table 7). By combining the three measured mass functions with the flux ratio of binary ‘B’ to ‘A’, and utilizing stellar evolution models, we have managed to extract reasonably accurate stellar parameters for all four stars in the two binaries (Fig. 7).
We have shown that if the RV and/or high-resolution imaging of this object are repeated 2-3 years from now, there is a good chance that changes in the motion of binary ‘A’ with respect to binary ‘B’ (i.e., within the binary ‘C’ system) will be detected. Additionally, in regard to detecting motion within the ‘C’ binary, we note that the eclipses of binary ‘A’ are deep enough to easily follow from ground-based observations. If the eclipses can be tracked with an accuracy of minutes, then we should be able to detect significant eclipse timing variations within a few years. This is physically equivalent to measuring via RV observations.
After analyzing and evaluating the quadruple system E1213, we argued that the stellar image E1234 is itself a single star (see Sect. 12), and is very likely bound to the quadruple system (see Sect. 11).
This system is one of only a relative handful of known quadruple or higher multiplicity stellar systems where the radial velocities have been measured for both binaries, and where, additionally, both binaries are eclipsing. It has a short enough outer orbital period so that motion of that longer period system can be detected within a few years.
- Ahn et al. (2012) Ahn, C.P., Alexandroff, R., Prieto, C.A., et al. 2012, ApJS, 203, 21
- Asplund (2009) Asplund, M., Grevesse, N., Sauval, A.J., & Scott, P. 2009, ARA&A, 47, 481
- Barstow et al. (2001) Barstow, M.A., Bond, H.E., Burleigh, M.R., & Holberg, J.B. 2001, MNRAS, 322, 891
- Batalha et al. (2011) Batalha, N. M., Borucki, W. J., Bryson, S. T., et al. 2011, ApJ, 729, 27
- Borkovits et al. (2015) Borkovits, T., Rappaport, S., Hajdu, T., & Sztakovics, J. 2015, MNRAS, 448, 946
- Borkovits et al. (2016) Borkovits, T., Hajdu, T., Sztakovics, J., Rappaport, S., Levine, A., Bíró, I.B., & Klagyivik, P. 2016, MNRAS, 455, 4136
- Borucki et al. (2010) Borucki, W.J., Koch, D., Basri, G., et al. 2010, Sci, 327, 977
- Buchhave et al. (2010) Buchhave, L.A., Bakos, G.Á., Hartman, J.D., et al. 2010, ApJ, 720, 1118
- Buchhave et al. (2012) Buchhave, L. A., Latham, D. W., Johansen, A., et al. 2012, Nature, 486, 375
- Conroy et al. (2014) Conroy, K.E., Prša, A., Stassun, K.G., Orosz, J.A., Fabrycky, D.C., & Welsh, W.F. 2014, AJ, 147, 45
- Cutri et al. (2013) Cutri, R.M., Wright, E.L., Conrow, T., et al. 2013, wise.rept, 1C.
- Derekas et al. (2016) Derekas, A., Plachy, E., Molnár, L., Sódor, Á., Benkő, J. M., Szabados, L., Bognár, Zs., Csák, B., Szabó, Gy. M., Szabó, R., Pál, A. 2016, submitted to MNRAS
- Di Folco et al. (2014) Di Folco, E., Dutrey, A., Le Bouquin, J.-B., et al. 2014, A&A, 565, 2
- Hadrava (1995) Hadrava, P. 1995, A&AS, 114, 393
- Hadrava (2004) Hadrava, P. 2004, in “Spectroscopically and Spatially Resolving the Components of the Close Binary Stars”, ASPC, eds. R.W. Hilditch, H. Hensberge, & K. Pavlovski, 318, 86
- Horch et al. (2009) Horch, E.P., Veillette, D.R., Baena Gallé, R., Shah, S.C., O’Rielly, G.V., & van Altena, W. 2009, AJ, 137, 5057
- Horch et al. (2011) Horch, E. P., Gomez, S. C., & Sherry, W. H. et al. 2011, AJ, 141, 45
- Howell et al. (2011) Howell, S. B., Everett, M. E., Sherry, W., Horch, E., & Ciardi, D. R. 2011, AJ, 142, 19
- Howell et al. (2014) Howell, S.B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398
- Huber et al. (2015) Huber, D., Bryson, S.T., Haas, M.R., et al. 2015, arXiv:1512.02643.
- Hut (1981) Hut, P., 1981, A&A, 99, 126
- Ilijic (2004) Ilijic, S. 2004, “Spectroscopically and Spatially Resolving the Components of the Close Binary Stars”, ASPC, eds. R.W. Hilditch, H. Hensberge, & K. Pavlovski, 318, 107
- Jiang & Tremaine (2010) Jiang, Y.-F. & Tremaine, S. 2010, MNRAS, 401, 977
- Kirk et al. (2016) Kirk, B., Conroy, K., Prsǎ, A., et al. 2016, AJ, 151, 68.
- Kupka et al. (1999) Kupka, F, Piskunov, N., Ryabchikova, T.A., Stemples, H.C., & Weiss, W.W. 1999, in “VALD-2: Progress of the Vienna Atomic Line Data Base”, 138, 119
- Kurucz (1992) Kurucz, R. L. 1992, The Stellar Populations of Galaxies: Proceedings of the 149th Symposium of the International Astronomical Union, held in Angra dos Reis, Brazil, August 5-9, 1991. Eds. B. Barbuy & A. Renzini, International Astronomical Union. Symposium no. 149, (Kluwer Academic Publishers; Dordrecht), 1992., p. 225
- LaCourse et al. (2015) LaCourse, D.M., Jek, K.J., Jacobs, T., et al. 2015, MNRAS, 452, 3561
- Lehmann et al. (2011) Lehmann, H., Tkachenko, A., Semaan, T., Gutiérrez-Soto, J., Smalley, B., Briquet, M., Shulyak, D., Tsymbal, V., & De Cat, P. 2011, A&A, 526, 124
- Lehmann et al. (20012) Lehmann, H., Zechmeister, M., Dreizler, S., Schuh, S., & Kanzler, R. 2012, A&A, 541, 105L
- Lehmann et al. (2016) Lehmann, H., Borkovits, T., Rappaport, S., Ngo, H, Mawet, D., Csizmadia, Sz., Forgács-Dajka, E. 2016, ApJ, 819, 33.
- Lohr et al. (2015) Lohr, M.E., Norton, A.J., Gillen, E., Busuttil, R., Kolb, U.C., Aigrain, S., McQuillan, A., Hodgkin, S.T., & González, E. 2015, A&A, 578, 103L
- Matijevic et al. (2012) Matijevič, G., Prša, A., Orosz, J.A., Welsh, W.F., Bloemen, S., & Barclay, T. 2012, AJ, 143, 123
- Moreno, Koenigsberger, Harrington (2008) Moreno, E., Koenigsberger, G., Harrington, D. M., 2011, A&A, 528, A48
- Ngo et al. (2015) Ngo, H., Knutson, H. A., Hinkley, S., Crepp, J. R., Bechter, E. B., Batygin, K., Howard, A. W., Johnson, J. A., Morton, T. D., Muirhead, P. S. 2015, ApJ, 800, 138
- Prša & Zwitter (2005) Prša & Zwitter 2005, ApJ, 628, 426
- Prša et al. (2016) Prša, A., et al., 2016 in preparation
- Raghavan et al. (2009) Raghavan, D., McAlister, H.A., Torres, G., et al. 2009, ApJ, 690, 394
- Rappaport et al. (2012) Rappaport, S., Deck, K., Levine, A., Borkovits, T., Carter, J., El Mellah, I., Sanchis-Ojeda, R., & Kalomeni, B. 2013, ApJ, 768, 33.
- Schutz et al. (2011) Schütz, O., Meeus, G., Carmona, A., Juhász, A., & Sterik, M.F. 2011, A&A, 533, 54
- Service et al. (2016) Service, M., Lu, J.R. , Campbell, R., Sitarski, B.N., Ghez, A.M., Anderson, J. 2016, submitted to PASP
- Shibahashi & Kurtz (2012) Shibahashi, H., & Kurtz, D.W. 2012, MNRAS, 422, 738
- Simon & Sturm (1994) Simon, K.P., & Sturm, E. 1994, A&A, 281, 286
- Skrutskie et al. (2006) Skrutskie, M.F., Cutri, R.M., Stiening, R., et al. 2006, AJ, 131, 1163.
- Slawson et al. (2011) Slawson, R.W., Prša, A., Welsh, W.F. 2011, AJ, 142, 160
- Shulyak et al. (2004) Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, C. & Weiss, W. W. 2004, A&A, 428, 993
- Still & Barclay (2012) Still, M., & Barclay, T. 2012, ascl.soft.08004
- Szentgyorgyi & Furesz (2007) Szentgyorgyi, A.-H., & Furész, G. 2007, RMxAC, 28, 1298
- Tkachenko et al. (2012) Tkachenko, A., Lehmann, H., Smalley, B., Debosscher, J., & Aerts, C. 2012, MNRAS, 422, 2960
- Tokovinin (2016) Tokovinin, A. 2016, arXiv, 160406399
- Torres (2006) Torres, G. 2006, AJ, 131,1702
- Tsymbal (1996) Tsymbal, V. 1996, in “M.A.S.S., Model Atmospheres and Spectrum Synthesis”, ASPC, 108, 198
- Yi et al. (2001) Yi, S., Demarque, P., Kim, Y.-C., et al. 2001, ApJS, 136, 417
- Zacharias et al. (2013) Zacharias, N., Finch, C.T., Girard, T.M., Henden, A., Bartlett, J.L., Monet, D.G., & Zacharias, M.I. 2013, ApJS, 145, 44
- Zahn (2008) Zahn, J.-P., 2008, in Tidal Effects in Stars, Planets and Disks, ed. M.-J. Goupil, & J.-P. Zahn, EAS Pub. Ser., 29, 67
- Zasche & Uhlar (2013) Zasche, P., & Uhlař, R. 2016, A&A, 588, 121