A PAIR OF GIANT PLANETS AROUND THE EVOLVED INTERMEDIATE-MASS STAR HD 47366: MULTIPLE CIRCULAR ORBITS OR A MUTUALLY RETROGRADE CONFIGURATION

Evolved stars (giants and subgiants) have been extensively searched for the last 10 years, aiming to explore planets around more massive stars than the Sun and to investigate the orbital evolution of planetary systems during the post-main-sequence phase. Precise radial-velocity surveys have discovered about 120 substellar companions around such evolved stars¹³ so far, and they are now known to exhibit statistical properties that are not necessarily similar to those around solar-type stars (e.g., Hatzes et al. 2005, 2006; Setiawan et al. 2005; Lovis & Mayor 2007; de Medeiros et al. 2009; Döllinger et al. 2009; Johnson et al. 2011; Wittenmyer et al. 2011; Omiya et al. 2012; Wang et al. 2012, 2014; Novak et al. 2013; Sato et al. 2013b; Lee et al. 2014; Trifonov et al. 2014; Jones et al. 2015; Reffert et al. 2015). Recently, very-high-precision photometry by the Kepler space telescope has succeeded in detecting planetary transits on evolved stars (Batalha et al. 2013). The discoveries include short-period planets and sub-Jupiter-mass ones around giants (Huber et al. 2013; Ciceri et al. 2015; Lillo-Box et al. 2014; Ortiz et al. 2015; Quinn et al. 2015; Sato et al. 2015), which had rarely been found by the ground-based radial-velocity surveys, and also candidates with longer periods, which could be counterparts of those found by the radial-velocity surveys.¹⁴

Among the planetary systems around evolved stars, multiple-planet systems are of particular interest in terms of the formation and evolution of planetary systems. Figure 1 shows the cumulative number of planet pairs in multiple-planet systems discovered by ground-based, radial-velocity, and transit surveys. The total mass of the planet pair, as well as the $\mathrm{log}\quad g$ and mass of the host star, are indicated by the color and size of the symbol, respectively. Interestingly, almost all of the giant-planet pairs with period ratios smaller than two are preferentially found around evolved intermediate-mass stars.¹⁵ For example, 24 Sex (G5 IV; Johnson et al. 2011) and η Cet (K2 III; Trifonov et al. 2014) host double giant-planet systems probably in mean-motion resonances of 2:1, which suggests that these systems have experienced differential convergent orbital migration (e.g., Lee & Peale 2002; Kley et al. 2004). In contrast, double giant-planet systems around HD 200964 (K0 IV; Johnson et al. 2011) and HD 5319 (G5 IV; Robinson et al. 2007; Giguere et al. 2015) are near the 4:3 resonance, which are considered to be difficult to form either via convergent migration, scattering, or in situ formation (Rein et al. 2012). BD+20 2457 (K2 III) hosts two brown dwarfs that are suggested to be near 3:2 resonance (Niedzielski et al. 2009). However, a detailed orbital stability analysis of the system revealed that the best-fit orbits derived by radial-velocity data are unstable, though they would be stable if the two brown dwarfs were mutually in retrograde configuration (Horner et al. 2014). Giant-planet pairs with such small period ratios are rarely found around solar-type stars, though many less massive ones including Neptunes and Super-Earths are found around them.

Zoom In Zoom Out Reset image size

Here we report the detection of a double giant-planet system around the intermediate-mass giant HD 47366 (K1 III, $M=1.8\quad {M}_{\odot }$), which adds to the growing population of multi-planet systems with small period ratios around evolved stars. The discovery was made by Okayama and Xinglong Planet Search Programs (e.g., Wang et al. 2012; Sato et al. 2013b) facilitated by joint precise radial-velocity observations from the Australian Astronomical Observatory (Sato et al. 2013a). The planetary system is intriguing for several reasons: the best-fit Keplerian orbit is unstable, it is near but less likely in 2:1 mean-motion resonance, and could be stable if the orbits are nearly circular or in retrograde configuration.

The rest of the paper is organized as follows. The stellar properties are presented in Section 2 and the observations are described in Section 3. The orbital solution and results of dynamical stability analysis are presented in Sections 4 and 5, respectively. Section 6 is devoted to discussion and summary.

HD 47366 (HIP 31674, BD−12 1566, HR 2437, TYC 5373-2001-1) is listed in the Hipparcos Catalogue (ESA 1997) as a K1 III star, with a visual magnitude of V = 6.11 and a color index of B − V = 0.994. The Hipparcos parallax π = 12.5 ± 0.42 mas (van Leeuwen 2007) corresponds to a distance of 80.0 ± 2.7 pc. The reddening E(B − V) = 0.028 was obtained from the Galactic dust map of Schlegel et al. (1998), with the correction given by Bonifacio et al. (2000) and a scaling factor of $1-\mathrm{exp}\;(-| d\mathrm{sin}b| /h)$, where d is the distance, b is the Galactic latitude, and h = 125 pc is the scale height of the reddening layer. The absolute magnitude M_V = 1.51 was derived from the distance and the interstellar extinction A_V = 3.1E(B − V). By adopting the broadband photometric color B − V and the estimated metallicity with the empirical calibration relation of Alonso et al. (1999, 2001), we derived the bolometric correction B.C. = −0.309 and the effective temperature T_eff = 4866 ± 100 K. We used a high signal-to-noise ratio (S/N ∼ 200), iodine-free spectrum taken with HRS to measure the equivalent widths (EWs) of ∼30 Fe i lines, to derive the iron abundance [Fe/H]. The line lists as well as their oscillation strengths ($\mathrm{log}{gf}$) were mainly taken from Hekker & Meléndez (2007), in which iron lines were carefully selected to avoid any blend by atomic or CN lines. The model atmosphere used in this work was interpolated from the line-blanketed, local thermodynamic equilibrium (LTE) ATLAS9-ODFNEW grid (Castelli & Kurucz 2004). The microturbulent velocity v_t was obtained by minimizing the trend between the abundances of different Fe i lines and their reduced equivalent widths ($\mathrm{log}({\rm{EWs}}/\lambda )$). The macroturbulent velocity was estimated with the empirical relations of v_macro versus T_eff given by Hekker & Meléndez (2007), and the projected rotational velocity ($v\mathrm{sin}i$) was determined with the method of Fekel (1997). The stellar mass, surface gravity ($\mathrm{log}g$), radius, and age were derived using a Bayesian approach with the Geneva database (Lejeune & Schaerer 2001), which includes the post-helium flash phases for stars with M ≥ 1.7 M_⊙. First, we interpolated an extensive grid of evolutionary tracks, with ΔM = 0.05 within 1.2 ≤ M/M_⊙ ≤ 3.6, Δ[Fe/H] = 0.02 within −0.4 ≤ [Fe/H] ≤ +0.3, and 500 points in each track. Then, for each data point, we calculated the likelihood functions of $\mathrm{log}L$, T_eff and [Fe/H] to match the observed values by assuming Gaussian errors. We adopted uniform prior probabilities of mass and [Fe/H]. It is noted that stars evolving more slowly have a higher probability of being observed. Without correcting this evolution effect, the resulting parameters would bias toward the rapid evolution phases. We therefore weighted the probability of each point along its evolutionary track by the normalized time-step (a_i+1,j − a_i,j)/(a_n,j − a_a,j), where a_i,j is the age of the ith interpolated point in the jth track, and n = 500 is the number of interpolated points in each track. Eventually, the probability distribution functions (PDFs) of the parameters yield M = 1.81 ± 0.13 M_⊙, R = 7.30 ± 0.33 R_⊙, $\mathrm{log}g\quad =2.97\pm 0.06$, and age = 1.61 ± 0.53 Gyr. The stellar mass is particularly important to derive the minimum masses of the orbiting planets detected with the Doppler technique. However, previous spectroscopic analyses gave discrepant results (1.87 M_⊙ by Liu et al. 2010; 1.2 M_⊙ by Mishenina et al. 2006) for HD 47366, which may be due to the different methods on finding T_eff and $\mathrm{log}g$. Our determinations were based on a similar method to that of Liu et al. (2010), but used the Geneva evolutionary tracks, instead of the Y² model (Yi et al. 2003), which does not account for the evolutionary phases after the helium-core flash. We found that the probability that the star has passed through the RGB tip and in-core helium burning phase is ∼88%. The stellar parameters of HD 47366 are listed in Table 1. In Figure 2, we plotted HD 47366 on the H-R diagram, together with the evolutionary tracks for stars with different masses and metal contents.

Zoom In Zoom Out Reset image size

Table 1.  Stellar Parameters for HD 47366

Parameter	This Work	Liu et al. (2010)	Mishenina et al. (2006)
Spectral Type	K1 III	⋯	⋯
Hipparcos Parallax π (mas)	12.5 ± 0.42	⋯	⋯
Distance (pc)	80.0 ± 2.7	⋯	⋯
Visual Magnitude V	6.11	⋯	⋯
Color Index B − V	0.994 ± 0.002	⋯	⋯
Redenning $E(B-V)$	0.028	⋯	⋯
Interstellar Extinction AV	0.087	⋯	⋯
Absolute Magnitude MV	1.51	1.467	1.044
Bolometric Correction B.C.	−0.309	−0.302	⋯
Bolometric Magnitude Mbol	1.20	1.165	⋯
Effective Temperature Teff (K)	4866 ± 100	4883 ± 100	4772
Surface Gravity $\mathrm{log}g$ (cgs)	2.97 ± 0.06	3.00 ± 0.1	2.60
Metallicity [Fe/H]	−0.02 ± 0.09	−0.10 ± 0.1	−0.16
Microturbulent Velocity vt (km s−1)	1.24	1.2 ± 0.2	1.2
Macroturbulent Velocity vmacro (km s−1)	5.54 ± 0.45	⋯	⋯
Projected Rotational Velocity $v\mathrm{sin}i$ (km s−1)	4.3 ± 0.8	⋯	⋯
Luminosity L (L⊙)	26.1 ± 1.8	26.2	38.5
Radius R (R⊙)	7.30 ± 0.33	7.16	9.09
Mass M (M⊙)	1.81 ± 0.13	1.87	1.2
Age (Gyr)	1.61 ± 0.53	⋯	⋯

Download table as:  ASCIITypeset image

The star shows no significant emission in the core of Ca ii HK lines, as shown in Figure 3, which suggests that the star is chromospherically inactive. We did not analyze flux variations of Ca ii HK lines for the star because of the poor S/N of the core of the lines.

Zoom In Zoom Out Reset image size

3.1. OAO Observations

Observations of HD 47366 at OAO were made with the 1.88 m reflector and the HIgh Dispersion Echelle Spectrograph (HIDES; Izumiura 1999) from 2006 December to 2014 April. We used both the conventional slit mode (HIDES-S) and the high-efficiency fiber-link mode (HIDES-F) of the spectrograph that became available in 2010 (Kambe et al. 2013). In the case of the HIDES-S, a slit width of the spectrograph was set to 200 μm (076) corresponding to a spectral resolution (R = λ/Δλ) of 67,000 by about a 3.3-pixel sampling. In the case of HIDES-F, a width of the image sliced by an image slicer is 105, corresponding to a spectral resolution of R = 55,000 by a 3.8-pixel sampling. Each observing mode uses its own iodine absorption cell for precise radial-velocity measurements (I₂ cell; Kambe et al. 2002, 2013), which provide a fiducial wavelength reference in a wavelength range of 5000–5800 Å. We have obtained a total of 50 and 7 data points of HD 47366 using HIDES-S and HIDES-F mode, respectively. The reduction of echelle data (i.e., bias subtraction, flat-fielding, scattered-light subtraction, and spectrum extraction) was performed using the IRAF¹⁶ software package in the standard way.

We performed a radial-velocity analysis following the method of Sato et al. (2002, 2012) and Butler et al. (1996), in which an I₂-superposed spectrum is modeled as a product of a high-resolution I₂ and a stellar template spectrum convolved with a modeled point-spread function (PSF) of the spectrograph. We model the PSF using multiple Gaussian profiles (Valenti et al. 1995) and obtain the stellar spectrum by deconvolving a pure stellar spectrum with the spectrograph PSF estimated from an I₂-superposed B-star spectrum. We obtained a typical measurement error in radial-velocity of about 4–5 m s⁻¹ for the star, which was estimated from an ensemble of velocities from each of the ∼300 spectral chunks (each ∼3 Å long) in every exposure. We listed the derived radial velocities for OAO data in Table 2 together with the estimated uncertainties.

Table 2.  OAO Radial Velocities for HD 47366

JD−2450000	Velocity (${\text{m s}}^{-1}$)	Uncertainty (${\text{m s}}^{-1}$)	Mode
4093.18674	11.36	4.59	HIDES-S
4144.11514	−19.66	4.69	HIDES-S
4172.01541	−18.29	3.85	HIDES-S
4467.12999	−4.46	4.31	HIDES-S
4553.97104	−31.25	4.14	HIDES-S
4755.29094	32.92	3.46	HIDES-S
4796.25711	36.33	4.17	HIDES-S
4818.17076	19.49	3.87	HIDES-S
4823.11554	14.96	4.17	HIDES-S
4834.16358	12.38	4.36	HIDES-S
4856.08725	−6.42	3.92	HIDES-S
4863.09739	−6.67	4.41	HIDES-S
4864.03798	−11.43	3.60	HIDES-S
4883.04016	−7.18	4.55	HIDES-S
4923.93523	−17.64	3.73	HIDES-S
5108.31507	40.77	4.18	HIDES-S
5131.32409	47.90	3.57	HIDES-S
5161.34584	39.99	4.64	HIDES-S
5182.16055	11.19	4.10	HIDES-S
5204.11564	−7.29	4.29	HIDES-S
5232.99361	−43.84	3.62	HIDES-S
5267.01163	−60.09	3.57	HIDES-S
5294.93562	−62.39	3.54	HIDES-S
5471.29879	24.76	3.61	HIDES-S
5502.23739	33.89	4.15	HIDES-S
5525.22640	33.98	4.85	HIDES-S
5545.21850	18.86	3.73	HIDES-S
5581.04001	−12.92	3.75	HIDES-S
5611.01311	−8.59	4.13	HIDES-S
5625.02350	−13.62	3.49	HIDES-S
5656.93188	−10.50	3.74	HIDES-S
5854.27845	33.84	3.69	HIDES-S
5877.25348	18.63	5.42	HIDES-S
5882.29779	12.23	3.83	HIDES-S
5922.13357	−12.73	3.88	HIDES-S
5938.09758	−34.44	3.79	HIDES-S
5976.03887	−49.42	4.57	HIDES-S
6010.93116	−64.77	5.74	HIDES-S
6032.96687	−55.43	4.48	HIDES-S
6215.30762	27.47	4.04	HIDES-S
6235.27265	41.83	3.68	HIDES-S
6249.30138	27.36	4.69	HIDES-S
6284.08997	−0.00	4.69	HIDES-S
6552.29701	42.87	5.53	HIDES-S
6577.28414	4.35	4.84	HIDES-S
6578.28739	14.21	4.48	HIDES-S
6578.30265	14.55	4.27	HIDES-S
6579.25810	12.00	5.63	HIDES-S
6618.19944	4.11	3.91	HIDES-S
6753.98409	−63.39	4.45	HIDES-S
6609.21551	21.82	4.04	HIDES-F
6616.27916	33.33	3.44	HIDES-F
6664.09442	−1.41	4.02	HIDES-F
6672.05938	−16.16	3.85	HIDES-F
6697.15561	−9.17	4.70	HIDES-F
6713.98988	−17.94	4.87	HIDES-F
6727.00570	−12.14	3.47	HIDES-F

Download table as:  ASCIITypeset image

3.2. Xinglong Observations

Observations of HD 47366 at the Xinglong Station started in 2006 November using the 2.16 m reflector and the Coudé Echelle Spectrograph (CES; Zhao & Li 2001). We used the blue-arm, middle focal length camera, and 1 K × 1 K CCD (pixel size of 24 × 24 μm²; hereafter CES-O) of the spectrograph, which covered a wavelength range of 3900–7260 Å with a spectral resolution of 40,000 by 2-pixel sampling. Although a wide wavelength range was obtained with a single exposure, only a wave band of Δλ ∼ 470 Å was available for radial-velocity measurements using an I₂ cell because of the small format of the CCD. Radial-velocity analysis for CES data was performed by the optimized method of Sato et al. (2002) for CES (Liu et al. 2008). We took five data points for HD 47366 using CES-O from 2006 November to 2007 December whose radial-velocity uncertainties were better than 25 m s⁻¹. In 2009 March, the old CCD was replaced by a new CCD having a smaller pixel size (2 K × 2 K, 13 × 13 μm²), and it gave a better radial-velocity precision than before, though the wavelength coverage was unchanged (Wang et al. 2012). We made observations of HD 47366 using the new CCD (hereafter CES-N) from 2009 December to 2010 December, and collected a total of 26 data points with a radial-velocity precision of 12–25 m s⁻¹ depending on weather conditions.

Since 2011 November, we have observed HD 47366 with the newly developed High Resolution Spectrograph (HRS) attached at the Cassegrain focus of the 2.16 m telescope. The fiber-fed spectrograph is the successor of the CES, giving a higher wavelength resolution and throughput. The single 4 K × 4 K CCD provides a simultaneous wavelength coverage of 3700–9200 Å, and the slit width of 190 μm gives a wavelength resolution of R = 45,000 with a 3.2-pixel sampling. An iodine cell is installed before the fiber entrance at the Cassegrain focus. We collected a total of 60 observations for the star from 2011 November to 2014 March. Radial-velocity analysis for the HRS data was performed by the same method as Sato et al. (2002, 2012), but optimized for HRS. The resulting measurement errors in radial velocity are ranging from 6 to 20 m s⁻¹ depending on weather conditions. The derived radial velocities for Xinglong data are listed in Table 3 together with the estimated uncertainties.

Table 3.  Xinglong Radial Velocities for HD 47366

JD−2450000	Velocity (${\text{m s}}^{-1}$)	Uncertainty (${\text{m s}}^{-1}$)	Mode
4046.34259	53.99	27.92	CES-O
4133.12441	53.23	21.39	CES-O
4398.37834	101.38	17.20	CES-O
4459.23135	117.08	30.81	CES-O
4459.24962	32.80	27.00	CES-O
5172.26505	96.73	19.05	CES-N
5172.28626	76.60	15.71	CES-N
5172.30759	51.67	14.37	CES-N
5173.27944	45.82	15.67	CES-N
5173.30065	56.52	16.41	CES-N
5173.32198	65.67	17.79	CES-N
5173.34310	65.70	14.76	CES-N
5226.05844	−45.55	24.43	CES-N
5226.10875	−20.17	22.59	CES-N
5227.08807	−26.57	11.76	CES-N
5227.10979	−46.69	13.18	CES-N
5227.13103	−12.38	14.99	CES-N
5281.99890	−14.95	13.80	CES-N
5282.00964	−22.57	15.56	CES-N
5282.02035	−61.02	17.25	CES-N
5282.98558	−45.16	16.54	CES-N
5283.00671	−59.43	16.13	CES-N
5283.02782	−52.03	15.03	CES-N
5283.98000	−51.62	13.79	CES-N
5284.00528	−29.21	20.65	CES-N
5284.02646	−48.95	20.95	CES-N
5284.04756	−6.26	21.05	CES-N
5520.33128	39.59	17.73	CES-N
5520.35299	21.56	16.20	CES-N
5552.23941	34.69	14.87	CES-N
5552.26341	24.93	21.94	CES-N
5878.34042	43.36	13.50	HRS
5878.36563	33.75	13.50	HRS
5878.38866	−0.73	16.05	HRS
5878.41177	42.81	16.18	HRS
5879.29416	4.27	19.96	HRS
5905.20117	23.16	11.56	HRS
5905.22422	14.68	12.04	HRS
5905.24720	36.90	14.27	HRS
5905.27024	48.98	15.64	HRS
5907.26067	1.38	8.18	HRS
5907.28402	−8.44	9.30	HRS
5907.30698	8.20	10.15	HRS
5907.32995	7.82	9.57	HRS
5936.13872	−29.36	10.11	HRS
5936.16225	−35.52	10.34	HRS
5936.18572	−32.01	10.67	HRS
5936.20863	−22.04	11.92	HRS
5959.11008	−30.57	22.97	HRS
5959.13303	−43.36	20.46	HRS
5959.15600	−43.22	18.84	HRS
5959.17896	−62.76	24.07	HRS
5960.09450	−49.14	13.67	HRS
5960.11747	−39.46	15.51	HRS
5966.04907	−39.47	19.32	HRS
5966.07211	−31.79	20.95	HRS
5966.09510	−79.17	17.80	HRS
5966.11809	−67.76	17.67	HRS
5997.97169	−66.24	10.18	HRS
5997.99495	−52.05	11.16	HRS
5998.01824	−47.66	11.05	HRS
5998.04152	−60.49	10.68	HRS
6199.33390	51.35	6.54	HRS
6199.35708	45.59	7.80	HRS
6232.33343	54.89	5.89	HRS
6232.35667	61.48	6.12	HRS
6287.18192	30.60	9.50	HRS
6287.20515	35.52	10.34	HRS
6287.22838	32.06	9.39	HRS
6287.25312	47.44	11.03	HRS
6318.08350	31.16	11.25	HRS
6318.10670	13.27	11.99	HRS
6318.13045	21.01	12.00	HRS
6318.15375	25.07	11.17	HRS
6344.02830	11.76	8.52	HRS
6344.05152	23.29	9.47	HRS
6344.07471	22.34	9.87	HRS
6583.31389	22.42	10.63	HRS
6583.33707	35.39	11.03	HRS
6583.36032	20.96	11.67	HRS
6611.28853	18.21	7.27	HRS
6611.31171	20.81	8.39	HRS
6611.33507	7.72	8.55	HRS
6645.18797	−18.03	12.37	HRS
6645.21133	9.17	12.68	HRS
6645.23463	−12.05	15.41	HRS
6645.25792	−22.92	16.31	HRS
6701.05113	−50.92	11.31	HRS
6701.07444	−55.62	15.95	HRS
6728.98753	−39.22	16.13	HRS
6729.01088	−50.79	19.71	HRS

Download table as:  ASCIITypeset images: 1 2

3.3. AAT Observations

Since the inner planet has a period near one year, and HD 47366 is quite far south for OAO, we obtained confirmation observations from the 3.9-m Anglo-Australian Telescope (AAT). This multi-telescope approach was also critical for the confirmation of the 360-day planet HD 4732b (Sato et al. 2013a). The UCLES echelle spectrograph (Diego et al. 1990) delivers a resolving power of R ∼ 45,000 with the 1-arcsecond slit, and has been used for 16 years by the Anglo-Australian Planet Search (e.g., Butler et al. 2001; Tinney et al. 2001; Jones et al. 2010; Wittenmyer et al. 2012b, 2014). Calibration of the spectrograph point-spread function is achieved using an iodine absorption cell temperature-controlled at 60.0 ± 0.1°C. The iodine cell superimposes a forest of narrow absorption lines from 5000 to 6200 Å, allowing for the simultaneous calibration of instrumental drifts as well as a precise wavelength reference (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 epoch spectra observed through the iodine cell. HD 47366 was observed at 13 epochs between 2013 March and 2014 April; the AAT velocities are given in Table 4.

Table 4.  AAT Radial Velocities for HD 47366

JD−2450000	Velocity (${\text{m s}}^{-1}$)	Uncertainty(${\text{m s}}^{-1}$)
6374.88737	−4.3	1.2
6375.91758	1.9	1.1
6376.93270	0.0	1.1
6377.98528	−2.0	1.2
6399.91709	13.0	1.4
6526.29116	70.0	1.4
6529.31433	85.4	1.5
6555.29760	63.2	1.5
6556.21605	67.5	1.5
6686.00004	−17.2	1.6
6744.92245	−42.0	1.5
6747.92675	−46.0	1.6
6764.87970	−40.9	1.3

Download table as:  ASCIITypeset image

The collected radial-velocity data for HD 47366 obviously show a variation with a period near one year. Although the one- year periodicity and the relatively low declination of the star allowed us to cover only half of the cycle, we can clearly see the periodic variation of the velocity amplitude every year, indicating the existence of a second period near two years. Thus we performed the least-squares orbital fitting by a double Keplerian model. The orbital parameters and the uncertainties were derived using the Bayesian Markov Chain Monte Carlo (MCMC) method (e.g., Ford 2005; Gregory 2005; Ford & Gregory 2007), following the analysis in Sato et al. (2013b). We took account of velocity offsets ΔRV of HIDES-F, CES-O, CES-N, HRS, and AAT data relative to HIDES-S data as free parameters in the orbital fitting. Extra Gaussian noises s for each of the six data sets were also incorporated as free parameters, though those for CES-O (s₃) and HRS (s₅) are fixed to 0 since their rms scatters to the best-fit orbit are comparable to their measurement errors. We generated five independent chains having 10⁷ points with acceptance rates of about 25%, the first 10% of which were discarded, and confirmed that each parameter was sufficiently converged based on the Gelman-Rubbin statistic (Gelman & Rubin 1992). We derived the median value of the merged posterior PDF for each parameter and set a 1σ uncertainty as the range between 15.87% and 84.13% of the PDF.

In Figures 4 and 5, we plot the derived Keplerian orbits together with the radial-velocity points and their measurement errors including the extra Gaussian noises. The inner planet (planet b) has orbital parameters of a period of ${P}_{b}={363.3}_{-2.4}^{+2.5}$ days, an eccentricity of ${e}_{b}={0.089}_{-0.060}^{+0.079}$, a minimum mass  of ${m}_{b}\mathrm{sin}i={1.75}_{-0.17}^{+0.20}\quad {M}_{{\rm{J}}}$, and a semimajor axis of ${a}_{b}={1.214}_{-0.029}^{+0.030}$ au, and the outer planet (planet c) has ${P}_{c}={684.7}_{-4.9}^{+5.0}$ days, ${e}_{c}={0.278}_{-0.094}^{+0.067}$, ${m}_{c}\mathrm{sin}i={1.86}_{-0.15}^{+0.16}\quad {M}_{{\rm{J}}}$, and ${a}_{c}={1.853}_{-0.045}^{+0.045}$ au. The obtained parameters are listed in Table 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Orbital Parameters for HD 47366

Parameter	HD 47366 b	HD 47366 c
Period P (days)	${363.3}_{-2.4}^{+2.5}$	${684.7}_{-4.9}^{+5.0}$
RV semiamplitude K1 (${\text{m s}}^{-1}$)	${33.6}_{-2.8}^{+3.6}$	${30.1}_{-2.0}^{+2.1}$
Eccentricity e	${0.089}_{-0.060}^{+0.079}$	${0.278}_{-0.094}^{+0.067}$
Longitude of periastron ω (deg)	${100}_{-71}^{+100}$	${132}_{-20}^{+17}$
Periastron passage Tp (JD−2450000)	${122}_{-55}^{+71}$	${445}_{-62}^{+55}$
Minimum mass ${m}_{2}\mathrm{sin}i$ (${M}_{{\rm{J}}}$)	${1.75}_{-0.17}^{+0.20}$	${1.86}_{-0.15}^{+0.16}$
Semimajor axis a (au)	${1.214}_{-0.029}^{+0.030}$	${1.853}_{-0.045}^{+0.045}$
Semimajor axis ratio ab/ac	0.6555${}_{-0.0043}^{+0.0041}$
Extra Gaussian noises for HIDES-S s1 (${\text{m s}}^{-1}$)	8.4${}_{-1.2}^{+1.4}$
Extra Gaussian noises for HIDES-F s2 (${\text{m s}}^{-1}$)	8.5${}_{-3.1}^{+4.7}$
Extra Gaussian noises for CES-O s3 (${\text{m s}}^{-1}$)	0 (fixed)
Extra Gaussian noises for CES-N s4 (${\text{m s}}^{-1}$)	13.1${}_{-7.1}^{+5.8}$
Extra Gaussian noises for HRS s5 (${\text{m s}}^{-1}$)	0 (fixed)
Extra Gaussian noises for AAT s6 (${\text{m s}}^{-1}$)	15.8${}_{-4.5}^{+5.6}$
Velocity offset of HIDES-F ΔRV2−1 (${\text{m s}}^{-1}$)	38.3${}_{-5.2}^{+5.2}$
Velocity offset of CES-O ΔRV3−1 (${\text{m s}}^{-1}$)	50${}_{-11}^{+11}$
Velocity offset of CES-N ΔRV4−1 (${\text{m s}}^{-1}$)	25.3${}_{-4.7}^{+4.7}$
Velocity offset of HRS ΔRV5−1 (${\text{m s}}^{-1}$)	14.5${}_{-2.5}^{+2.4}$
Velocity offset of AAT ΔRV6−1 (${\text{m s}}^{-1}$)	20.5${}_{-5.3}^{+5.6}$
Number of data of HIDES-S N1	50
Number of data of HIDES-F N2	7
Number of data of CES-O N3	5
Number of data of CES-N N4	26
Number of data of HRS N5	60
Number of data of AAT N6	13
rms (${\text{m s}}^{-1}$)	14.7

Note. Velocity offsets are the values relative to HIDES-S data.

Download table as:  ASCIITypeset image

The Hipparcos satellite made a total of 158 photometric observations for HD 47366 from 1990 March to 1993 March, and revealed photometric stability for the star down to σ_HIP = 0.007 mag, though the Hipparcos observations were not contemporaneous with the radial-velocity observations. Figure 6 shows the generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) of the Hipparcos photometric data of the star. We did not find any clue that the radial-velocity variations correlate with the brightness variations of the star. Furthermore, the rotational period of the star can be estimated to be shorter than about 86 days based on its radius $R=7.3\quad {R}_{\odot }$ and its projected rotational velocity $v\mathrm{sin}i\;=4.3\;{\text{km s}}^{-1}$ (see Table 1). This is much shorter than either the 363-day or the 685-day periods observed in radial velocity. Thus the rotational modulation of spots on the stellar surface is not considered to be a viable explanation for the observed radial-velocity variations.

Zoom In Zoom Out Reset image size

We also performed a spectral line shape analysis for the star, following the method in Sato et al. (2007). Cross-correlation profiles between pairs of stellar templates, each of which were extracted from five I₂-superposed spectra at phases of velocity $\sim 30\;{\text{m s}}^{-1}$ (phase 1), $\sim -10\;{\text{m s}}^{-1}$ (phase 2), and $\sim -50\;{\text{m s}}^{-1}$ (phase 3) by using the method of Sato et al. (2002), were derived for about 110 spectral segments (4–5 Å width each). Three bisector quantities of the cross-correlation profiles, BVS, BVC, and BVD, were then calculated, which are the velocity difference between the two flux levels of the bisector, the difference of the BVS of the upper and lower halves of the bisector, and the average of the bisector at three flux levels, respectively. The flux levels of 25%, 50%, and 75% of each cross-correlation profile were used to calculate the above three bisector quantities. As a result, we obtained BVS$\;=6.6\pm 5.5\;{\text{m s}}^{-1}$, BVC$\;=\;-2.5\pm 3.0\;{\text{m s}}^{-1}$, and BVD$\;=-86.1\pm 10.8\;{\text{m s}}^{-1}$ for the cross-correlation profile between phases 1 and 3, and BVS$\;=\;-4.1\pm 4.9\;{\text{m s}}^{-1}$, BVC$\;=3.1\pm 2.6\;{\text{m s}}^{-1}$, and BVD$\;=\;-35.1\pm 11.6\;{\text{m s}}^{-1}$ for the one between phases 2 and 3. The BVD values are consistent with the velocity differences between the phases, and the BVS and BVC values are much smaller than the BVD values. Therefore, we conclude that the observed radial-velocity variations did not originate from the distortion of the spectral lines, but from their parallel shifts, as expected for orbital motion.

The orbital parameters obtained in Section 4 were derived by fitting a double Keplerian model to the radial-velocity data, which does not guarantee that the orbits are stable over long periods of time. In order to investigate the orbital stability of the system and to further constrain orbital parameters, we performed a dynamical analysis of the system. We used a fourth-order Hermite scheme for the numerical integration (Kokubo et al. 1998).

Figure 7 shows the lifetime of the system, calculated by the 10⁶-year integrations. The lifetime here is defined as the time elapsing before the semimajor axis of one planet deviates by 10% from its initial value. All of the initial orbital parameters for the orbital integrations are fixed to the best-fit ones derived in Section 4 except for those shown in each axis of each panel. Figures 7(a)–(c) shows the lifetime of the system in the (a_c, e_c), (a_c, ω_c), and (ω_c, e_c) planes, respectively. In the panels, prograde and edge-on orbits (i_b = i_c = 90°) are assumed. The dashed line in panel (a) shows the orbit-crossing boundary where the apocenter distance of planet b is equal to the pericenter distance of planet c. The planets b and c actually have almost crossing orbits with each other in the case of the best-fit orbits (marked by a red cross in each panel), and, as shown in the figure, they are unstable. Orbits above the boundary line easily become unstable within 1000 years. Stable orbits appearing above the boundary are in the 2:1 mean-motion resonance. If a_c is near 1.93 au with an appropriate argument of pericenter, the orbits are stable with up to e_c ∼ 0.7 over 10⁶ years thanks to the resonance, though this a_c value is beyond the 3σ range of ${a}_{b}/{a}_{c}(={0.6555}_{-0.0043}^{+0.0041})$ that is well determined by our radial-velocity data (panel (b)). On the other hand, if e_c is smaller than ∼0.15 (1.4σ away from the best-fit value), the orbits are also stable even with the best-fit a_c (1.85 au) in prograde configuration (panel (c)).

Zoom In Zoom Out Reset image size

If the mutual inclination i_mut between the two planets is higher than 160°, the orbit becomes stable in a wider parameter range. Figure 7(d) shows the stability map in the plane of inclination (i_c) and the ascending node (Ω_c) of planet c (ascending node of planet b is set to Ω_b = 0). The mutual inclination of planets depends on i_b, i_c, and Ω_c − Ω_b as $\mathrm{cos}{i}_{{\rm{mut}}}=\mathrm{cos}{i}_{b}\mathrm{cos}{i}_{c}+\mathrm{sin}{i}_{b}\mathrm{sin}{i}_{c}\mathrm{cos}({{\rm{\Omega }}}_{c}-{{\rm{\Omega }}}_{b})$, and its contours are shown by black or purple lines in the figure. Since the absolute masses of the planets depend on the angle from the line of sight, the cases of three different inclinations for planet b, i_b = 30°, 60°, and 90°, are shown. The absolute mass of planet c is also changed according to i_c. The absolute masses of planets b and c are inversely proportional to $\mathrm{sin}{i}_{b}$ and $\mathrm{sin}{i}_{c}$, respectively. As seen in the figure, the orbits with i_mut ∼ 180° (i.e., mutually retrograde) are stable in all cases. This is because planet encounters occur in a shorter time with higher mutual velocities in the retrograde system compared to the prograde system. This is similar to the case of the two brown-dwarf candidates orbiting BD+20 2457, in which the system is unstable for the prograde orbits, while it is stable for the retrograde orbits (Horner et al. 2014).

Figure 8 shows the 7 × 10⁴ year evolution of the system. The upper panel shows the evolution of the coplanar retrograde system with the best-fit orbital parameters derived in Section 4. Although the outer planet violates the orbit of the inner planet, they are stable over 10⁶ years. The bottom panel shows a prograde case with the same orbital parameters except assuming e_c,ini = 0. As seen in the figure, if the eccentricity of the outer planet is small, the system can keep stable orbits, while the system with the outer planet in eccentric orbit quickly causes orbital crossing and becomes unstable (see Figure 7).

Zoom In Zoom Out Reset image size

Figure 9 shows maps of the stability index $D=| \langle {n}_{2}\rangle -\langle {n}_{1}\rangle |$ for the system. From a 5000-year numerical integration, we obtain an average of mean motion $\langle {n}_{1}\rangle$ of planet c and subtract it from an average of mean motion $\langle {n}_{2}\rangle$ obtained in the next 5000-year orbit (see Couetdic et al. 2010 for the details of the stability analysis). The left and right panels are the maps for prograde and retrograde configuration, respectively. From a comparison with panel (a) of Figure 7, we find that the system is regular when ${\mathrm{log}}_{10}D\leqslant -3$. In the case of the retrograde orbit, the orbit is stable even if the eccentricity of the planet c exceeds the orbit-crossing boundary (dashed line) by ∼0.1.

Zoom In Zoom Out Reset image size

We have reported the detection of a double giant-planet system around the K1 giant HD 47366 (M_⋆ = 1.81 ± 0.13 M_⊙, R_⋆ = 7.30 ± 0.33 R_⊙) from precise radial-velocity measurements at OAO, Xinglong, and AAT. The inner (planet b) and outer planets (planet c) have minimum masses of ${1.75}_{-0.17}^{+0.20}\quad {M}_{{\rm{J}}}$ and ${1.86}_{-0.15}^{+0.16}\quad {M}_{{\rm{J}}}$, semimajor axes of ${1.214}_{-0.029}^{+0.030}$ au and ${1.853}_{-0.045}^{+0.045}$ au, and eccentricities of ${0.089}_{-0.060}^{+0.079}$ and ${0.278}_{-0.094}^{+0.067}$, respectively, which were derived by fitting a double Keplerian model to the radial-velocity data. The period ratio of the two planets is slightly smaller than 2, and thus the system adds to the growing population of such multi-giant-planet systems with small orbital separation around evolved intermediate-mass stars.

The dynamical analysis of the system revealed, however, that the best-fit orbits in prograde configuration are unstable. We found that it is stable in the following cases: (1) the two planets are in the 2:1 mean-motion resonance, (2) the eccentricity of planet c is less than ∼0.15, (3) the mutual inclination the two planets is larger than 160° (i.e., retrograde configuration).

If the two planets are in the stable 2:1 resonance, the ratio of the semimajor axis should be a_b/a_c ∼ 0.63 (a_c ∼ 1.93 au) and the eccentricity of planet c should be 0.5 ≲ e_c ≲ 0.7 or 0 ≤ e_c ≲ 0.2 (see Section 5 and the darkest-blue regions in panel (a) of Figure 7). On the other hand, a_b/a_c is well determined from our observations to be ${a}_{b}/{a}_{c}={0.6555}_{-0.0043}^{+0.0041}$ and e_c is determined to be ${e}_{c}={0.278}_{-0.094}^{+0.067}$. Both of these values deviate from those for the high-eccentricity case with more than 3σ, and thus from the view point of the observational data, it is less likely that the planets are in the high-eccentricity resonance. We cannot completely reject the possibility of the low-eccentricity case. However, considering the semimajor-axis ratio, it is also less likely that the planets are in the resonance.

It is more plausible that the eccentricity of planet c is less than ∼0.15. The value is just 1.4σ from the best-fit value, and the best-fit orbit is consistent with a circular one to the 3σ level. Thus, it is reasonable to think that the two planets actually have nearly circular orbits and thus they are stable. In this case, the period ratio of the planets is slightly smaller than 2, and they are not in 2:1 resonance.

The third possibility is the most extreme one. Basically, we cannot determine the inclination or the ascending node of a planetary orbit by radial-velocity measurements, and thus we cannot discriminate between prograde and retrograde orbits observationally. Gayon & Bois (2008) and Gayon-Markt & Bois (2009) studied mutually retrograde configurations for some multi-planet systems around solar-type stars and showed that the configurations are dynamically stable and consistent with the radial-velocity observations. Recently, Horner et al. (2014) studied the dynamical stability of a double brown-dwarf system around the giant BD+20 2457 (Niedzielski et al. 2009), and showed that the prograde orbits best-fitted to the radial-velocity data are unstable, while those in retrograde configuration are stable. The HD 47366 system is similar to the case of the BD+20 2457 system. Our dynamical analysis showed that the coplanar retrograde configuration for the HD 47366 system is stable in a wider range of orbital parameters including the eccentric ones which are unstable for prograde configuration (see Figure 8). Confirmation of the high eccentricity of the outer planet by continuous observations will strongly support the retrograde hypothesis for the HD 47366 system, though the formation mechanism of such a mutually retrograde system is largely unknown. In any of the above cases, the planetary orbits can be altered by the close encounter of the planets in a relatively short period of time, and thus more detailed dynamical modeling for the observed radial-velocity variations is necessary in order to provide a definitive orbital solution.

The stellar and planetary parameters of HD 47366 system are similar to those of the η Cet system; two massive planets (${m}_{b}\mathrm{sin}i=2.6\quad {M}_{{\rm{J}}}$, ${m}_{c}\mathrm{sin}i=3.3\quad {M}_{{\rm{J}}}$), with periods of P_b = 407 days and P_c = 740 days and eccentricities of e_b = 0.12 and e_c = 0.08 orbit around a K giant star with a mass of 1.7 M_⊙ and a radius of 14.3 R_⊙ (Trifonov et al. 2014). Trifonov et al. (2014) revealed that, for a coplanar configuration, the system can be stable by being trapped in an anti-aligned 2:1 mean-motion resonance in a region with moderate e_b eccentricity, which lies about 1σ away from the best-fit Keplerian orbits. A much larger non-resonant stable region also exists in low-eccentricity parameter space, though it appears to be much farther from the best fit than the 2:1 resonant region. This is in contrast with the case of HD 47366 in which low eccentricities of the planets are more plausible than 2:1 resonant configuration. Another possible 2:1 resonant system around an evolved star is 24 Sex, in which two massive planets (${m}_{b}\mathrm{sin}i=1.99\quad {M}_{{\rm{J}}}$, ${m}_{c}\mathrm{sin}i=0.86\quad {M}_{{\rm{J}}}$), with periods of P_b = 452.8 days and P_c = 883.0 days and eccentricities of e_b = 0.09 and e_c = 0.29 orbit around a G-type subgiant with a mass of 1.54 M_⊙ and a radius of 4.9 R_⊙ (Johnson et al. 2011). The period ratio of the planets is the closest to two of the three systems. Wittenmyer et al. (2012a) revealed that the best-fit orbit of the system is mutually crossing and it is only dynamically feasible if the planets are in 2:1 resonance, which can protect them from close encounters. As described, all of the above three systems are near the 2:1 resonance, but they show different properties in their orbital configuration. Further analysis and comparison of the systems would help us understand the dynamics of planetary systems in more detail.

It is unknown why multi-giant-planet systems with small orbital separations are preferentially found around evolved intermediate-mass stars (see Figure 1). It may be a primordial property of planets around intermediate-mass stars that could be an outcome of planet formation or an acquired one as a result of orbital evolution caused by stellar evolution (stellar tide and mass loss) of central stars. In the case of HD 47366, if the star is a helium-core burning star that has passed through the RGB tip (see Section 2), the planetary orbits could have been migrated outward by ∼10% (∼0.1 au) during the RGB phase because of the effect of mass loss of the central star (Kunitomo et al. 2011; Villaver et al. 2014). However, the effect works for both of the two planets and thus it could not result in making the orbital separation smaller. The star might have harbored a third planet that was engulfed and possibly partially responsible for the current orbital configuration, though the star does not exhibit a possible signature of the planet engulfment such as an overabundance of lithium (e.g., Siess & Livio 1999; Adamow et al. 2012), Liu et al. (2014) obtained the lithium abundance of¹⁷ A(Li) = 0.38 for HD 47366, suggesting that the star is an Li-depleted giant. Considering also that many of the systems with small orbital separation are found around less evolved subgiants, the effect of stellar evolution may be less significant compared to that of stellar mass. Unfortunately, the ongoing Doppler planet searches cannot clearly discriminate between the above two factors, stellar mass and stellar evolution, because the targets of the searches are both evolved and intermediate-mass stars. Investigating planets around intermediate-mass, main-sequence stars or those around evolved low-mass (∼1 M_⊙) stars would help clarify each effect separately.

This research is based on data collected at the Okayama Astrophysical Observatory (OAO), which is operated by National Astronomical Observatory of Japan, Xinglong Station, which is operated by the National Astronomical Observatory of China, and at the Australian Astronomical Observatory. We are grateful to all the staff members of the observatories for their support during the observations. This work was partially Supported by the Open Project Program of the Key Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese Academy of Sciences. We thank students of Tokyo Institute of Technology and Kobe University for their kind help with the observations at OAO. B.S. was partially supported by MEXT's program "Promotion of Environmental Improvement for Independence of Young Researchers" under the Special Coordination Funds for Promoting Science and Technology, and by Grant-in-Aid for Young Scientists (B) 17740106 and 20740101 and Grant-in-Aid for Scientific Research (C) 23540263 from the Japan Society for the Promotion of Science (JSPS). M.N. is supported by Grant-in-Aid for Young Scientists (B) 21740324 and H.I. is supported by Grant-In-Aid for Scientific Research (A) 23244038 from JSPS. This work was partially funded by the National Natural Science Foundation of China under grants 111173031, 1233004, and 11390371, as well as the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences, grant No. XDB09000000, and by the JSPS under Grant-in-Aid for Scientific Research (B) 17340056 (H.A.) and grant 08032011-000184 in the framework of the Joint Research Project between China and Japan. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.