94 Ceti: a triple star with a planet and dust disc
94 Ceti is a triple star system with a circumprimary gas giant planet and far-infrared excess. Such excesses around main sequence stars are likely due to debris discs, and are considered as signposts of planetary systems and, therefore, provide important insights into the configuration and evolution of the planetary system. Consequently, in order to learn more about the 94 Ceti system, we aim to precisely model the dust emission to fit its observed SED and to simulate its orbital dynamics.
We interpret our APEX bolometric observations and complement them with archived Spitzer and Herschel bolometric data to explore the stellar excess and to map out background sources in the fields. Dynamical simulations and 3D radiative transfer calculations were used to constrain the debris disc configurations and model the dust emission.
The best fit dust disc model for 94 Ceti implies a circumbinary disc around the secondary pair, limited by dynamics to radii smaller than 40 AU and with a grain size power-law distribution of . This model exhibits a dust-to-star luminosity ratio of . The system is dynamically stable and N-body symplectic simulations results are consistent with semi-analytical equations that describe orbits in binary systems. In the observations we also find tentative evidence of a circumtertiary ring that could be edge-on.
keywords:stars: individual: 94 Ceti (HD 19994, HIP 14954) – stars: planetary systems – stars: binaries: general – submillimetre: stars – infrared: stars
Approximately half of the stars belong to binary or multiple systems (Duquennoy & Mayor, 1991; Lada, 2006; Eggleton & Tokovinin, 2008; Raghavan et al., 2010; Duchêne & Kraus, 2013). Stellar companions may have a large impact on planetary formation processes. In particular, they are found to truncate circumstellar discs and put a limit on the extent of material available for planetary formation (Artymowicz & Lubow, 1994; Jang-Condell et al., 2008; Andrews et al., 2010; Jang-Condell, 2015). They are also expected to stir and increase the eccentricities and relative velocities of planetesimals and planetary embryos that might have formed around the primary (Quintana et al., 2007; Thébault et al., 2009; Rafikov & Silsbee, 2015). Consequently, multiple star systems generate significant collisional activity which may shorten the lifetime of a disc and make it difficult to form planets. Observations tend to support this, since it was found that protoplanetary discs are half as common in young (a few Myr) binary stars as compared to single stars (Cieza et al., 2009; Kraus et al., 2012). In addition, these discs are much less massive compared to those found around single stars (by a factor , Harris et al., 2012) and were found to be more short-lived as well (Duchêne, 2010).
However, debris discs have been found around old binary stars, first with Spitzer (Trilling et al., 2007) and then with Herschel (Rodriguez et al., 2015). They were detected through their host star infrared excess emission, signature of the presence of micron-sized dust grains. These grains are expected to acquire unstable orbits due to both gravitational and non-gravitational perturbations. The latter scenario includes radiation pressure and Poynting-Robertson drag that effectively clear dust clouds of small and large grains (m and between 1 m to 1 mm, respectively). As such, dust in circumstellar discs must be continuously replenished from collisional processes in rings of parent bodies (planetesimals) akin to that of the asteroid and Edgeworth-Kuiper belts in our Solar System (Artymowicz, 1997; Krivov et al., 2008; Wyatt, 2008; Moro-Martin, 2013). These underlying reservoirs of large bodies show that the presence of a stellar companion does not necessarily hinder the formation of the building blocks necessary to form planets.
The presence of at least one stellar companion does not necessarily inhibit planetary formation either, since numerous planets were found in binary or multiple systems. However, the majority of these planets were found in wide binary systems (Mugrauer et al., 2005; Daemgen et al., 2009; Muterspaugh et al., 2010; Bergfors et al., 2013). In the cases where the companion was not expected to significantly affect planetary formation processes, it is not surprising that these systems were found to be hosts to planets (Eggenberger et al., 2007; Desidera & Barbieri, 2007; Duchêne, 2010) and that the properties of circumstellar protoplanetary and debris discs in binaries are found to be similar to those hosted by single stars for separations larger than several tens of AU (Kraus et al., 2012; Harris et al., 2012; Rodriguez et al., 2015).
According to the Open Exoplanet Catalogue
In addition, investigating associated debris discs and tracing the interactions between discs, planets and/or stellar companions can considerably aid our understanding of the dynamical history of a planetary system. This kind of study was carried out for the Solar System with the Nice model of Gomes et al. (2005), but also for the systems of e.g. Pictoris (see e.g. Dent et al., 2014; Nesvold & Kuchner, 2015; Millar-Blanchaer et al., 2015; Apai et al., 2015), Fomalhaut (see e.g. Beust et al., 2014; Faramaz et al., 2015; Cataldi et al., 2015; Lawler et al., 2015), HR 8799 (see e.g. Moore & Quillen, 2013; Contro et al., 2015; Booth et al., 2016), Ceti (Lawler et al., 2014), and HD 69830 (Payne et al., 2009). Therefore, studying in detail close binary or hierarchical systems, where both planets and circumstellar emission have been found, is crucial.
Five systems, out of the 23 found in triple star systems mentioned above, are associated with circumstellar dust emission. These are HD 178911 (Saffe & Gómez, 2004; Kóspál et al., 2009), Fomalhaut (Aumann, 1985; Kalas et al., 2008; Mamajek et al., 2013; Su et al., 2016), HD 40979 (e.g. Kóspál et al. 2009; Dodson-Robinson et al. 2011), 51 Eridani (e.g. Riviere-Marichalar et al. 2014), and 94 Ceti (HIP 14954, HD 19994). Another possible candidate is L1551 IRS 5, which is a young embedded binary system and potentially a triple star system (Lim & Takakuwa, 2006). Due to the low number of known triple star systems with both planets and dust emission, it is important to study each of these in great detail. We will focus here on the 94 Cet system.
The third star of 94 Cet was recently discovered by Röll et al. (2011, 2012) with astrometry, speckle interferometry, and radial velocity measurements. The system is associated with at least one circumprimary planet (Eggenberger et al., 2003), and also shows a far-infrared excess emission (Trilling et al., 2008; Eiroa et al., 2013) that possibly originates from circumstellar dust. This is a hierarchical triple star system, at a distance of pc from the Sun, where 94 Cet A is an F8 V star, and 94 Cet B and C are both M dwarfs that form a binary pair that together orbits the primary on a 2029 year long orbit.
The infrared excess of 94 Cet was detected by Trilling et al. (2008) with Spitzer/MIPS at 70 m and later confirmed, at 100 and 160 m by Eiroa et al. (2013) as part of the key project DUNES (DUst around NEarby Stars) of Herschel (Pilbratt et al., 2010). Infrared excesses at such long wavelengths originate from the thermal emission of dust particles surrounding the star. To improve the coverage of the spectral energy distribution (SED) we observed this system with APEX (the Atacama Pathfinder EXperiment) at 870 m.
Eiroa et al. (2013) found that the inferred dust temperature corresponds to a black-body radial distance from the primary star that corresponds to a dynamical unstable region due to the companion stars. Another hint that this emission may not be associated with 94 Cet A is that there is a marginally significant offset between the expected and observed position of the source, if it is co-moving with 94 Cet A. It may thus be a background source or be associated with the other stellar members of this system.
The structure of this paper is as follows; we present the stellar and system properties in Section 2, and we summarize the observations, data reduction, and the observational results in Section 3. In Sections 4 and 5, we discuss the nature of the excess, the extended emission, and background sources. Assuming that the excesses originate from disc(s), we apply both dynamical simulations and radiative transfer simulations and present these results in Section 6, and summarize our conclusions in Section 7.
2 Stellar and system properties
The binary nature of 94 Cet was first discovered by Admiral Smyth in 1836 (Raghavan et al., 2006; Smyth, 1844) who was able to resolve the stars (they share proper motion, mas yr and mas yr). The orbital parameters have been constrained and refined by Hale (1994), and more recently by Roberts et al. (2011). The companion star is a binary with two M dwarfs on a one year orbit (Röll et al., 2011, 2012) that together orbits 94 Cet A on a 2029 year long orbit. Fig. 5 shows a schematic overview of this system and orbital parameters are summarized in Table 1.
|Long. of asc.|
In the CORALIE survey, that is based on radial velocity measurements, a planetary companion around the primary star was also discovered (Queloz et al., 2004; Mayor et al., 2004), designated 94 Cet Ab. This planet has a mass of , period of 536 days, and a semi-major axis of 1.4 AU.
The properties of the two companion stars are not well-known as previous works were based on the assumption that 94 Cet is a binary. However, their masses have been estimated by Röll et al. (2011, 2012) with radial velocity measurements. These masses suggest them to be M class main sequence stars (the primary star’s age is between 1 and 5 Gyr) and so we may use other studies to infer additional properties. From Cox (2000, pp. 388-389) we find that the masses and of main sequence stars correspond to the spectral classes M0 and M3, respectively. Bessell (1991) showed that these spectral classes have effective temperatures of 3700 K and 3300 K and luminosities of and 0.01- , respectively (see also Rajpurohit et al. 2013). With these data we approximate their radii (Stefan-Boltzmann’s law) to and and estimate their surface gravity to and respectively, (see Table 2).
|Names||94 Cet A||94 Cet B||94 Cet C|
|HIP 14954 A, HD 19994|
|ICRS (J2000) R.A.||3124643|
|ICRS (J2000) Dec.||-1114596|
|Spectral class||F8 V - G0 IV||M0 V||M3 V|
|Effective temperature (K)||6187||3700||3300|
|Age (Gyr)||1.17 - 5.67||1.17 - 5.67||1.17 - 5.67|
|Proper motion (mas yr), R.A.|
|Proper motion (mas yr), Dec.|
General reference is Eiroa et al. (2013, and references therein) except when stated otherwise. Estimate from tables (Cox, 2000, pp. 388-389) using the known dynamical mass. Estimate from tables (Bessell, 1991) using the spectral class. Estimate from Stefan-Boltzmann law . Approximated from radius and mass to find a spectral model in a PHOENIX grid. Age range based on both X-ray luminosities and activity index (Eiroa et al., 2013). 94 Cet’s system proper motion is mas yr and mas yr (Eiroa et al. 2013 and SIMBAD Astronomical Database).
3 Observational data
In this work we made use of photometric data to obtain properties of the dust conforming the disc surrounding 94 Cet, and of imaging data used to make a spatial analysis of the whole system. In this section we explain how we obtained and processed the Herschel and APEX-LABOCA data plus we present the supplementary data used in our analysis. The complete observational data are presented in Table 3.
|Instrument/||Obs/Pgm ID/||Beam||Observing Date||Field centre||Offset|
|Mode||AORKEY||(m)||FWHM ()||year-mo-day||(s)||coordinates (J2000)||()|
|APEX-LABOCA||090.F-9302(A)||870||19.5||2012-08-15 –||28905||3124644, 11460|
is the on-source integration time.
Offset between measured main source position and source tabulated coordinates (J2000) of the primary star compensated for proper motion, average offset for PACS is 24 (Sánchez-Portal et al., 2014) at this observing period (OD 661).
The Herschel/PACS (Poglitsch
et al., 2010) observations were taken as part of the DUNES Open Time Key Programme and were published by Eiroa
et al. (2013)
The observations were reduced using the Herschel Interactive Processing Environment, hipe (Ott, 2010), version 13.0.0, using the PACS calibration version 69. PACS images are affected by shot noise (frequency dependent, ), consequently, we used a high-pass filter to reduce these effects. We chose radii of 20 and 25 frames for 100 and 160 m, respectively, allowing us to eliminate background structures larger than 82 and 102 arcseconds. By the nature of the procedure, high-pass filtering results in a flux loss which we rectify by masking all pixels ten times brighter than the standard deviation of non-zero flux pixels. Deglitching was done by using the second level spatial deglitching task in hipe. PACS observations were carried out at two different observational angles to decrease striping effects. In the final maps we combined individual scans for each band to increase the signal-to-noise ratio (SNR). Finally we used a drizzling method to obtain a final image scale of 1 per pixel at 100 m and 2 per pixel at 160 m compared to the native instrument pixel sizes of 32 for 100 m and 64 for 160 m.
Aperture photometry was performed using radii of 5 and 8 for 100 m and 160 m, respectively; these particular radii were used because they were found to maximise the SNR by Eiroa et al. (2013). We performed aperture and colour corrections following Balog et al. (2014). To carry out the aperture corrections we divided the integrated flux by 0.521 for 100 m and by 0.527 for 160 m. The colour corrections were done by dividing the aperture corrected flux by 1.033 for 100 m and 1.074 for 160 m, these values are correct for black body temperatures around 5000 K (Table 1 of Herschel documentation PICC-ME-TN-038, 2011).
To check whether or not our source is extended we compared the full width half maximum (FWHM) of 2D Gaussian fits with the beam FWHM. We found that only the central source was marginally extended as we can see in the radial profiles shown in Fig. 2. The profile error in this plot are estimated from the standard deviation of the radial profiles to the east, west, south, and north of the central source, and from the observed uncertainty. FWHM of Gaussians fitted on the central source was 82 and 109 at 100 and 160 m, respectively. These correspond to 185 and 246 AU at the system’s distance. The point spread function (PSF) FWHM is 68 at 100 m and 114 at 160 m with a scan speed of 20 sec. Thus we used an aperture with a radius of 6 at 100 m, and the aperture correction used on the central source was 0.595 (see also Marshall et al. 2014 for details).
The integrated flux uncertainty was estimated in a process similar to the one explained in Eiroa et al. (2013). We began by measuring the fluxes inside of forty squared apertures; each one of the same area as the aperture used on the 94 Cet. They were put at random positions, in annulii between 10 to 20″ from the source aperture, making sure to avoid all of the nearby sources. Subsequently, the standard deviation of all of these integrated fluxes was multiplied by , where is the number of background squares. Lastly, for the final error estimate, we incorporated the calibration error (5 per cent, Balog et al., 2014) with a quadratic sum.
We used the bolometer camera LABOCA (Large APEX BOlometer CAmera, Siringo et al. 2009) to observe 94 Cet, project id 090.F-9302(A). The LABOCA operating wavelength is 870 m, centred on a 150 m wide window. The 295 bolometers yield a circular field of view of 114 and the beam has a FWHM of 195. It uses different mapping modes to fill the undersampled parts of the array; we used the on-the-fly (OTF) mode, in which the maps are scanned back-and-forth linearly row-by-row or column-by-column. Skydips determined the atmospheric opacity with values between 0.19 and 0.40, with an average of and a precipitable water vapour of PWV mm.
The observations were centred on the J2000 coordinates R.A. 3124644 and Dec. 11460 and produced a 114127 sized map. Calibrations yield Jy beam V with Uranus and Neptune as calibrators
The final reduced image has a pixel resolution of 4 pixel and the flux density is given in Jy beam. Point sources’ flux densities are measured directly from their peak pixel, and the error estimate is given by the RMS of the background around each source, i.e. in a region extending between two beam radii to 60 from the source.
3.3 Ancillary Data
We took optical data from two different sources Strömgren uvby photometry from Hauck & Mermilliod (1997) and Johnson BV plus Cousins I from Perryman et al. (1997). 2MASS data was taken from Cutri et al. (2003). For the mid-infrared regime we took data from the Wide-field Infrared Survey Explorer (WISE) and the Akari satellites; the WISE was taken from the All-Sky data Release Catalogue (Wright et al., 2010) and the Akari data was taken from the Akari/IRC mid-IR all-sky Survey (Ishihara et al., 2010). We took data from the Infrared Astronomical Satellite IRAS from the IRAS Point Source Catalogue (Helou & Walker, 1988). Finally, we complemented the SED using Spitzer data. Spitzer/MIPS (Rieke et al., 2004) observations are published and summarized in more detail by Trilling et al. (2008). A Spitzer/IRS (Houck et al., 2004) spectrum (PID 102, PI: Werner, Rebull et al. 2008) extracted from the DUNES database, and not reprocessed by us, is also shown with the SED.
3.4 Observational results
Since the angular distance between 94 Cet A and BC during the observations was 24 and the 100 m beam FWHM is 68, the photometry presented here includes the fluxes from all 94 Cet components.
The photometric data is presented in Table 4 and the stellar SED is presented in Fig. 3. Excesses are clearly detected, both at 100 and 160 m (12 and 10 , respectively) and marginally at 70 m (3.4 , where denotes the error estimate at each wavelength). The photospheric spectrum was extracted from the high-resolution PHOENIX/GAIA grid (Brott & Hauschildt, 2005) by Eiroa et al. (2013) with the parameters summarized in Table 2. For the companion stars we used synthetic spectra from the PHOENIX library by Husser et al. (2013) (assuming similar metallicities as 94 Cet A, [Fe/H]0.2). We approximated the long wavelength part of the spectra (m for B and C, m for A) with Rayleigh-Jeans tails. The spectra were normalised using the V-band magnitude 11.5 for 94 Cet BC and the B and C flux density ratio of 0.29 (Röll et al., 2011, 2012) was used.
Our PACS photometry results are similar to those of Eiroa et al. (2013) (colour corrected; mJy at 100 m, and mJy at 160 m), this is caused by the differences in hipe versions, in PACS calibration trees, and in aperture corrections. Our 100 m image has a RMS between 2.0 – 2.4 mJy beam, our 160 m image has a RMS between 3.3 – 4.0 mJy beam, and they contain some background sources labeled as numbers 1-5 (1-4 are visible in Fig. 1).
|3.353||WISE (W1) (4)|
|11.561||WISE (W3) (4)|
|22.088||WISE (W4) (4)|
Total flux density of both central source and extended emission.
(1) Hauck & Mermilliod (1997). (2) Hipparcos (Perryman et al., 1997). (3) 2MASS Point Source Catalogue (Cutri et al., 2003). (4) WISE All-Sky data Release Catalogue (Wright et al., 2010). (5) Akari/IRC mid-IR all-sky Survey, II297 in Vizier, colour corrected (Ishihara et al., 2010). (6) IRAS Point Source Catalogue, II/125 in Vizier, colour corrected (Helou & Walker, 1988). (7) Spitzer/MIPS, Trilling et al. (2008). (8) Herschel/PACS, this work, colour corrected. (9) Herschel/PACS, estimates by Eiroa et al. (2013), colour corrected. (10) APEX-LABOCA, this work.
The emission that appears around the central source at 100 and 160 m may be associated with 94 Cet. We estimate the total flux density of the extended source by separating it into an eastern, southern, and western regions designated simply as E, S, and W in Fig. 1. The angular distances, from the observed stellar position, for each of these regions are 180 for E, 171 for S and 351 for W; these correspond to projected distances from the main source of 407, 387 and 794 AU, respectively. The eastern and western extensions are aligned along the projected plane of the companion star’s orbit around the primary. The southern extension is along the direction to the companion stars.
To estimate the flux densities for each region, we began by subtracting a standard PACS PSF from the main source. This PSF was normalised to the flux density of the nearby point source 94 Cet 2 to avoid underestimating the fluxes by subtracting the extended central emission. The flux densities for each region were then estimated as point sources and the results are presented in Table 5. The total flux density of the whole 94 Cet system is listed in Table 4 and plotted in Fig. 3 as red diamonds.
Among the standard PACS PSFs listed in the Herschel technical document PICC-ME-TN-033 (2015), we chose Boötis because it has the highest SNR. Additionally, Kennedy et al. (2012) found that the PSF widths varies about 2 – 4 per cent at 100 m and 1 per cent at 160 m so the effect of chosing another reference star is negligible.
4 Surrounding sources
It is imperative that we first distinguish which sources are not part of the 94 Cet system and which are likely part of it. We then focus on the latter and attempt to provide first constraints on their positions in the system, i.e., if they are due to circumstellar, circumbinary, or even circumtertiary dust emission.
4.1 Background contamination
Background confusion is an important issue at 100 and 160 m. When studying cold circumstellar dust, 100 m is the preferred wavelength because of its sensitivity to the emission of dust with temperatures of 100 K and because of the high contrast between the dust and the host star. However, dust emission from high-redshift galaxies may also have significant flux densities at this wavelength which could be misinterpreted as part of the disc emission.
Hence we can estimate the probability of contamination by background galaxies. In figure 7 of Berta et al. (2011) it is shown that the number of galaxies with a flux density greater than 2 mJy beam is about – deg, or – beam. This number decreases significantly with increased flux. For example, at 6 mJy beam (3 times the background RMS) the number of expected galaxies decreases to – deg, or – beam. In the case of the 94 Cet excess at 100 m (25 mJy) the number count has dropped off to – beam; this gives the probability of up to 2.3 per cent of coincidental alignment for a given source at these flux levels. As comparison, our estimated probability is lower than the one reported in Krivov et al. (2013), in which they calculated a 4.8 per cent probability of confusion between background galaxies and the coldest debris discs yet found.
4.2 Pointing accuracy
The average Herschel pointing offset () for this observing period (OD 661) is 236 (Sánchez-Portal et al., 2014). The last column in Table 3 refers to the observed offset, i.e. the difference between the observed position and J2000 position compensated for 94 Cet’s system proper motion.
The observed offsets of 94 Cet are and at 100 and 160 m, respectively, and while they are marginally significant, there are three scenarios which could explain them. First, the observed emission does not belong to the 94 Cet system and actually originates from a background galaxy but, as we showed earlier, the probability of this is lower than 3 per cent. Second, the Herschel pointing offset, , could have been larger for these observations, implying that our observed offsets are smaller and the observed emission actually matches the position of the primary star of the system. Third, the observed emission is actually on the expected position for the BC companion pair and originates from a circumbinary debris disc around them.
The second scenario is weakened by the analysis done in Eiroa et al. (2013), in which they found that dust emission temperature corresponds to a distance to the primary star where no stable orbits can exist due to influences from the secondary. Meanwhile, the third scenario is strengthened by the orbital measurements from Roberts et al. (2011), which locate the BC companion pair at observed offsets of and for 100 m and 160 m, respectively. As such, we continue with our analysis assuming the third scenario is correct.
5 Stable orbits and disc sizes
Multiple star systems give rise to interesting dynamics. In this system we can expect stable orbits around each of the stars, circumbinary orbits around the two secondary stars, and circumtertiary orbits around the whole system. By using the semi-analytical expression based on simulations by Holman & Wiegert (1999) we can estimate the sizes of stable regions, expressed as critical semi-major axes, around each star in a binary (see also Wiegert et al. 2014 and references therein).
We proceeded by calculating the critical semi-major axes for 94 Cet B and C independently ( and ), followed by coupling them and forming the BC binary component (). The BC component was then coupled with 94 Cet A to form the A-(BC) tertiary system ( and , see Fig. 5).
We found the following maximum radii for stable regions in the 94 Cet system, with error estimates in parentheses, and the corresponding angular sizes based on the distance to the system,
The fifth critical semi-major axis, i.e. for circumtertiary orbits, can be approximated to be at least three times the semi-major axis (Wiegert & Holman, 1997; Holman & Wiegert, 1999) which corresponds to 660 AU from the system barycentre. The western source is at an angular distance of 351 which corresponds to 794 AU at this distance. Its position angle and distance fits well with a possible circumtertiary dust ring.
Although this disc would extend on scales of several hundreds of AU, it is not excluded to have a collisionally active belt. Firstly, because significantly extended debris discs have already been observed around solar-like stars, for instance around HD 202628 where the disc extends up to AU (Krist et al., 2012), HD 207129 with a belt at AU (Löhne et al., 2012), or Ret (HIP 15371, Faramaz et al. 2014), which has a wide stellar binary companion and a resolved, asymmetric disc at AU. Secondly, the presence of a massive gravitational perturber, such as the BC pair orbiting the primary, is expected to increase the eccentricities of planetesimals inside the belt and raise the collisional activity, all the more since it is on an eccentric orbit.
Due to its large extent, this system may have suffered from perturbations from galactic tides and passing stars. Indeed, this type of perturbations is expected to induce pseudo-random losses of angular momentum on bodies orbiting their host star at large separations, whether these bodies are planetesimals (Heisler & Tremaine, 1986) or stellar binary companions (Kaib et al., 2013). Consequently, the orbits of these bodies tend to acquire small pericentres, which, in the case of planetesimals, will potentially lead to cometary activity (Kaib & Quinn, 2009). On the other hand, if it is a stellar binary companion, the survival of material around the primary may be more endangered and the architecture of any formed planetary system may be dramatically reshaped (Kaib et al., 2013). These effects are expected to be significant for old (Gyr) systems and for bodies orbiting at very large separations (between and AU). Since the separations involved in the case of the 94 Cet system are smaller than these values, this system is most probably not affected significantly by galactic tides and passing stars.
As mentioned before, Eiroa et al. (2013), using a black-body fit, attempted to retrieve the dust location. They found a disc radius of 95 AU, far inside the unstable region. Using the comparisons between black-body radii and resolved radii by Pawellek & Krivov (2015) we finally estimate a range of ‘true’ radii for this black body radius, i.e. 267 – 455 AU when taking error bars into account.
6 Numerical study
The resulting ‘true’ radii are not sufficiently accurate and therefore we resort to full radiative transfer modelling to derive more accurate values. We also use N-body symplectic simulations as a complementary approach to our analytical estimates of the locations of stable orbits in this system.
6.1 Radiative transfer modelling
We use the Monte-Carlo based program radmc-3d presented in Dullemond (2012)
|Grain size range, (m)||8 to|
|Blow out radius, (m)||1.3|
|Grain density, (g cm)||2.5|
|Size distr., (cm)||to 4|
|Surface density, (g cm)||to|
|Vertical distribution, (AU)|
|Absorption coeff., (cm g)|
|Scatter coeff., (cm g)|
|Inner disc radius, (AU)|
|Outer disc radius, (AU)|
The mass absorption coefficient, , describes how well the dust grains absorb and re-emit radiation. We obtained it from the extinction coefficients () of Miyake & Nakagawa (1993), who studied the effects of different grain sizes and size distributions of compact spherical silicate grains. The minimum grain size, that dominates the emission, is assumed to be . We will be using the 100 m data because of its better resolution and can assume minimum observed sizes between 10 and 100 m. The extinction coefficients are normalised with a gas-to-dust ratio of 100 (see Liseau et al. 2015, and references therein).
We may compare these opacities with other more recent studies: Weingartner & Draine (2001), Zubko et al. (2004) and Draine (2006) for instellar dust, and Kataoka et al. (2014) for dust aggregates in protoplanetary discs. All these works show mass absorption coefficients between 1 and 10 cm g at wavelengths around 300 m. However, they assume smaller grains than we expect in a circumstellar environment. Grains smaller than 10 m tend to exhibit up to two orders of magnitude higher absorption at shorter wavelengths ( m) and stronger silicon features than grains of, e.g., 1 mm. The final grain size range and size distribution exponent, , are also significant contributors to the inferred total dust disc mass. was first suggested by Dohnanyi (1969) but we also vary it between 2 and 4 in steps of 0.5.
The lower size limit is based on an inferred blow-out radius of the grains (), i.e., the smallest possible grains allowed due to stellar pressure forces (Plavchan et al., 2009; Wiegert et al., 2014). Simulations (Wyatt et al., 2011; Löhne et al., 2012; Thébault, 2016) show that the lower cut-off is smooth and a good approximation of the smallest allowed grain size should be around 6 times the blow-out radius. We can compute the blow-out radius by assuming a mass grain density of 2.5 g cm and obtain 1.3 m for 94 Cet A. Kirchschlager & Wolf (2013) found that only stars with effective temperatures higher than 5250 K have a determined blow-out radius and, since the companion stars are of M-type with effective temperatures under 4000 K, it is not expected for them to have a blow-out radius. Nevertheless, we assume the much brighter primary star generates a blow-out radius for the complete system. Since the larger grains do not contribute to the thermal emission of the dust, we fix the upper grain size to 1 mm; the same value as several studies which allows us to have comparable results.
The albedo () may vary widely depending on grain constituents and possible ice covering. Miyake & Nakagawa (1993) and Inoue et al. (2008) both studied the albedo of grains with different sizes and with/without ice covering. They show that the albedo for small silicate particles ( 10 m) is relatively stable between 0.5 and 0.6 and then drops down to zero at longer wavelengths between 100 to 500 m.
We estimate the optical depth at three different wavelengths: 0.5, 100 and 160 m. We use the chosen extinction coefficients and cloud densities of – g cm to obtain values of – at optical wavelengths and – at the far-infrared (FIR) wavelengths. From these values we can assume that the discs are optically thin.
We use radmc-3d, which simulates the whole wavelength range and takes scattering into account, to simulate the heating processes at optical wavelengths, where the disc is less optically thin than at FIR. With more extreme disc models (e.g. higher densities) it might even be optical thick, and we would risk overestimating the emission at FIR and underestimate the mass if we were not using radiative transfer simulations.
Similarly as in Wiegert et al. (2014) we start the simulations with a ‘standard’ disc using the grain size power-law distribution of index (homogeneously distributed in the disc) and radial disc surface density . Keep in mind that a is in reality a ring at the outer edge of the allowed radius, and a is instead a dust ring just outside the vaporisation radius. There were only small differences in the resulting SEDs with outside the range of to 2.
The accuracy of each model was inspected by eye and quantified with a reduced expressed as , where is the observed flux density, is the corresponding model flux density, and is the number of data points computed for each wavelength, that is, at 70, 100 and 160 m, as this is where the emission was found. Through an iterative process we explored the parameter space to test the degeneracy of each model.
The dust emission models are quantified by the luminosity ratio between the stellar total luminosity and the dust emission luminosity, i.e. .
Circumstellar disc SEDs
The results of radiative transfer simulations of circumstellar discs are shown in the top panel of Fig. 6. The circumprimary disc and companion star circumstellar discs were simulated separately but are shown together. The models were normalised at the 100 m flux density.
The radius of the circumprimary disc may be located at up to 54 AU. The PACS beam radius at the stellar distance corresponds to 81 and 130 AU for 100 and 160 m, respectively, which implies that the disc would not be resolved. The inner radius was set to 0.1 AU. There is a planet at 1.4 AU that could form gaps in such a disc, however the available data is insufficient for this to be visible and the models presented here do not include these effects.
The simulated circumprimary disc is in general too warm to fit the data (average temperatures are around K). Rings at the outermost edge of the stable region were too warm.
In the case of the companion pair BC, the small orbit significantly limits the radii of circumstellar discs to less than 0.2 AU. The simulated SEDs in the top panel of Fig. 6 (red curves) are based on discs with inner radius of 0.008 AU and outer radii of 0.21 and 0.16 AU for 94 Cet B and C, respectively. The circumsecondary dust is too warm to fit the observations, with average dust temperature of K.
Circumbinary disc SEDs
The circumbinary disc refers to dust orbiting both the secondary stars, 94 Cet B and C. The radial limits from the companion barycentre are 3 AU to 40 AU (01 to 18, Equation 1).
The resulting SEDs from radmc-3d are shown in the middle panel of Fig. 6. Those with were too warm. Reducing the disc into a ring in the outermost parts of the stable region did not reduce the temperature significantly as this did not increase the average distance to the primary star. However, reducing the number of small grains and increasing the number of large grains, i.e. reducing the size distribution significantly cools the disc. Thus we also show simulated SEDs with and 2.5 in the same figure.
The best fit dust model () corresponds to a disc with , , average dust temperature of K, and a luminosity ratio of . The inferred total dust mass would be . Error estimates are based on the average error of 7.9 per cent of the flux densities at 70, 100, and 160 m.
A formally better fit could probably be found with a slightly less than 3.5. However, we have too few SED measurements to constrain this better and decrease the .
Circumtertiary disc SEDs
The eastern and western extensions appear where a circumtertiary ring is expected. Fig. 7 shows a comparison between the observed emission and photometry from a simulated circumtertiary disc at 160 m based on symplectic particle simulations (see Section 6.2).
We estimate the circumtertiary flux density from the combined fluxes of the eastern and western extensions to mJy at 100 m and mJy at 160 m as plotted in the lower panel of Fig. 6.
Initially the combined east plus west flux density appear to follow the black body of dust ring with a temperature corresponding to the inner dynamical stable radius, 660 AU. The emission is however unconstrained on the Rayleigh-Jeans part of the SED and a ‘standard’ dust model with and is a sufficient input to radmc-3d.
The best fit model corresponds to a dust ring with fractional luminosity and average temperature of K. The temperature is unconstrained though and an upper limit of 30 K is a better estimate. With the total dust mass would be .
A possible upper limit of the ring radius could however be measured directly from the maximum distance between 94 Cet and the outer part of W, i.e. ″ or 860 AU.
In Fig. 7 we see that a ring would exhibit peak flux densities at the outer edges to the west and east. This is also where we see the western emission, however at the eastern edge there is no detectable emission. Furthermore, the eastern extension corresponds better with where the gap would be under the assumption that the ring is in the same orbital plane as the stars. However, an edge-on ring would emit along that region and a very clumpy ring could possibly be able to fit.
Another explanation could come from the pericentre glow effect. Pericentre glow appears when a circumstellar ring is eccentric so that the pericentre is more heated than the other parts of the ring, resulting in a horseshoe shaped source. This effect may occur due to the influence of a companion as, e.g., a planet (Wyatt et al. 1999, compare with the Fomalhaut ring, Acke et al. 2012).
6.2 N-body simulations
Symplectic integration techniques take advantage of the fact that in a planetary system, the mass of the central body is much larger than all the other ones; however, they fail if all massive bodies have comparable masses, such as in multiple stellar systems. Therefore, we used the symplectic integrator HJS
The initial conditions are summarized in Table 7. All simulations contain 10,000 test particles with semi-major axes uniformly and randomly distributed between the initial inner and outer values. Their eccentricities are randomly distributed between 0 and 0.05, in order to mimic the low eccentricity-orbits one can expect at the end of the protoplanetary phase, as well as low inclinations randomly distributed between -3 and 3 relative the outer orbital plane. The remaining initial angles, longitude of ascending node, longitude of periastron and mean anomaly, are randomly distributed between 0 and . We use values corresponding to a cold disc, which is not necessarily true given the age of the system, however, we wish here to test our analytical predictions on the location of stable orbits in the current configuration of the system, so that this simple approach should be sufficient for our purpose.
|Central star(s)||94 Cet A||94 Cet BC||94 Cet ABC|
|Initial outer semi-major axis (AU)||60||50||750|
|Final outer semi-major axis (AU)||50||40||770|
|Initial inner semi-major axis (AU)||10||3||250|
|Final inner semi-major axis (AU)||10||2.5||590|
|Number of particles|
|Length of simulation (Myr)||20||20||20|
We ran three separate simulations, circumprimary (around 94 Cet A), circumbinary (around 94 Cet B and C), and circumtertiary. In each case, a timestep of 1/20 of the smallest orbital period was used, which ensures a conservation of energy with a typical error of (Beust, 2003). The results are shown in Fig. 8. With a 20 Myr simulation length ( of the age of the system), that corresponds to 10,000 orbits of the secondary around the primary, we can assume that the system is stable.
The results are consistent with the estimates in Equation 1 (and the discs used for radmc-3d). The low density component between 600 and 670 AU probably consists of particles in the process of being ejected from the system.
In conclusion, the semi-analytical expression from Holman & Wiegert (1999) is useful also for triple systems where the component’s semi-major axes differs enough so that one can approximate the system as two separate binaries.
Multi-component systems pose many difficult questions concerning planet formation processes. 94 Cet provides one of few unique opportunities of a case study of a triple stellar system (Röll et al., 2011, 2012) with at least one planet (Eggenberger et al., 2003) and FIR excess (Trilling et al., 2008; Eiroa et al., 2013). For these reasons and a marginally significant offset in the PACS data we aimed to model the stellar excess and simulate the orbital dynamics of the system.
The central source fits well with that of a circumbinary disc around the companion pair (inside the dynamically stable radii) at K, with a dust grain size distribution between 3 and 3.5, and fractional luminosity . The disc extends from 3 AU to 40 AU from the companion stars’ barycentre with a surface density distribution that was . Assuming and the grain size range of 8 m – 1 mm this corresponds to .
The eastern and western extensions have the flux densities mJy at 100 m and mJy at 160 m, which corresponds well with the possibility of a circumtertiary ring of dust () with the fractional luminosity and temperature of K. However the uncertainty is significant and the temperature is unconstrained. The ring was assumed to extend from an inner edge at 600 – 650 AU to AU, and the initial surface density distribution of .
The system and disc configurations are stable after 20 Myr. The particle discs were simulated with the symplectic integrator HJS (Beust, 2003) using particles for each disc.
Our models provide evidence for the possibility of both circumsecondary and circumteriary dust. It is possible that the lack of circumprimary (hotter) dust emission is due to additonal planets emptying the circumprimary neighbourhood and our evidence for the circumtertiary ring is only tentative. As such it would be useful to further constrain the nature of this system through additional observations.
Even without Herschel there are a few possibilities available. ALMA, for example, could probably not observe any circumtertiary dust, but may be able to confirm our findings on the circumsecondary dust, and the LMT (Large Millimeter Telescope) in Mexico can observe 94 Cet during autumn and could reach a of 0.3 mJy in 36 hours, possibly enough to reach the dust emission.
We would like to thank R. Liseau, and also C. Eiroa, G. Kennedy, A. Krivov, J.P. Marshall, and K. Torstensson for their many comments, insights, and support for this project. We appreciate the continued support of the Swedish National Space Board (SNSB) for our Herschel projects. V. Faramaz acknowledges the support from FONDECYT Postdoctoral Fellowships, project no 3150106, and the support from the Millenium Nucleus RC130007 (Chilean Ministry of Economy). F. Cruz-Saenz de Miera is supported by CONACyT research grant SEP-2011-169554. Many thanks to Hervé Beust for his guidance using the HJS code.
This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org, the SIMBAD database, operated at CDS, Strasbourg, France, the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation, and the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
Appendix A Background galaxies
Five background sources with S/N were found in the PACS field. For the positions of the sources 94 Cet-2 and 5, we found the catalogued sources SSTSL2 J031244.01-011111.8 and SSTSL2 J031248.04-010850.9, respectively, in the Spitzer Heritage Archive
The field appears to lack previously detected sources, there is however a great number of sources to the west of the field. The lack of catalogued sources is not surprising though as the stars in the DUNES catalogue were deliberately chosen to avoid regions of noisy background as e.g. the galactic plane. In the POSS infrared images, 94 Cet-3 possibly appears, however very weakly. The other sources are either covered by the large PSF from the star or too faint. We do see the sources denoted as 94 Cet-2, 3, and 4 in the MIPS 70 m data, however, all of these have S/N . And as mentioned, the stellar proper motion between Spitzer and Herschel observations are just slightly more than 1 arcsec.
|R.A. (J2000)||Dec. (J2000)||Spectral|
|h m s||index|
|94 Cet||3 12 46.43||-1 11 51.8|
|94 Cet-1||3 12 41.56||-1 12 12.6|
|94 Cet-2||3 12 44.02||-1 11 12.5|
|94 Cet-3||3 12 44.08||-1 10 33.5|
|94 Cet-4||3 12 47.08||-1 11 22.8|
|94 Cet-5||3 12 48.05||-1 08 51.1|
Spectral index for and 160 m, and defined as . Observed coordinates, not from references. Coordinates correspond with SSTSL2 J031244.01-011111.8. Coordinates correspond with SSTSL2 J031248.04-010850.9.
|Source||R.A. (J2000)||Dec. (J2000)||(mJy)|
|h m s|
|94 Cet-L1||3 12 39.08||-1 17 06.8|
|94 Cet-L2||3 12 42.30||-1 15 32.2|
|94 Cet-L3||3 12 46.02||-1 16 16.3|
|94 Cet-L4||3 12 46.56||-1 07 00.0|
|94 Cet-L5||3 12 47.37||-1 13 20.3|
|94 Cet-L6||3 12 51.89||-1 09 36.0|
|94 Cet-L7||3 12 58.34||-1 13 02.5|
In Table 8 we clearly see that all sources have similar upper limit spectral indices (upper limits due to non-detections in the LABOCA data) as the central source. We also compare these data with larger samples of FIR galaxies from GOODS-Herschel (Elbaz et al., 2011, 2013) and Herschel-ATLAS (Rigby et al., 2011) in Fig. 9. We define the spectral index as with wavelengths from to 350, 500, or 870 m depending on available data.
However, we also see that it is difficult to distinguish the stellar excess from the background sources and other FIR galaxies by looking at the spectral indices. This is not strange as we are comparing extragalactic dust with circumstellar dust at these wavelengths.
The similarities are clearer if we compare with model SEDs of FIR galaxies (redshifts of to 2) from Chary &
Most of the PACS background sources are likely FIR galaxies. However, the central source is extended in nature and fits the positions and position angles of the stellar components, its orbit, and the expected size of a circumtertiary ring. With a 2 per cent probability of coincidental alignment we find it very improbable that the dust emission seen at 94 Cet could be attributed to a FIR background galaxy.
- thanks: This publication is based on observations with Herschel which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.
- thanks: It is also based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.
- pubyear: 2016
- pagerange: 94 Ceti: a triple star with a planet and dust disc
thanks: This publication is based on observations with Herschel which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. thanks: It is also based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.–A
- thanks: This publication is based on observations with Herschel which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.
- thanks: It is also based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.
- The DUNES data archive can be found at http://sdc.cab.inta-csic.es/dunes/jsp/masterTableForm.jsp
- Primary calibrators are Mars, Uranus, and Neptune, and secondary calibrators are listed at http://www.apex-telescope.org/bolometer/laboca/calibration/
- Comprehensive Reduction Utility for sharc-2 http://www.submm.caltech.edu/~sharc/crush/
- Available at http://ipag.osug.fr/~beusth/hjs.html
- Searched through the NASA/IPAC Infrared Science Archive, Caltech/JPL. https://irsa.ipac.caltech.edu/
- The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. http://ned.ipac.caltech.edu/
- Templates and instructions available at http://david.elbaz3.free.fr/astro_codes/chary_elbaz.html
- Acke B., et al., 2012, A&A, 540, A125
- Andrews S. M., Czekala I., Wilner D. J., Espaillat C., Dullemond C. P., Hughes A. M., 2010, ApJ, 710, 462
- Apai D., Schneider G., Grady C. A., Wyatt M. C., Lagrange A.-M., Kuchner M. J., Stark C. J., Lubow S. H., 2015, ApJ, 800, 136
- Artymowicz P., 1997, Annual Review of Earth and Planetary Sciences, 25, 175
- Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651
- Aumann H. H., 1985, PASP, 97, 885
- Balog Z., et al., 2014, Experimental Astronomy, 37, 129
- Bergfors C., et al., 2013, MNRAS, 428, 182
- Berta S., et al., 2011, A&A, 532, A49
- Bessell M. S., 1991, AJ, 101, 662
- Beust H., 2003, A&A, 400, 1129
- Beust H., et al., 2014, A&A, 561, A43
- Booth M., et al., 2016, MNRAS,
- Brott I., Hauschildt P. H., 2005, in Turon C., O’Flaherty K. S., Perryman M. A. C., eds, ESA Special Publication Vol. 576, The Three-Dimensional Universe with Gaia. p. 565 (arXiv:astro-ph/0503395)
- Cataldi G., Brandeker A., Olofsson G., Chen C. H., Dent W. R. F., Kamp I., Roberge A., Vandenbussche B., 2015, A&A, 574, L1
- Chary R., Elbaz D., 2001, ApJ, 556, 562
- Cieza L. A., et al., 2009, ApJL, 696, L84
- Contro B., Wittenmyer R. A., Horner J., Marshall J. P., 2015, preprint, (arXiv:1505.03198)
- Cox A. N., 2000, Allen’s astrophysical quantities
- Cutri R. M., et al., 2003, VizieR Online Data Catalog, 2246
- Daemgen S., Hormuth F., Brandner W., Bergfors C., Janson M., Hippler S., Henning T., 2009, A&A, 498, 567
- Dent W. R. F., et al., 2014, Science, 343, 1490
- Desidera S., Barbieri M., 2007, A&A, 462, 345
- Dodson-Robinson S. E., Beichman C. A., Carpenter J. M., Bryden G., 2011, AJ, 141, 11
- Dohnanyi J. S., 1969, \jgr, 74, 2531
- Draine B. T., 2006, ApJ, 636, 1114
- Duchêne G., 2010, ApJL, 709, L114
- Duchêne G., Kraus A., 2013, ARAA, 51, 269
- Dullemond C. P., 2012, RADMC-3D: A multi-purpose radiative transfer tool (ascl:1202.015)
- Duquennoy A., Mayor M., 1991, A&A, 248, 485
- Eggenberger A., Udry S., Mayor M., 2003, in Deming D., Seager S., eds, Astronomical Society of the Pacific Conference Series Vol. 294, Scientific Frontiers in Research on Extrasolar Planets. pp 43–46
- Eggenberger A., Udry S., Chauvin G., Beuzit J.-L., Lagrange A.-M., Ségransan D., Mayor M., 2007, A&A, 474, 273
- Eggleton P. P., Tokovinin A. A., 2008, MNRAS, 389, 869
- Eiroa C., et al., 2013, A&A, 555, A11
- Elbaz D., et al., 2011, A&A, 533, A119
- Elbaz D., et al., 2013, VizieR Online Data Catalog, 353, 39119
- Faramaz V., et al., 2014, A&A, 563, A72
- Faramaz V., Beust H., Augereau J.-C., Kalas P., Graham J. R., 2015, A&A, 573, A87
- Fuhrmann K., 2008, MNRAS, 384, 173
- Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Nature, 435, 466
- Hale A., 1994, AJ, 107, 306
- Harris R. J., Andrews S. M., Wilner D. J., Kraus A. L., 2012, ApJ, 751, 115
- Hauck B., Mermilliod M., 1997, VizieR Online Data Catalog, 2215, 0
- Heisler J., Tremaine S., 1986, \icarus, 65, 13
- Helou G., Walker D. W., eds, 1988, Infrared astronomical satellite (IRAS) catalogs and atlases. Volume 7: The small scale structure catalog Vol. 7
- Holman M. J., Wiegert P. A., 1999, AJ, 117, 621
- Houck J. R., et al., 2004, ApJS, 154, 18
- Husser T.-O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6
- Inoue A. K., Honda M., Nakamoto T., Oka A., 2008, \pasj, 60, 557
- Ishihara D., et al., 2010, A&A, 514, A1
- Jang-Condell H., 2015, ApJ, 799, 147
- Jang-Condell H., Mugrauer M., Schmidt T., 2008, ApJL, 683, L191
- Kaib N. A., Quinn T., 2009, Science, 325, 1234
- Kaib N. A., Raymond S. N., Duncan M., 2013, Nature, 493, 381
- Kalas P., et al., 2008, Science, 322, 1345
- Kataoka A., Okuzumi S., Tanaka H., Nomura H., 2014, A&A, 568, A42
- Kennedy G. M., Wyatt M. C., Sibthorpe B., Phillips N. M., Matthews B. C., Greaves J. S., 2012, MNRAS, 426, 2115
- Kirchschlager F., Wolf S., 2013, A&A, 552, A54
- Kóspál Á., Ardila D. R., Moór A., Ábrahám P., 2009, ApJL, 700, L73
- Kraus A. L., Ireland M. J., Hillenbrand L. A., Martinache F., 2012, ApJ, 745, 19
- Krist J. E., Stapelfeldt K. R., Bryden G., Plavchan P., 2012, AJ, 144, 45
- Krivov A. V., Müller S., Löhne T., Mutschke H., 2008, ApJ, 687, 608
- Krivov A. V., et al., 2013, ApJ, 772, 32
- Lada C. J., 2006, ApJL, 640, L63
- Lawler S. M., et al., 2014, MNRAS, 444, 2665
- Lawler S. M., Greenstreet S., Gladman B., 2015, ApJL, 802, L20
- Lim J., Takakuwa S., 2006, ApJ, 653, 425
- Liseau R., et al., 2015, A&A, 578, A131
- Löhne T., et al., 2012, A&A, 537, A110
- Mamajek E. E., et al., 2013, AJ, 146, 154
- Marshall J. P., et al., 2014, A&A, 565, A15
- Mayor M., Udry S., Naef D., Pepe F., Queloz D., Santos N. C., Burnet M., 2004, A&A, 415, 391
- Millar-Blanchaer M. A., et al., 2015, ApJ, 811, 18
- Miyake K., Nakagawa Y., 1993, \icarus, 106, 20
- Montesinos B., et al., 2016, preprint, (arXiv:1605.05837)
- Moore A., Quillen A. C., 2013, MNRAS, 430, 320
- Moro-Martin A., 2013, Dusty Planetary Systems. p. 431, doi:10.1007/978-94-007-5606-9˙9
- Mugrauer M., Neuhäuser R., Seifahrt A., Mazeh T., Guenther E., 2005, A&A, 440, 1051
- Muterspaugh M. W., et al., 2010, AJ, 140, 1657
- Nesvold E. R., Kuchner M. J., 2015, ApJ, 815, 61
- Ott S., 2010, in Mizumoto Y., Morita K.-I., Ohishi M., eds, Astronomical Society of the Pacific Conference Series Vol. 434, Astronomical Data Analysis Software and Systems XIX. p. 139 (arXiv:1011.1209)
- Pawellek N., Krivov A. V., 2015, MNRAS, 454, 3207
- Payne M. J., Ford E. B., Wyatt M. C., Booth M., 2009, MNRAS, 393, 1219
- Perryman M. A. C., et al., 1997, A&A, 323
- Pilbratt G. L., et al., 2010, A&A, 518, L1
- Plavchan P., Werner M. W., Chen C. H., Stapelfeldt K. R., Su K. Y. L., Stauffer J. R., Song I., 2009, ApJ, 698, 1068
- Poglitsch A., et al., 2010, A&A, 518, L2
- Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., Roush T., Fong W., 1994, ApJ, 421, 615
- Queloz D., Mayor M., Naef D., Pepe F., Santos N. C., Udry S., Burnet M., 2004, in Penny A., ed., IAU Symposium Vol. 202, Planetary Systems in the Universe. p. 106
- Quintana E. V., Adams F. C., Lissauer J. J., Chambers J. E., 2007, ApJ, 660, 807
- Rafikov R. R., Silsbee K., 2015, ApJ, 798, 69
- Raghavan D., Henry T. J., Mason B. D., Subasavage J. P., Jao W.-C., Beaulieu T. D., Hambly N. C., 2006, ApJ, 646, 523
- Raghavan D., et al., 2010, ApJS, 190, 1
- Rajpurohit A. S., Reylé C., Allard F., Homeier D., Schultheis M., Bessell M. S., Robin A. C., 2013, A&A, 556, A15
- Rebull L. M., et al., 2008, ApJ, 681, 1484
- Rieke G. H., et al., 2004, ApJS, 154, 25
- Rigby E. E., et al., 2011, MNRAS, 415, 2336
- Riviere-Marichalar P., et al., 2014, A&A, 565, A68
- Roberts Jr. L. C., Turner N. H., ten Brummelaar T. A., Mason B. D., Hartkopf W. I., 2011, AJ, 142, 175
- Rodriguez D. R., Duchêne G., Tom H., Kennedy G. M., Matthews B., Greaves J., Butner H., 2015, MNRAS, 449, 3160
- Röll T., Seifahrt A., Neuhäuser R., Köhler R., Bean J., 2011, in EAS Publications Series. pp 429–432, doi:10.1051/eas/1045076
- Röll T., Neuhäuser R., Seifahrt A., Mugrauer M., 2012, A&A, 542, A92
- Saffe C., Gómez M., 2004, A&A, 423, 221
- Sánchez-Portal M., et al., 2014, Experimental Astronomy, 37, 453
- Siringo G., et al., 2009, A&A, 497, 945
- Skrutskie M. F., et al., 2006, AJ, 131, 1163
- Smyth W. H., 1844, A cycle of celestial objects
- Su K. Y. L., Rieke G. H., Defrére D., Wang K.-S., Lai S.-P., Wilner D. J., van Lieshout R., Lee C.-F., 2016, ApJ, 818, 45
- Thébault P., 2016, A&A, 587, A88
- Thébault P., Marzari F., Scholl H., 2009, MNRAS, 393, L21
- Trilling D. E., et al., 2007, ApJ, 658, 1264
- Trilling D. E., et al., 2008, ApJ, 674, 1086
- Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296
- Wiegert P. A., Holman M. J., 1997, AJ, 113, 1445
- Wiegert J., et al., 2014, A&A, 563, A102
- Wright E. L., et al., 2010, AJ, 140, 1868
- Wyatt M. C., 2008, ARAA, 46, 339
- Wyatt M. C., Dermott S. F., Telesco C. M., Fisher R. S., Grogan K., Holmes E. K., Piña R. K., 1999, ApJ, 527, 918
- Wyatt M. C., Clarke C. J., Booth M., 2011, Celestial Mechanics and Dynamical Astronomy, 111, 1
- Zubko V., Dwek E., Arendt R. G., 2004, ApJS, 152, 211