THE PAN-PACIFIC PLANET SEARCH. IV. TWO SUPER-JUPITERS IN A 3:5 RESONANCE ORBITING THE GIANT STAR HD 33844

We have used our iodine-free template spectrum ($R\quad \sim$ 60,000, S/N ∼ 200) to derive spectroscopic stellar parameters. In brief, the iron abundance [Fe/H] was determined from the equivalent widths of 32 unblended Fe lines and the LTE model atmospheres adopted in this work were interpolated from the ODFNEW grid of ATLAS9 (Castelli & Kurucz 2004). The effective temperature (T_eff) and bolometric correction (BC) were derived from the color index B − V and the estimated metallicity using the empirical calibration of Alonso et al. (1999, 2001). Since the color-T_eff method is not extinction-free, we corrected for reddening using $E(B-V)=0.0290$ (Schlegel et al. 1998). The stellar mass and age were estimated from the interpolation of Yonsei-Yale (${{\rm{Y}}}^{2}$) stellar evolution tracks (Yi et al. 2003). The resulting stellar mass of $1.78\pm 0.18\;{M}_{\odot }$ was adopted for calculating the planet masses. Our derived stellar parameters are given in Table 4 and are in excellent agreement with the results of Jones et al. (2011) who found a mass of $1.74\pm 0.18\;{M}_{\odot }$ and a radius of 5.33 ± 0.51 ${R}_{\odot }$.

Table 4.  Stellar Parameters for HD 33844

Parameter	Value	References
Spec. Type	K0 III	Houk & Smith-Moore (1988)
Distance (pc)	100.9 ± 6.5	van Leeuwen (2007)
$(B-V)$	1.040 ± 0.009	Perryman et al. (1997)
$E(B-V)$	0.0290	⋯
AV	0.0903	⋯
Mass (${M}_{\odot }$)	1.78 ± 0.18	This work
	1.74 ± 0.18	Jones et al. (2011)
V sin i (km s−1)	<1	This work
	1.65	Jones et al. (2011)
[Fe/H]	+0.27 ± 0.09	This work
	+0.17 ± 0.10	Jones et al. (2011)
	+0.19 ± 0.12	Luck & Heiter (2007)
Teff (K)	4861 ± 100	This work
	4890	Jones et al. (2011)
	4710	Massarotti et al. (2008)
	4886	Luck & Heiter (2007)
log g	3.24 ± 0.08	This work
	3.05	Jones et al. (2011)
	3.1	Massarotti et al. (2008)
vt (km s−1)	1.00 ± 0.15	This work
	1.17	Jones et al. (2011)
	1.42	Luck & Heiter (2007)
Radius (${R}_{\odot }$)	5.29 ± 0.41	This work
	5.33 ± 0.51	Jones et al. (2011)
Luminosity (${{\rm{L}}}_{\odot }$)	14.1 ± 1.8	This work
	14.4	Jones et al. (2011)
	12.6	Massarotti et al. (2008)
Age (Gyr)	${1.88}_{-0.48}^{+0.76}$	This work

Download table as:  ASCIITypeset image

Particularly for giant stars, where spots and pulsations can induce spurious radial velocity shifts with periods of hundreds of days, any claim of orbiting planets must be carefully examined to rule out intrinsic stellar signals (e.g., Hatzes & Cochran 2000; Reffert et al. 2015; Trifonov et al. 2015). For HD 33844, the periods of the two signals (551 and 916 days) are nowhere near the window function peaks at 384 and 8.1 days (AAT) or 30 and 364 days (Keck). Spurious periods in observational data commonly arise at those periods due to sampling (imposed by bright-time scheduling and the yearly observability of a given target).

To check whether the observed velocity variations could be due to intrinsic stellar processes, we examined the All-Sky Automated Survey (ASAS) V band photometric data for HD 33844 (Pojmanski & Maciejewski 2004). A total of 596 epochs were obtained from the ASAS All Star Catalog.¹⁴ We computed the mean magnitude per epoch over the five apertures then subjected the time series to an iterative sigma-clipping process. We removed points more than $3\sigma$ from the grand mean then recalculated the mean and its standard deviation. This process was performed three times, after which 511 epochs remained with a mean value of 7.28 ± 0.02. The generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) is shown in Figure 3 with the periods of the planets marked as dashed lines. While there are significant periodicities at 730 and 1250 days, there is little power near the periods of the candidate planets (551 and 916 days).

Zoom In Zoom Out Reset image size

We also checked for correlation between the radial velocities and the equivalent width of the Hα absorption line, which has been used as an activity indicator for giants (Hatzes et al. 2015) as well as for M dwarfs (Robertson et al. 2013). The equivalent widths were measured in a 2 Å window centered on Hα to avoid contamination by telluric lines. A generalized Lomb–Scargle periodogram of the Hα equivalent widths from the AAT spectra (Figure 4) shows no significant periodicities and there are no correlations with the velocities. Furthermore, the bisector velocity spans (defined as the velocity difference between line bisectors from the upper and lower part of an absorption line) computed from the FEROS spectra show no correlation with the radial velocities.

Zoom In Zoom Out Reset image size

The HD 33844 system appears to contain two super-Jupiter planets in orbits relatively close to each other. Given their mass and proximity, it is clearly important to consider whether the planets are dynamically feasible. That is, could planets on such tightly packed orbits be dynamically stable on timescales comparable to the lifetime of the system? A first estimate of the system's stability can be garnered by simply assessing the dynamical separation of the two planets, considering the separation of their orbits compared to their mutual Hill radius. Following Gladman (1993), we can calculate the mutual Hill radius of the two planets as follows:

$$\begin{eqnarray}&&{R}_{{\rm{H}}}={\left[\displaystyle \frac{({m}_{1}+{m}_{2})}{3{M}_{\odot }}\right]}^{1/3}\left[\displaystyle \frac{({a}_{1}+{a}_{2})}{2}\right],\end{eqnarray} \tag{ 1 }$$

where the symbols have their usual meaning and the subscripts refer to the inner (1) and outer (2) planets respectively. Following this formulism, we find that the best-fit orbits for the two candidate planets are separated by 3.8 times their mutual Hill radius (R_H = 0.167 AU). For low-eccentricity orbits Gladman (1993) found that orbits became unstable at separations smaller than $\sim 2\sqrt{3}=3.46\;{R}_{{\rm{H}}}$. The HD 33844 system is therefore close to this critical separation and as the proposed orbits are somewhat eccentric, it is clearly important to subject the proposed planets to further scrutiny. In contrast to widely separated systems such as HD 121056, (where the planets orbit at 0.4 and 3.0 AU—more than 9 mutual Hill radii apart), for which N-body simulations were not necessary, here we must rigorously test the HD 33844 system stability as in our previous work (e.g., Marshall et al. 2010; Horner et al. 2013; Wittenmyer et al. 2013a).

Most interesting are those systems (e.g., Robertson et al. 2012a, 2012b; Wittenmyer et al. 2012b) for which the planets prove stable and dynamically feasible across just a small fraction of the potential orbital solutions. In these cases, which typically feature planets moving on or close to mutually resonant orbits, dynamical simulations serve a dual purpose. First, they provide evidence that supports the existence of the planets and secondly, they provide a strong additional constraint on the potential orbits followed by those planets, helping us to better tie down their true orbits than can be achieved on the basis of the observations alone.

Here, we study the stability of the candidate planets orbiting HD 33844 following a now well-established route. We created a suite of 126,075 copies of the HD 33844 system. In each of these cloned systems, the initial orbit of HD 33844b was the same, located at its nominal best-fit values (Table 5). For each system, we systematically varied the initial semimajor axis (a), eccentricity (e), argument of periastron (ω), and mean anomaly (M) of HD 33844c. The masses of the planets were held fixed at their minimum values (m sin i, Table 5). We note that changing the mass of the planets could alter the stability of the system. This can be illustrated by examination of Equation (1): it is immediately apparent that if the masses of the planets are increased, so too is the size of their mutual Hill radius and thereby the strength of their mutual interaction. However, the effect is actually relatively small when compared with the influence of their orbital elements. As such, in this work we solely explore the influence of "element space" and leave the exploration of "mass space" for future work, once the orbital elements of the planets have been better constrained through followup observations. Once the uncertainties in those elements are sufficiently small, it might be possible to use "mass space" to constrain the maximum masses of the planets and thereby obtain some constraints on the inclination of the system to our line of sight, but such calculations are beyond the scope of this work. We do, however, note that a tentative upper limit on the masses of the planets can be obtained using the resonance overlap criterion (e.g., Wisdom 1980; Deck et al. 2013). This analytic estimate sets an upper bound of ∼10 ${M}_{{\rm{Jup}}}$ for each planet. As such, we can be fairly confident that the two bodies are planetary in nature rather than being brown dwarfs.

Since previous studies have shown that the stability of a system is most strongly dependent on semimajor axis and eccentricity, we tested 41 discrete values of each of these variables spanning the full $\pm 3\sigma$ uncertainty ranges. At each of the 1681 a − e pairs created in this way, we tested 15 unique values of the argument of periastron and five of the mean anomaly, distributed in each case evenly across the $1\sigma$ uncertainty ranges in these variables. This gave us a total of 126,075 unique potential orbits for HD 33844c.

We then used the Hybrid integrator within the n-body dynamics package Mercury (Chambers 1999) to follow the evolution of each of the test planetary systems for a period of 100 Myr. Simulations were stopped early if either of the planets were ejected from the system (upon reaching a barycentric distance of 5 AU, which would require significant strong instability between the two planets). They were also halted if either of the planets fell into the central star whose mass was set at $1.78\;{M}_{\odot }$ or if they collided with one another. If any of these events happened, the time of collision/ejection was recorded and the simulation was brought to a close.

As a result of these simulations, we are able to examine the dynamical stability of the HD 33844 system as a function of the initial orbit on which HD 33844c was placed. Figure 5 shows the stability of the system as a function of the initial semimajor axis and eccentricity of that planet's orbit. In Figure 6 we show how the orbital solutions tested in our dynamical simulations fit to the observed data, expressed in terms of the difference in total ${\chi }^{2}$ relative to the best fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is immediately clear from Figure 5 that the proposed orbital solution for the system lies in a region of complex dynamical behavior with both extremely stable and unstable solutions being possible. It is reassuring, however, to note that broad regions of dynamical stability lie comfortably within the 1-σ uncertainties on the proposed solution, particularly toward lower eccentricities. We note that least-squares radial velocity fitting routines are well-known to inflate eccentricities (e.g., Shen & Turner 2008; O'Toole et al. 2009; Wittenmyer et al. 2013c). That stability is due to orbits in that region being trapped in mutual 3:5 mean motion resonance with the orbit of HD 33844b. From Table 5 the Keplerian solution gives a period ratio of 1.661 and 1.687 for the dynamical fit. Within their uncertainties these solutions agree with each other and are wholly consistent with the 3:5 resonance (period ratio 1.667).

In addition to the central region of stability, orbits at smaller semimajor axes fall into a broad region of stability that extends across the full span of tested orbital eccentricities. This feature is the result of the mutual 2:3 mean motion resonance between HD 33844c and HD 33844b which is centered on 2.102 AU.

We can also see evidence of unstable resonant behavior through the plot. Most strikingly, there is a band of unstable solutions centered at 2.33 AU. This band is the result of the 4:7 mean motion resonance between the two planets. A further unstable resonant region can be seen around 2.195 AU, associated with the 5:8 resonance between the planets. Finally, the 5:9 resonance can be found at 2.375 AU, which is likely the cause of the sculpting of the stability of the system in that region.

As such, we can conclude that the candidate planets orbiting HD 33844 are dynamically feasible but that they most likely move on mutually resonant, low-eccentricity orbits. As a further check, we investigated the behavior of the resonant angles for a number of key resonances in the vicinity of the complicated "stability terrain" around the best-fit orbit. We found that the best-fit orbital solution is strongly influenced by its proximity to the 3:5 resonance. In particular, we found that the resonant angle $\phi =5{\lambda }_{2}-3{\lambda }_{1}-{\omega }_{1}+{\omega }_{2}$ alternates between libration and circulation in a regular manner, completing three full cycles (one libration + one circulation) per 500 years of integration. A similar but noisier behavior was observed for a scenario in which the planets were located on orbits whose periods were in 5:8 commensurability. Here the resonant angle for the 5:8 mean motion resonance switched chaotically between periods of libration (lasting up to 1000 years) and period of smooth, slow circulation. The influence of resonant interaction for these solutions was unmistakable, and given the proximity of the best-fit solution to the location of the 3:5 resonance, it seems most likely that the planets are trapped within it (although we note that the 5:8 and 7:12 mean motion resonances also fall within the 1σ uncertainty in semimajor axis) together with an abundance of higher-order weaker resonances.