A DETAILED ANALYSIS OF THE HD 73526 2:1 RESONANT PLANETARY SYSTEM^*

AAPS Doppler measurements are made with the UCLES echelle spectrograph (Diego et al. 1990). An iodine absorption cell provides wavelength calibration from 5000 to 6200 Å. The spectrograph point-spread function (PSF) and wavelength calibration are derived from the iodine absorption lines embedded in every pixel of the spectrum by the cell (Valenti et al. 1995; Butler et al. 1996). The result is a precision Doppler velocity estimate for each epoch, along with an internal uncertainty estimate, which includes the effects of photon-counting uncertainties, residual errors in the spectrograph PSF model, and variation in the underlying spectrum between the iodine-free template and the epoch spectra observed through the iodine cell. All velocities are measured relative to the zero-point defined by the template observation. A total of 36 Anglo-Australian Telescope (AAT) observations have been obtained since 1999 February 2 (Table 1) and were used in the following analysis, representing a data span of 4836 days. The exposure times range from 300 to 900 s, and the mean internal velocity uncertainty for these data is 4.1 m s⁻¹.

Since HD 73526 is among the faintest AAPS targets (V = 9.0), we have obtained supplemental observations with the Planet Finder Spectrograph (PFS; Crane et al. 2006, 2008, 2010) on the 6.5 m Magellan II (Clay) telescope. The PFS is a high-resolution (R ∼ 80, 000) echelle spectrograph optimized for high-precision, radial velocity measurements (e.g., Albrecht et al. 2011, 2012; Anglada-Escudé et al. 2012; Arriagada et al. 2013). The PFS also uses the iodine cell method, as described above, to obtain precise radial velocities. The 20 measurements of HD 73526 are given in Table 2. The data span 856 days and have a mean internal uncertainty of 1.2 m s⁻¹.

New radial velocity observations of exoplanetary systems can sometimes result in substantial modification of the best-fit planetary orbits. For example, the two planets in the HD 155358 system were initially reported to be in orbital periods of 195 and 530 days (Cochran et al. 2007). A further five years of observations revealed that the outer planet actually has an orbital period of 391.9 days and is trapped in the 2:1 MMR (Robertson et al. 2012a). In light of the possibility that the best-fit orbits of the two planets may be significantly different than initially presented in Tinney et al. (2006), we begin our orbit fitting process with a wide-ranging search using a genetic algorithm. This approach is often used when the orbital parameters of a planet candidate are highly uncertain (e.g., Tinney et al. 2011; Wittenmyer et al. 2012a; Horner et al. 2012b) or when data are sparse (Wittenmyer et al. 2011). We allowed the genetic algorithm to search a wide parameter space and it ran for 50,000 iterations, testing a total of about 10⁷ possible configurations. We then fit the two datasets simultaneously using GaussFit (Jefferys et al. 1987), a generalized least-squares program used here to solve a Keplerian radial velocity orbit model. The GaussFit model has the ability to allow the offsets between multiple datasets to be a free parameter. The parameters of the best two-planet solution obtained by the genetic algorithm were used as initial inputs to GaussFit and a jitter of 3.3 m s⁻¹ was added in quadrature to the uncertainty of each observation (following Tinney et al. 2006). The best-fit Keplerian solutions are given in Table 3; planetary minimum masses m sin i are derived using a stellar mass of 1.014 ± 0.046 M_☉ (Takeda et al. 2007). This fit has a reduced χ² of 1.63 and a total rms of 6.32 m s⁻¹ (AAT–7.67 m s⁻¹; PFS–2.75 m s⁻¹).

Because the two planets are massive enough and orbit close enough to each other to be interacting, we also apply a full dynamical model to these data. This model includes the effects of planet–planet interactions and can be used to place constraints on the system's inclination to the sky plane, i, a quantity that cannot be determined from Keplerian fitting alone. The system inclination then sets the true masses of the planets. The technique is described fully in Tan et al. (2013) for the HD 82943 two-planet system. The Levenberg–Marquardt (Press et al. 1992) method is adopted as our fitting method. Using the Keplerian best fit as an initial guess, assuming coplanar edge-on orbits, the Levenberg–Marquardt algorithm converges to a local minimum with a $\chi ^2_{\nu }$ of about 1.70 and an rms of about 6.54 m s⁻¹. Based on this local minimum, we conduct a parameter grid search (Lee et al. 2006; Tan et al. 2013) to ensure a global search for the best fit. This minimum is indeed a global dynamical best fit assuming coplanar edge-on orbits; two other local minima with slightly larger $\chi ^2_{\nu }$ have been found. The coplanar edge-on best-fit parameters are listed in Table 4, with their error bars determined by the covariance matrix. This fit and its residuals are shown in the left panel of Figure 1. The right panel of Figure 1 shows that both resonance angles are librating, with amplitudes of ±40° (θ₁) and ±60° (θ₂).

Zoom In Zoom Out Reset image size

Table 4. Dynamical Fit Solutions

Parameter	Planet b	Planet c
K (m s−1)	85.4 ± 2.3	62.3 ± 1.8
Period (days)	189.65 ± 0.21	376.93 ± 0.69
Eccentricity	0.265 ± 0.021	0.198 ± 0.029
ω (deg)	198.3 ± 3.6	294.5 ± 11.3
Mean anomaly (deg)	105.0 ± 5.0	153.4 ± 9.0
a (AU)	0.65 ± 0.01	1.03 ± 0.02
i (deg)	90.0
Mass (MJup)	2.35 ± 0.12	2.19 ± 0.12
$\chi ^{2}_{\nu }$	1.70
rms (m s−1)	6.54
K (m s−1)	83.0 ± 2.1	61.5 ± 1.6
Period (days)	189.01 ± 0.27	379.32 ± 0.92
Eccentricity	0.292 ± 0.022	0.244 ± 0.026
ω (deg)	202.3 ± 3.2	285.3 ± 10.6
Mean anomaly (deg)	102.8 ± 3.9	163.2 ± 8.3
a (AU)	0.65 ± 0.01	1.03 ± 0.02
i (deg)	40.2
Mass (MJup)	3.50 ± 0.17	3.32 ± 0.17
$\chi ^{2}_{\nu }$	1.72
rms (m s−1)	6.59
K (m s−1)	81.4 ± 2.3	63.1 ± 1.6
Period (days)	189.43 ± 0.82	378.29 ± 2.81
Eccentricity	0.308 ± 0.020	0.293 ± 0.021
ω (deg)	205.7 ± 3.4	284.3 ± 9.9
Mean anomaly (deg)	99.5 ± 4.0	165.8 ± 6.8
a (AU)	0.649 ± 0.012	1.029 ± 0.021
i (deg)	20.8
Mass (MJup)	6.22 ± 0.33	6.10 ± 0.31
$\chi ^{2}_{\nu }$	1.80
rms (m s−1)	6.76

Download table as:  ASCIITypeset image

Assuming the planets are in coplanar orbits, we then allow the inclination to the sky plane to vary along with other fitting parameters. Figure 2 shows χ² and rms as a function of sin i and Figure 3 shows best-fit parameters as a function of sin i. The χ² curve is shallow in the range of sin i ≳ 0.6, but then shows a clear local minimum at sin i ∼ 0.36 (i = 208). Two further local minima were found, at inclinations of i = 900 and 402. The parameters of these three solutions are given in Table 4; the planetary masses scale accordingly as 1/sin  i, resulting in more massive planets for the low-inclination solutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The χ² curve and fitting parameters (Figure 3) show discontinuities along the sin i axis, especially those near sin i ∼ 0.43. To understand these discontinuities, we explore grids in different fixed sin i, to see the evolution of the parameter space along different inclinations. Figure 4 shows K₁–K₂ grids for different sin i. Initially, when the orbits are at sin i ∼ 0.43, there is only one χ² minimum (K₁ ∼ 84). Then a new local minimum (K₁ ∼ 82) appears around sin i = 0.427, whose fitting parameters are significantly different from the original minimum (K₁ ∼ 84). When sin i drops down to 0.425, the original minimum vanishes and the new one becomes a single minimum in parameter space. The appearance of the additional $\chi ^2_{\nu }$ minimum results in the big "jump" of fitting parameters at about sin i = 0.43 (see Figure 3).

Zoom In Zoom Out Reset image size

In summary, we have four possible configurations for this system (one Keplerian fit and three dynamical fits). The four solutions are not substantially different from one another, apart from the sin i factor for the three solutions in Table 4, which serves to increase the planetary masses relative to the Keplerian scenario in which we have assumed the planets to be at their minimum masses (m sin i). As a first-order check of dynamical stability, the best-fit system configuration at each inclination was integrated for 10⁸ yr. For all fits so tested, at inclinations of 267, 300, 334, 369, 400, and 900 (sin i = 0.45, 0.5, 0.55, 0.6, 0.64, and 1.0), the systems remained stable for 10⁸ yr. However, since dynamical stability is highly dependent on the initial conditions, we expand on these tests in the next section to obtain a more robust and complete picture of the stability of the various configurations.

When analyzing any multiple-planet system, it is prudent to investigate the long-term dynamical stability of the system. As more multi-planet systems are discovered, the announcement of planetary systems that turn out to be dynamically unfeasible has become increasingly common. Detailed N-body simulations can be used to test the veracity of planet claims. Sometimes, the results of such tests have shown that some systems simply cannot exist in their proposed configuration on astronomically relevant timescales (e.g., Horner et al. 2011; Wittenmyer et al. 2012a; Horner et al. 2012a; Wittenmyer et al. 2013). In other cases, dynamical testing can place additional constraints on planetary systems, particularly when the planets are in or near resonances (e.g., Robertson et al. 2012a, 2012b; Wittenmyer et al. 2012b). In this section, we examine the various solutions for the HD 73526 system, performing detailed dynamical tests of the planetary system configurations as given in Tables 3 and 4. Given that the four solutions are not substantially different from each other in terms of goodness of fit, these dynamical stability tests can serve to discern which scenario is most plausible; a solution that is favored by the fitting process may prove to be unstable, and hence unfeasible.

As in our previous dynamical work (e.g., Marshall et al. 2010; Wittenmyer et al. 2012c; Horner et al. 2012b), we used the Hybrid integrator within the N-body dynamics package Mercury (Chambers 1999) to perform our integrations. We held the initial orbit of the inner planet fixed at its best-fit parameters, as given in Table 3, and then created 126,075 test systems. In those test systems, the initial orbit of the outer planet was varied systematically in semi-major axis a, eccentricity e, periastron argument ω, and mean anomaly M, resulting in a 41 × 41 × 15 × 5 grid of "clones" spaced evenly across the 3σ range in those parameters. We assumed the planets were coplanar with each other and, for the Keplerian case, we assigned masses equivalent to their minimum mass, msin i (Table 3). We then followed the dynamical evolution of each test system for a period of 100 Myr and recorded the times at which either of the planets was removed from the system. Planets were removed if they collided with one another, hit the central body, or reached a barycentric distance of 10 AU.

We performed these dynamical simulations for the Keplerian fit (i = 90°), the dynamical fit at i = 90°, and the lowest-inclination dynamical fit: the configuration given in Table 4 at i = 208. For the last scenario, the planet masses were scaled according to the derived system inclination i. Clearly, the masses of the planets are a proxy for the expected dynamical stability—systems containing more massive planets are likely to be less stable. Hence, the three scenarios we have tested, at i = 90° and i = 208, represent the extremes of dynamical stability (or instability) for the HD 73526 system.

To explore the effects of mutual inclinations between the planets, we performed five additional N-body simulations, for scenarios in which the two planets were inclined with respect to each other. These simulations were set up exactly as described above, using the parameters of the Keplerian solution (Table 3), except at a lower resolution due to computational limitations: a 21 × 21 × 5 × 5 grid in a, e, ω, and M. Five runs were performed, at mutual inclinations of 5°, 15°, 45°, 135°, and 180°. The latter two cases represent retrograde scenarios, which can sometimes allow for a larger range of dynamically stable orbits (Eberle & Cuntz 2010; Horner et al. 2011).

The results of our dynamical stability simulations for the Keplerian solution are shown in Figure 5. We show six panels, for the coplanar and five mutually inclined scenarios as described above. For the coplanar and 5° cases (panels (a) and (b)), the best-fit set of parameters (shown by the open square with 1σ crosshairs) lies in a region of moderate stability, with mean system survival times of ∼10⁶ yr. The stability rapidly degrades as the inclination between the planets becomes significant and even for retrograde cases (panels (e) and (f)), the nominal best-fit system destabilizes within 10⁴ yr. From these simulations, we can conclude that the HD 73526 planets are most likely coplanar with each other. Panels (a) and (b) also show that the stability of the system increases as the outer planet takes on lower eccentricities. For e ≲ 0.2, the mean survival times exceed 10⁷ yr. This is not a surprising result, as high eccentricities generally increase the possibility of strong interactions or even orbit crossings although systems in protected resonances may remain stable for some values of M and ω. Indeed, the statistics of multi-planet systems show that planets in multiple systems tend to have lower eccentricities (Wright et al. 2009; Wittenmyer et al. 2009).

Zoom In Zoom Out Reset image size

As shown in Figure 5, the Keplerian best-fit solution is stable on Myr timescales. However, the colored squares in Figure 5 represent the mean survival times across the range of mean anomalies and ω tested. As the best-fit solution for the HD 73526 system places the planets on resonant orbits, their stability will naturally be highly sensitive to the values of these angles. Hence, Figure 6 shows the outcomes of the 75 individual simulations performed at the best-fit a and e spanning a 5 × 15 grid in M and ω. We see that the nominal solution (where a and e are fixed at the best-fit values) lies at the point of maximum stability. While this is reassuringly consistent with our expectations of enhanced stability within the resonance, we caution that each colored square in Figure 6 represents only a single run and dynamical evolution is known to be a chaotic process (e.g., Horner et al. 2004a, 2004b).

Zoom In Zoom Out Reset image size

The long-term stability results for the dynamical fit with i = 90° are shown in Figure 7. It is immediately apparent that this solution results in a higher degree of stability, with the entire 1σ region exhibiting mean lifetimes exceeding 10⁷ yr. The right panel of Figure 7 shows the results from the 75 individual runs in the central best-fit square, as in Figure 6. For this case, when we examine the dependence on M and ω, we see that the entire diagonal region (including the best fit) is stable for 10⁸ yr. In contrast, Figure 8 shows the dynamical stability of the i = 208 solution from Table 4. Although this fit is formally almost as good as the Keplerian fit, the increased masses implied by the inclination render the system unstable on short timescales (1000 yr).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size