THE THREE-DIMENSIONAL ARCHITECTURE OF THE υ ANDROMEDAE PLANETARY SYSTEM

The υ Andromedae (υ And) planetary system was the first discovered multi-exoplanet system around a main-sequence star and is possibly still the most studied multi-planet system other than our solar system. υ And b was discovered using the radial velocity (RV) technique by Butler et al. (1997) at the Lick Observatory. Two years later, combined data from Lick and the Advanced Fiber-Optic Echelle spectrometer (AFOE) at the Whipple Observatory revealed the presence of two additional planets, υ And c and d (Butler et al. 1999). Follow-up by François et al. (1999) confirmed that the existence of planets was the best explanation for the RV variations. Even at the time of Butler et al. (1999), the semi-major axes of the planets (0.059, 0.83, and 2.5 au), the minimum masses (0.71, 2.11, and 4.61 M_Jup), and eccentricities (0.034, 0.18, and 0.41) made it clear that this system was very unlike our solar system, and it has presented a challenge to planet formation models that explain the solar system. Stepinski et al. (2000) confirmed the presence of these three RV signatures in the existing data using two different fitting algorithms, but stressed that the eccentricity of planet c was poorly constrained by the existing data.

Thus far, astrometry is one of the few techniques that can be used to break the msin i degeneracy in the RV method for nontransiting planets (the other is a relatively new technique that uses high-resolution spectra to directly observe the RV of the planet; for example, see Brogi et al. 2012; Rodler et al. 2012). Astrometry is the process of measuring a star's movement on the plane of the sky, and hence provides two-dimensional information that is orthogonal to RV. Because this measurement is made relative to other objects in the sky, it is extremely difficult to obtain the high precision necessary to detect planets. For small, close-in planets, the necessary precision is in the μas range (Quirrenbach 2010), since astrometry is more sensitive to planets with relatively large masses, low inclinations, and large semi-major axes. Nonetheless, Mazeh et al. (1999) reported a small, positive detection in the Hipparcos (HIP) data of an astrometric signal at the period of planet d and derived an inclination of 156° and a mass of 10.1 M_Jup. However, Pourbaix (2001) demonstrated that astrometric fits to the HIP data for 42 stars, including υ And, were not significantly improved by the inclusion of a planetary orbit, and that the inclinations for planets c and d could be statistically rejected. Reffert & Quirrenbach (2011) reanalyzed the HIP data and placed an upper mass limit on both planets c and d of 8.3 M_Jup and 14.2 M_Jup, but did not claim true masses.

υ And became the first multi-planet system to have a positive astrometry detection above 3σ when McArthur et al. (2010) detected the orbits of planets c and d using the Hubble Space Telescope (HST). Their orbital fit included all previously obtained RVs (including re-reduced Lick data), and added RVs from the Hobby–Eberly Telescope at the McDonald Observatory. The combined astrometry+RV fit did not converge when planet b was included, indicating that planet b presents no astrometric signal. Indeed, using their Equation (8), which relates the RV and astrometry, planet b would be expected to have a signal of α ∼ 40 μas at an inclination of ∼3°, well below HST's detection limit of 0.25 mas. This nondetection puts a weak upper mass limit on planet b of ∼78 M_Jup, as the planet would be astrometrically detectable by HST at inclinations below ∼05.

McArthur et al. (2010) found that planets c and d have inclinations of 7868 ± 1003 and 23758 ± 1316, respectively, relative to the plane of the sky, with corresponding masses of $13.98^{+2.3}_{-5.3} \,M_{{\rm Jup}}$ and $10.25^{+0.7}_{-3.3} \,M_{{\rm Jup}}$. The mutual inclination between the two planets is 29917 ± 1°. This value is quite unlike any mutual inclination found among the planets of our solar system. Subsequent dynamical studies suggest that this may be the result of a three-body planet–planet scattering scenario in which one planet is ejected from the system (Barnes et al. 2011; Libert & Tsiganis 2011).

By looking at the infrared excess in the stellar spectrum attributed to planet b, Harrington et al. (2006) attempted to chart the phase offset, i.e., the angle between the hottest point on the surface and the sub-stellar point. Assuming a radius of <1.4 R_Jup, the observed amplitude of flux variation demanded that the planet's inclination must be >30° from the plane of the sky. Unfortunately, this study was based on only five epochs of data over a single orbit, and thus does not provide tight constraints. The infrared phase curve of υ And b was revisited with seven additional short epochs and one continuous ∼28 hr observation by Crossfield et al. (2010). The picture presented in Crossfield et al. (2010) was consistent with Harrington et al. (2006), but they allowed larger radii in their model, finding that the inclination must be >28° for a 1.3 R_Jup planet and >14° for a 1.8 R_Jup.

There is marginal evidence for a fourth planet orbiting υ And. McArthur et al. (2010) found an improvement in their fit when a linear trend indicative of a longer period planet was included. Later, Curiel et al. (2011) found a signal at 3848.9 days using the Lick (Fischer et al. 2003; Wright et al. 2009) and ELODIE radial velocities (Naef et al. 2004). These authors have taken this to be a fourth planet in the system, as suggested by McArthur et al. (2010). The McArthur et al. (2010) analysis used re-reduced Lick data (D. Fischer & M. Giguere 2009, private communication) for their combined RV and astrometry orbital fit. As explained in McArthur et al. (2014), these data include updated γ values (constant velocity offsets in the RVs) that removed this signal from the Lick data. Later, Tuomi et al. (2011) analyzed the older published RV data sets (Fischer et al. 2003; Wright et al. 2009) and also found a period for this fourth planet (that was an artifact of the missing γ) of 2860 days, but noted that the data sets seemed inconsistent. Tuomi et al. (2011) performed fits and calculated the Bayesian inadequacy criterion for the individual data sets. They found that the Wright et al. (2009) Lick data produced a significantly different period for planet e (3860 days) and the Bayesian inadequacy criterion indicates that this data set has a >0.999 probability of being inconsistent with the other data sets. While a longer period planet, indicated by a small slope in the radial velocities, may exist in this system, the fourth planet signal reported by Curiel et al. (2011) was a product of the earlier reduction of the Lick data, which did not account for an instrument change that caused a shift in the γ. For this reason, we do not include this planet in our study.

The rotation and obliquity of υ And A are also of interest in this study (see Section 3.1). Measurements of vsin i = 9.6 ± 0.5 km s⁻¹ (Valenti & Fischer 2005) and stellar radius $R_{\star } = 1.64^{+0.04}_{-0.05} \,R_{\odot }$ (Takeda et al. 2007) limit the rotation period to ≲ 8 days for physical values of i; however, the only measured period consistent with these data is 7.3 days (Simpson et al. 2010). This period suggests an obliquity i ∼ 60° (measured from the sky plane), but the signal at this period is very weak and it is impossible to distinguish this obliquity from i ∼ 120° (spinning in the opposite sense).

Numerous dynamical studies of υ Andromedae have been performed, both numerical and analytical. Early studies focused on the stability of the system using N-body models and analytic theory, prior to the astrometric detection of planets c and d (McArthur et al. 2010). These studies showed that the stability of the system is highly sensitive to the eccentricities of planets (d in particular; Laughlin & Adams 1999; Barnes & Quinn 2001, 2004), the relative inclinations of the planets (Rivera & Lissauer 2000; Stepinski et al. 2000; Chiang et al. 2001; Lissauer & Rivera 2001; Ford et al. 2005; Michtchenko et al. 2006), their true masses (since only minimum masses were known prior to McArthur et al. 2010; Rivera & Lissauer 2000; Stepinski et al. 2000; Ito & Miyama 2001), the effect of general relativity on planet b's eccentricity (Nagasawa & Lin 2005; Adams & Laughlin 2006; Migaszewski & Goździewski 2009), and the accuracy of the RV data (Lissauer 1999; Stepinski et al. 2000; Goździewski et al. 2001), and also found that there were stable regions only between planets b and c and exterior to planet d (Rivera & Lissauer 2000; Lissauer & Rivera 2001; Barnes & Raymond 2004; Rivera & Haghighipour 2007). The dynamical study of McArthur et al. (2010) found stable configurations for all three planets on a timescale of 10⁵ yr and constrained the inclination of planet b to i < 60° or i > 135°. It has also been demonstrated that analytical or semi-analytical theory does not adequately describe the dynamics of the system (Veras & Armitage 2007), unless taken to very high order in the eccentricities (Libert & Henrard 2007), though these studies assumed coplanarity since the inclinations of the planets were not known at the time.

Other studies dealt with the evolution and formation of certain features of the system, in particular, the apparent alignment of the pericenters (or "apsidal alignment," Δϖ, noted by Rivera & Lissauer 2000; Chiang et al. 2001) and large eccentricities of planets c and d. Some have investigated whether the present day system could have been produced by interactions with a dissipating disk (Chiang & Murray 2002; Nagasawa et al. 2003), by planet–planet scattering (Nagasawa et al. 2003; Ford et al. 2005; Barnes & Greenberg 2007; Ford & Rasio 2008; Barnes et al. 2011), by interactions with the stellar companion, υ And B (Barnes et al. 2011), by secular or resonant orbital evolution (Jiang & Ip 2001; Malhotra 2002; Ford & Rasio 2008; Libert & Tsiganis 2011), or by accelerations acting on the host star (Namouni 2005). Michtchenko & Malhotra (2004) found that Δϖ can be in a state of circulation, libration or a "nonlinear secular resonance," all within the observational uncertainty. In short, the dynamical evolution of the system appears to be highly sensitive to the initial conditions, and planet–planet scattering appears to be the most promising explanation for its current state.

In particular, McArthur et al. (2010) noted that the true masses of planets c and d naturally resolve a difficulty in explaining the system's formation. When only the minimum masses of the two planets were known, it was generally assumed in dynamical analysis that planet d was the larger (msin i_c = 1.8898 M_Jup and msin i_d = 4.1754 M_Jup); however, mechanisms that can excite the eccentricities of the planets, such as resonance crossing (Chiang et al. 2002) or close encounters (Ford et al. 2001), tend to result in the smaller mass planet having the larger eccentricity. Some authors noted that because planet d is observed to have the larger eccentricity (e_c = 0.245 ± 0.006, e_d = 0.316 ± 0.006; McArthur et al. 2010), the formation of the system may have required the presence of a gas disk or the ejection of an additional low-mass planet (Chiang et al. 2002; Rivera & Lissauer 2000; Ford et al. 2005; Barnes & Greenberg 2007).

Finally, Burrows et al. (2008) produced theoretical pressure–temperature profiles and spectra for several close-in giant planets, including υ And b, and compared these results with the phase curve in Harrington et al. (2006). Unfortunately, the phase curve data were too sparse to make any conclusions about the size or structure of the atmosphere or the planet's inclination, only that a range of radii and inclinations are consistent with the data, and that a temperature inversion in the atmosphere is also consistent. The authors suggested that more frequent observations and multi-wavelength data may break the degeneracies in their model.

Observations account for the inclinations and true masses of planets c and d, but the inclination and true mass of planet b remain undetermined. Early dynamical studies found stable regions of parameter space for the coplanar, three-planet system; however, stable configurations have not previously been identified for three planets over long timescales since the large mutual inclination between planets c and d was discovered by astrometry. Additionally, it is unclear whether the phase curve observations (Harrington et al. 2006; Crossfield et al. 2010) are consistent with the RV+astrometry observations (McArthur et al. 2010). Here we explore all of the above issues using dynamical models.

This work is a sweep through parameter space for stable configurations for all three planets, following up on the dynamical analysis in McArthur et al. (2010). An acceptable configuration in this study is one that (1) satisfies the RV+astrometry fit (McArthur et al. 2010), (2) is dynamically stable, and (3) is consistent with the IR phase curve measurements (Crossfield et al. 2010). Our results satisfy all three criteria; however, reconciliation with the phase curve measurements requires planet b to have an inflated radius and motivates us to include tidal heating in this study.

Section 2 describes the methods that we use to explore parameter space, model the dynamics, and estimate tidal heating in planet b. Section 3 describes the results we obtain for stability and system evolution. Section 4 focuses on tidal heating and reconciliation with Crossfield et al. (2010). In Section 5, we discuss our results in the context of previous studies, and summarize our conclusions in Section 6.

As in McArthur et al. (2010), all RV data sets are re-examined using the published errors (which are large for the AFOE data, in particular). The Lick data set has been re-reduced since the original publications, resulting in new γ offsets. This updated set, published recently in Fischer et al. (2014), was used in McArthur et al. (2010; D. Fischer & M. Giguere 2009, private communication), and similarly we use it here. The RV fit is performed on each data set individually and compared with the other data sets for consistency, and large outliers are examined and removed, if necessary.

Trials are generated by drawing randomly from within the uncertainties of the χ² best fit to the RV and astrometric data. Table 1 shows the parameter space explored. The astrometric constraints on the semi-major axes of planets c and d are ignored as the RV derived periods provide much stricter constraints on these quantities. Trials are not generated taking into account interdependencies of the various parameters; however, because of the small size of the uncertainties, all trials should be faithful to the data. Nevertheless, we include χ² values for our stable cases to confirm consistency. This work is not meant to represent a complete analysis of all possible configurations. We are merely establishing the existence of stable cases within the observational constraints.

The nominal eccentricity of planet d is 0.316; however, we find that very few trials are stable above e_d ∼ 0.3 (see Figure 1), so we apply a much looser constraint to the lower bound of planet d's eccentricity, drawing from a uniform distribution across the domain 0.246 < e_d < 0.322, rather than a Gaussian distribution.

Zoom In Zoom Out Reset image size

For the stability analysis, we use HNBody (Rauch & Hamilton 2002), which contains a symplectic integrator for central-body-like systems, i.e., systems in which the total mass is dominated by a single object. This symplectic scheme alternates between Keplerian motion and Newtonian perturbations at each time step (see Wisdom & Holman 1991). During one half-step (the "kick" step), all gravitational interactions are calculated and the momenta are updated accordingly. During the other half-step (the "drift" step), the system is advanced along Keplerian (two-body) motion, using Gauss's f and g equations (see Danby 1988). The entire integration is done in Cartesian coordinates. Unlike Mercury (Chambers 1999), post-Newtonian (general relativistic) corrections are included as an optional parameter in HNBody, and we utilize them here.

While HNBody is fast and its results compare well with results from Mercury, its definitions of the osculating elements used at input and output differ slightly from Mercury's. The mass factor used in the definitions of semi-major axis and eccentricity (for astrocentric or barycentric elements) does not include the planet's mass; in other words, the planet is treated as a zero mass particle during input and output conversions between Cartesian coordinates and osculating elements. Because of this, the use of osculating elements during input can result in incorrect periods and Cartesian velocities. For most planetary systems, which have poorer constraints on the periods and semi-major axes of the planets, this aspect makes little difference, but for υ Andromedae the periods are known with the high precision that comes with >15 yr of RVs, and so all of our input and output from HNBody is done in Cartesian positions and velocities, to ensure proper orbital frequencies.

A second important conversion must be done to account for a difference in units. HNBody and Mercury enforce the relationship G M_☉ D² au⁻³ = k², where k is the Gaussian gravitational constant (defined to be 0.01720209895 au$^{3/2} \,M_{\odot }^{-1/2} \,D^{-1}$), D is the length of the day, and au is the astronomical unit, based on the IAU definitions prior to 2012. The current accepted IAU units fix the au to be exactly 149, 597, 870, 700 m and the constant k is no longer taken to be a constant value. This redefinition was done for ease of use and to allow for reconciliation with the length contraction and time dilation of Einstein's relativity. Because these N-body models were developed prior to this redefinition, k was used as a fundamental constant, as that allowed for better accuracy in solar system integrations when the au was uncertain. We believe these integrators still accurately represent the dynamics of planetary systems, since these models do not attempt to account for length contraction and time dilation and in any case these effects will be very small. Note that HNBody does account for relativistic precession, which is important for the motion of planet b, as mentioned in Section 1.2.

The units utilized by the integrator are then au, days, and solar masses. For observations of exoplanets, SI units are more sensible, but the values of G, M_☉, and au need not obey the above constraint. Hence the simulated orbital periods will not be equal to the measured orbital periods unless we first convert from SI units, which do not obey this constraint, to IAU units, which do. We accomplish this by choosing a value for the au that satisfies G M_☉ D² au⁻³ = k², given the G, M_☉, and D used in the model of the observations. We then check that this definition of the au correctly reproduces the orbital periods of the planets when calculated using k and Kepler's third law,

$$\begin{equation} T^2 = \frac{4 \pi ^2}{k^2(M_{\star }+m_p)} a^3, \end{equation} \tag{ 1 }$$

where T is the planet's period in days, M_⋆ and m_p are the masses of the star and planet in solar masses, and a is the planet's semi-major axis in au. Finally, to verify that HNBody "sees" the correct periods, we run two-body integrations for each planet at a high time resolution and performed fast Fourier transforms on the planets' velocities, confirming that this approach is the most accurate.

For the reasons described above, we take the orbital elements from the RV+astrometry, convert to Cartesian coordinates for dynamical modeling, then convert back to orbital elements for stability analysis. For analysis of the orbital evolution, we convert from Cartesian "line-of-sight" coordinates to the Cartesian invariable plane (the plane perpendicular to the total angular momentum of the system) coordinates prior to the conversion back to orbital elements. The invariable plane of the system is inclined by ∼15° from the sky, so inclinations measured from this plane are similar to those measured from the sky plane; see the Appendix for a potential pitfall in modeling astrometrically measured orbits.

The initial parameter space includes inclinations of planet b from 0° to 180°, which is to say that any inclination is consistent with the observations. Extremely low inclination, i_b < 1° or i_b > 179°, places the mass of planet b in the brown dwarf range. Inclinations ≳ 90° are retrograde orbits with respect to the orbits of the outer two planets (note that the outer system's invariable plane is inclined only ∼15° from the sky plane). The observations allow for such configurations, and we cannot rule them out based on dynamical stability.

We use a constant phase-lag (CPL) model (Darwin 1880) to estimate the amount of tidal energy that could be dissipated in the interior of υ And b. This model is described in detail in Appendix E of Barnes et al. (2013; see also Ferraz-Mello et al. 2008). In short, the model treats the tidal distortion of the planet as a superposition of spherical harmonics with different frequencies, which sum to create a tidal bulge that lags the rotation by a constant phase. The strength of tidal effects is contained in the parameter Q, called the tidal quality factor, estimated to be ≳ 10⁴ for the planet Jupiter (Goldreich & Soter 1966; Yoder & Peale 1981; Aksnes & Franklin 2001). For close in Jovians like υ And b, Ogilvie & Lin (2004) suggest Q could be as high as 5 × 10⁷, hence we model a range of Q values from 10⁴ to 10⁸.

The CPL model, commonly used in the planetary science community, is only an approximate representation of tidal evolution (for an in-depth discussion of the limitations of the model, see Efroimsky & Makarov 2013). However, at the eccentricities we explore here (e ≲ 0.15), results from this model are qualitatively similar to other tidal models. The CPL model has the advantage of fast computation and is accurate enough for our purpose here, which is merely to establish the possibility of tidal heating in the interior of planet b. Given the uncertainty in the model and in the properties of the planet, our results should not be taken to be a precise calculation of the planet's tidal conditions.

We run our initial set of 1000 trials for 1 Myr and flag as unstable those trials in which one or more planets were lost. The resulting stability maps are shown in Figure 1. Here we show the most relevant parameters, i.e., those with the largest uncertainty (Ω_b and i_b) and those that have a large effect on stability (e_c and e_d). The coordinate system used here is the RV+astrometry ("line-of-sight") coordinate system, in which i is measured from the plane of the sky and Ω is measured counterclockwise from north.

We note regions of greater stability concentrate around inclinations for planet b of ≲ 40° and ≳ 140° and ascending nodes of ∼0°. A chasm of instability lies across inclinations of ∼60° to ∼140°.

Next, we examine the eccentricity and inclination evolutions for all trials in which three planets survived 1 Myr. Many of these "stable" cases exhibit chaotic evolution, with planet b even reaching eccentricities of ∼0.9. We assume trials with chaotic evolution are in the process of destabilizing and should be discarded. This leaves us with ∼30 trials which we then ran for 100 Myr. Trials in which all planets survive 100 Myr integrations with no chaotic evolution or systematic regime changes are considered robustly stable. These 10 cases are plotted as x's in Figure 1.

υ Andromedae is estimated to be ∼3 Gyr old (Takeda et al. 2007); our ideal goal would thus be to demonstrate stability over this time span; however, because a very small time step is necessary to resolve the 4.6 day orbit of planet b, simulations of the system on spans of gigayears are computationally prohibitive. Hence, we limit ourselves to a domain of 100 Myr (1/30th of the system's lifetime), and note that this length of time corresponds to ∼8 billion orbits of planet b and that no previous study has been able to show stability for all three planets on this timescale.

The χ² results for our stable cases are shown in Table 2. Here, χ² represents the goodness of fit of each configuration to the data, and would ideally be equal to the number of degrees of freedom (DOF) in the model of the data. Configuration PRO1 is chosen as our nominal case because it has the lowest χ² value of the prograde (orbit of planet b) cases.

For our four prograde, stable trials, we generate the stability maps shown in Figure 2 (initial conditions shown in Table 3). Keeping all other parameters constant, we vary the inclination and ascending node of planet b to further explore these "islands" of stability. In order to keep all of our cases consistent with the observations of planet b, we adjust its mass with changes in inclination and subsequently must adjust its semi-major axis to maintain the observed period. Thus changes in inclination imply not only a different mass via the msin i degeneracy, but also a change in the semi-major axis.

Zoom In Zoom Out Reset image size

In all cases, we see that our solutions occupy stable regions of phase space. PRO2 and PRO3 appear to be perched on the edges of two large stable regions, while PRO1 occupies a very narrow stable "inclination stripe" at ∼5°. As in Figure 1, inclination in Figure 2 is measured from the sky plane, not the invariable plane of the system.

An additional complication to the dynamics of the system is the oblateness of the host star. Migaszewski & Goździewski (2009) showed the importance of J₂ (the leading term in the gravitational quadrupole moment) of the star in their secular analysis, though they used a stellar radius of R = 1.26 R_☉, significantly smaller than the current best measurement of R = 1.631 ± 0.014 R_☉ (Baines et al. 2008). To verify the importance of oblateness, we simulated our best prograde χ² case (PRO1), varying J₂ from 10⁻⁵ to 10⁻² and R_⋆ from 1.26 R_☉ to 1.63 R_☉. We find that values of J₂ ≳ 10⁻³ cause the system to become unstable, and lower values significantly change the eccentricity evolution of planet b. Unfortunately, the J₂ value of the star is not known and there exists some disagreement regarding its radius, and thus a detailed exploration of parameter space including these two additional parameters is beyond the scope of this work. Therefore, in our primary analysis, the quadrupole moment is ignored (J₂ = 0).

The eccentricity and inclination evolutions for our prograde trials are shown in Figures 3–6. In these figures, inclination is measured from the invariable plane, that is, a plane perpendicular to the total angular momentum vector of the system. We see that for all cases, the evolution of the eccentricities and inclinations is periodic for at least 100 Myr (∼8 billion orbits of planet b). At a glance, one might expect for PRO1 to lose planet b; however, as shown in Figure 3, the pattern seen in its eccentricity evolution is repeated reliably over many orbits, and therefore the configuration is robustly stable.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For planets c and d, Figures 7 and 8 show Δϖ = ϖ_d − ϖ_c (the difference of the longitudes of pericenter) and their mutual inclination Ψ_cd for our cases in which planet b is in prograde motion. We see that Δϖ undergoes circulation in cases PRO1 and PRO2, and librates around anti-alignment in cases PRO3 and PRO4, although it is very close to the separatrix in both. The amplitudes of libration are ∼240° and ∼210°, for PRO3 and PRO4, respectively, and rms values about the libration center (180°) are 55° and 47°, respectively. The mutual inclination between planets c and d oscillates between ∼30° and ∼40° in cases PRO2, PRO3, and PRO4. The angle explores a slightly wider range in case PRO1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size