THE MASS OF HD 38529c FROM HUBBLE SPACE TELESCOPE ASTROMETRY AND HIGH-PRECISION RADIAL VELOCITIES^*

We first model RV data, a significant fraction of which comes from the Hobby–Eberly Telescope (HET). Measurements from the California–Carnegie exoplanet research group (Wright et al. 2009) and a few from the McDonald Harlan J. Smith telescope (Wittenmyer et al. 2009) were also included. The California–Carnegie data were particularly valuable, increasing the time span from 4 to over 11 yr. Our astrometry covers only ∼70% of the orbit of HD 38529c, and in the absence of a multi-period span of RVs, would not be sufficient to establish accurate perturbation elements, particularly period, eccentricity, and periastron passage. All RV sources are listed in Table 2, along with the rms of the residuals to the combined orbital fits described below. The errors for all published RV and our new HET RV have been modified by adding in quadrature the expected RV jitter from stellar activity determined in Section 2.2.

All the RVs were obtained using I₂ cell techniques. The HET data were obtained with the HET High-resolution Spectrograph, described in Tull (1998) and processed with the I₂ pipeline described in Bean et al. (2007), utilizing robust estimation to combine the all velocities from the individual chunks. Typically three HET observations are secured within 10–15 minutes. These are combined using robust estimation to form normal points for each night. The HET normal points and associated errors are listed in Table 3.

Combining RV observations from different sources is possible in the modeling environment we use. GaussFit (Jefferys et al. 1988) has the capability to simultaneously solve for many separate velocity offsets (because velocities from different sources are relative, having differing zero points), along with the other orbital parameters. Relative offsets (γ) and associated errors are listed in Table 11.

Orbital parameters derived from a combination of HET, HJS, Lick, and Keck RVs and HST astrometry will be provided in Section 3.5. Figure 1 shows the entire span of data along with the best-fit multiple-Keplerian orbit. We note that there is sufficient bowing in the residuals to justify continued low-cadence RV monitoring, particularly given the prediction of Moro-Martín et al. (2007) of dynamical stability for planets with periods as short as ∼70 yr. Subtracting in turn the signature of first one, then the other known companion from the original velocity data we obtain the component b and c RV orbits shown in Figure 2, each phased to the relevant periods. Re-iterating, all RV fits were modeled simultaneously with the astrometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared to the typical perturbation RV curve (e.g., Hatzes et al. 2005; McArthur et al. 2004; Cochran et al. 2004), our original orbits for components b and c exhibited significant scatter, much due to the identified stellar noise source discussed in Section 2.2 below. Periodogram analysis of RV residuals to simultaneous fits of components b and c indicated a significant peak with a period near 197 days. The existence of the signal is fairly secure. A bootstrap analysis carried out by randomly shuffling the RV residual values 200,000 times (keeping the times fixed), and determining if the random data periodogram had peaks higher than the real data periodogram in the frequency range 0 < ν < 0.02 day⁻¹, yielded a false alarm probability, FAP = 5 × 10⁻⁴. This motivated the addition of a third Keplerian component, resulting in the fit shown in the bottom panel of Figure 2. Even though the amplitude of this signal is about that expected from stellar noise, including this component (five additional parameters) in the combined modeling improved both the reduced χ², and the rms scatter as shown in Table 2. Identification of the cause of the signal will be discussed further in Section 5.

There are a number of sources of RV noise intrinsic to HD 38529: pulsations and velocity perturbations introduced by star spots and/or plages. The velocity effects caused by the latter two are modulated by stellar rotation. Valenti & Fischer (2005) measure a rotation of HD 38529, V_rot sin  i = 3.5 ± 0.5 km s⁻¹. HD 38529 is subgiant star, evolving toward the giant branch of the Hertszprung–Russell diagram, and is expected to have a higher level of pulsational activity than a main sequence star (Hatzes & Zechmeister 2008). The pulsational amplitude can be estimated using the scaling relationship of Kjeldsen & Bedding (1995)V_amp = (L/L_☉)/(M/M_☉) × 0.234 m s⁻¹. The luminosity and mass of HD 38529 yield a pulsational amplitude of ∼1 m s⁻¹, so this alone cannot account for the excess RV scatter.

Several relationships between the amplitude of RV noise and the fraction of star spot coverage have been developed. Saar & Donahue (1997) obtain A_RV = 6.5 × v sin  i × f^0.9, where f is the spot filling factor in percent. Hatzes (2002) obtained A_RV = (8.6v sin  i − 1.6) × f^0.9. We can estimate the spot filling factor from FGS photometry of HD 38529. The FGS has been shown to be a photometer precise at the 2 mmag level (Benedict et al. 1998). We flat-fielded the HD 38529 FGS photometry, using an average of the counts from the astrometric reference stars listed in Table 5, and plotted it against time. Clearly not constant at a level ten times our internal precision, a Lomb–Scargle periodogram showed a significant period at P = 31.6 days (FAP = 4.3 × 10⁻⁴). A sine wave fit to the photometry yielded P = 31.65 ± 017 with an amplitude = 1.5 ± 0.2 mmag. Figure 3 is a plot of these photometric data phased to that period.

Zoom In Zoom Out Reset image size

The Valenti & Fischer (2005) V_rot sin  i = 3.5 km s⁻¹ and a stellar radius from Baines et al. (2008b), R = 2.44 ± 0.22 R_☉, would predict a minimum P_rot = 32 ± 5  days. Interpreting the modulation period of 31.6 days as the stellar rotation period, we ascribe the photometric variation (0.15%) to rotational modulation of star spots. The photometric amplitude suggests an RV noise level of 4–5 m s⁻¹. Taking the HET velocity rms as closer to the true RV variation, we identify the remaining RV scatter as a combination of the three effects identified.

The astrometric reference frame for HD 38529 consists of four stars. Any prior knowledge concerning these four stars eventually enters our modeling as observations with error and yields the most accurate parallax and proper motion for the prime target, HD 38529. These periodic and non-periodic motions must be removed as accurately and precisely as possible to obtain the perturbation inclination and size caused by HD 38529c. Of particular value are independently measured proper motions. This particular prior knowledge comes from the UCAC3 catalog (Zacharias et al. 2010). Figure 5 shows the distribution in FGS1r pickle coordinates of the 23 sets of four reference star measurements for the HD 38529 field. At each epoch we measured each reference stars two to four times, and HD 38529 five times.

Because the parallax determined for HD 38529 is measured with respect to reference frame stars which have their own parallaxes, we must either apply a statistically derived correction from relative to absolute parallax (van Altena et al. 1995, Yale Parallax Catalog, YPC95), adopt an independently derived parallax (e.g., Hipparcos), or estimate the absolute parallaxes of the reference frame stars. In principle, the colors, spectral type, and luminosity class of a star can be used to estimate the absolute magnitude, M_V, and V-band absorption, A_V. The absolute parallax for each reference star is then simply,

$$\begin{equation} \pi _{\rm abs} = 10^{-(V-M_V+5-A_V)/5}. \end{equation} \tag{ 1 }$$

3.2.1. Reference Star Photometry and Spectroscopy

Our band passes for reference star photometry include BVRI photometry of the reference stars from the NMSU 1 m telescope located at Apache Point Observatory and JHK (from the Two Micron All Sky Survey (2MASS) ¹⁰). Table 7 lists the visible and infrared photometry for the HD 38529 reference stars. The spectra from which we estimated spectral-type and luminosity class were obtained on 2009 December 9 using the RCSPEC on the Blanco 4 m telescope at CTIO. We used the KPGL1 grating to give a dispersion of 0.95 Å pixel⁻¹. Classifications used a combination of template matching and line ratios. The spectral types for the higher S/N stars are within ±1 subclass. Classifications for the lower S/N stars are ±2 subclasses. Table 7 lists the spectral types and luminosity classes for our reference stars.

Table 7. Astrometric Reference Star Adopted Spectrophotometric Parallaxes

ID	Sp. T.a	V	MV	m − M	AV	πabs(mas)
2	F2V	14.12	3	11.12	1.302	1.1 ± 0.3
4	F0V	13.05	2.7	10.35	1.147	1.4 0.3
5	K0III	13.87	0.7	13.17	1.395	0.4 0.1
6	K1III	14.34	0.6	13.74	1.271	0.3 0.1

Note. ^aSpectral types and luminosity class estimated from classification spectra, colors, and reduced proper motion diagram (Figures 6 and 7).

Download table as:  ASCIITypeset image

Figure 6 contains a (J − K) versus (V − K) color–color diagram of the reference stars. Schlegel et al. (1998) find an upper limit A_V ∼ 2 toward HD 38529, consistent with the absorptions we infer comparing spectra and photometry (Table 7).

Zoom In Zoom Out Reset image size

The derived absolute magnitudes are critically dependent on the assumed stellar luminosity, a parameter impossible to obtain for all but the latest type stars using only Figure 6. To confirm the luminosity classes obtained from classification spectra we abstract UCAC3 proper motions Zacharias et al. (2010) for a one-degree-square field centered on HD 38529, and then iteratively employ the technique of reduced proper motion (Yong & Lambert 2003; Gould & Morgan 2003) to discriminate between giants and dwarfs. The end result of this process is contained in Figure 7. Reference stars ref-5 and ref-6 are confirmed to have luminosity class III (giant).

Zoom In Zoom Out Reset image size

3.2.2. Adopted Reference Frame Absolute Parallaxes

The HD 38529 reference frame contains only four stars. From positional measurements we determine the scale, rotation, and offset "plate constants" relative to an arbitrarily adopted constraint epoch for each observation set. As for all our previous astrometric analyses, we employ GaussFit (Jefferys et al. 1988) to minimize χ². The solved equations of condition for the HD 38529 field are

$$\begin{equation} x^\prime = x + lc_x(\it B-V) - \Delta XF_x \end{equation} \tag{ 2 }$$

$$\begin{equation} y^\prime = y + lc_y(\it B-V) - \Delta XF_y \end{equation} \tag{ 3 }$$

$$\begin{equation} \xi = Ax^\prime + By^{\prime } + C - \mu _\alpha \Delta t - P_\alpha \pi - \mbox{ORBIT}_x \end{equation} \tag{ 4 }$$

$$\begin{equation} \eta = Dx^\prime + Ey^\prime + F - \mu _\delta \Delta t - P_\delta \pi - \mbox{ORBIT}_y \end{equation} \tag{ 5 }$$

for FGS1r data. Identifying terms, x and y are the measured coordinates from HST; (B − V) is the Johnson (B − V) color of each star; and lc_x and lc_y are the lateral color corrections, applied only to FGS1r data. Here XF_x and XF_y are cross filter corrections (see Benedict et al. 2002b) in x and y applied to the observations of HD 38529. A, B, D and E are scale and rotation plate constants and C and F are offsets; μ_α and μ_δ are proper motions; Δt is the epoch difference from the constraint epoch; P_α and P_δ are parallax factors; and π is the parallax. We obtain the parallax factors from a JPL Earth orbit predictor (Standish 1990), upgraded to version DE405. Orientation to the sky for the FGS1r data is obtained from ground-based astrometry (2MASS catalog) with uncertainties of 001. ORBIT_x and ORBIT_y are functions (through Thiele–Innes constants, e.g., Heintz 1978) of the traditional astrometric and RV orbital elements listed in Table 11.

The Optical Field Angle Distortion (OFAD) calibration (McArthur et al. 2002) reduces as-built HST telescope and FGS1r distortions with magnitude ∼1'' to below 2 mas over much of the FGS1r field of regard. These data were calibrated with a revised OFAD generated by McArthur in 2007. From histograms of the FGS astrometric residuals (Figure 8) we conclude that we have obtained correction at the ∼1 mas level. The reference frame "catalogs" from FGS1r in ξ and η standard coordinates (Table 8) were determined with 〈σ_ξ〉 = 0.15 and 〈σ_η〉 = 0.15 mas.

Zoom In Zoom Out Reset image size

Table 9. Final Reference Star Proper Motions

ID	V	μαa	$\sigma _{\mu _\alpha }$	μδa	$\sigma _{\mu _\delta }$
2	14.12	−13.32	0.15	15.45	0.12
4	13.05	2.32	0.09	7.87	0.09
5	13.87	−12.26	0.11	13.13	0.09
6	14.34	3.76	0.10	−4.91	0.09

Note. ^aμ_α and μ_δ are relative motions in mas yr⁻¹.

Download table as:  ASCIITypeset image

We linearly combine unperturbed Keplerian orbits, simultaneously modeling the RVs and astrometry. The period (P), the epoch of passage through periastron in years (T), the eccentricity (e), and the angle in the plane of the true orbit between the line of nodes and the major axis (ω), are the same for an orbit determined from RV or from astrometry. The remaining orbital elements come only from astrometry. We force a relationship between the astrometry and the RV through this constraint (Pourbaix & Jorissen 2000)

$$\begin{equation} \displaystyle {{\alpha \sin i \over \pi _{\rm abs}} = {P K (1 - e^2)^{1/2}\over 2\pi \times 4.7405}}, \end{equation} \tag{ 6 }$$

where quantities derived only from astrometry (parallax, π_abs, primary perturbation orbit size, α, and inclination, i) are on the left, and quantities derivable from both (the period, P and eccentricity, e), or RVs only (the RV amplitude of the primary, K), are on the right. In this case, given the fractional orbit coverage of the HD 38529c perturbation afforded by the astrometry, all right-hand side quantities are dominated by the RVs.

For the parameters critical in determining the mass of HD 38529 we find a parallax, π_abs = 25.11 ± 0.19 mas and a proper motion in R.A. of −78.60 ± 0.15 mas yr⁻¹ and in decl. of −141.96 ± 0.11 mas yr⁻¹. Table 10 compares values for the parallax and proper motion of HD 38529 from HST and both the original Hipparcos values and the recent Hipparcos re-reduction (van Leeuwen 2007). We note satisfactory agreement for parallax and proper motion. Our precision and extended study duration have significantly improved the accuracy and precision of the parallax and proper motion of HD 38529.

We find a perturbation size, α_c = 1.05 ± 0.06 mas, and an inclination, i = 483 ± 37. These, and the other orbital elements are listed in Table 11 with 1σ errors. Figure 9 illustrates the Pourbaix & Jorrissen relation (Equation 6) between parameters obtained from astrometry (left side) and RVs (right side) and our final estimates for α_c and i. In essence, our simultaneous solution uses the Figure 9 curve as a quasi-Bayesian prior, sliding along it until the astrometric and RV residuals are minimized. Gross deviations from the curve are minimized by the high precision of all of the right-hand side terms in Equation (6) (Tables 10 and 11).

Zoom In Zoom Out Reset image size

At this stage we can assess the reality of the HD 38529c astrometric perturbation by plotting the astrometric residuals from a model that does not include a component c orbit. Figure 10 shows the R.A. and decl. components of the FGS residuals plotted as small symbols. We also plot normal points formed from those smallest symbols within each of the 23 data sets listed in Table 4. The largest symbols denote the final normal points formed for each of our (effectively) five epochs. Each plot contains as a dashed line the R.A. and decl. components of the perturbation we find by including an orbit in our modeling. Finally, Figure 11 shows the perturbation on the sky with our normal points.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The planetary mass depends on the mass of the primary star, for which we have adopted M_* = 1.48 M_☉ (Takeda et al. 2007). For that M_* we find M_c = 17.6^+1.5_−1.2 M_Jup, a significant improvement over the Reffert & Quirrenbach (2006) estimate, M_c = 38⁺³⁶₋₁₉ M_Jup, but agreeing within the errors. HD 38529c is likely a brown dwarf, but only about 3σ from the "traditional" planet–brown dwarf dividing line, 13 M_Jup, the mass above which deuterium is thought to burn. In Table 11 the mass value, M_c, incorporates the present uncertainty in M_*. However, the dominant source of error is in the inclination estimate. Until HD 38529c is directly detected, its radius is unknown. Comparing to the one known transiting brown dwarf, CoRot-Exo-3b (Deleuil et al. 2008), a radius of R ∼ 1R_JUP seems reasonable.

Given the adopted Table 11 errors, HD 38529c is either one of the most massive exoplanets or one of the least massive brown dwarfs. We can compare our true mass to the (as of 2010 January) 69 transiting exoplanetary systems, each also characterized by true mass, not Msin i. As shown in the useful Exoplanetary Encyclopedia (Schneider 2009) only one companion (CoRoT-Exo-3 b, Deleuil et al. 2008) has a mass in excess of 13 M_Jup. Whether this "brown dwarf desert" (Grether & Lineweaver 2006) in the transiting sample is due to the difficulty of migrating high-mass companions (bringing them in close enough to increase the probability of transit), or to inefficiencies in gravitational instability formation is unknown.

The age of the host star, ∼3.3 By, would suggest that HD 38529c has not yet cooled to an equilibrium temperature. An estimated temperature and self-luminosity for a 17 M_Jup object that is 3.3 By old can be found from the models of Hubbard et al. (2002). Those models predict that HD 38529c has an effective temperature, T_eff ≃ 400 K, and L = 2.5e − 7 L_☉, about 20 times brighter than what we estimated using these same models for Eri b (Benedict et al. 2006). Unfortunately HD 38529 has about 16 times the intrinsic brightness of Eri, erasing any gain in contrast. We note that due to the eccentricity of the orbit, HD 38529c is actually within the present-day habitable zone for a fraction of its orbit. As HD 38529 continues to evolve and brighten, the habitable zone will move outward and HD 38529c will be in that zone for some period of time.

If the inner known companion HD 38529b is a minimum mass exoplanetary object (assuming M = Msin i = 0.8 M_Jup), our 1 mas astrometric per-observation precision precludes detecting that 2 mas of arc signal. Invoking (with no good reason) coplanarity with HD 38529c similarly leaves us unable to detect HD 38529b. However, with the motivation of our previous result for HD 33636 (Bean et al. 2007), we can test whether or not HD 38529b is also stellar by establishing an upper limit from our astrometry. To produce a perturbation, α_b > 0.2 mas (a 3σ detection, given σ_α = 0.06 mas from Table 11), and the observed RV amplitude, K_b = 59 m s⁻¹, requires M_b ∼ 0.1 M_☉ in an orbit inclined by less than 05. Our limit is lower than that established with the CHARA interferometer (Baines et al. 2008a), who established a photometric upper limit of G5 V for the b component. While it might be possible to use 2MASS and SDSS (Ofek 2008) photometry of this object to either confirm or eliminate a low-mass stellar companion by backing out a possible contribution from an M, L, or T dwarf, using their known photometric signatures (Hawley et al. 2002; Covey et al. 2007), we lack precise (1%) knowledge of the intrinsic photometric properties of a sub-giant star in the Hertzsprung gap with which to compare.

In summary, RVs from four sources, Lick and Keck (Wright et al. 2009), HJS/McDonald (Wittenmyer et al. 2009), and our new high-cadence series from the HET, were combined with HST astrometry to provide improved orbital parameters for HD 38529b and HD 38529c. Rotational modulation of star spots with a period P = 31.66 ± 0.17 days produces 0.15% photometric variations, spot coverage sufficient to produce the observed residual RV variations. Our simultaneous modeling of RVs and over three years of HST FGS astrometry yields the signature of a perturbation due to the outermost known companion, HD 38529c. Applying the Pourbaix & Jorrissen constraint between astrometry and RVs, we obtain for the perturbing object HD 38529c a period, P = 2136.1 ± 0.3 days, inclination, i = 483 ± 37, and perturbation semimajor axis, α_c = 1.05 ± 0.06 mas. Assuming for HD 38529a stellar mass M_* = 1.48 ± 0.05 M_☉, we obtain a mass for HD 38529c, M_b = 17.6^+1.5_−1.2 M_Jup, within the brown dwarf domain. Our independently determined parallax agrees within the errors with Hipparcos, and we find a close match in proper motion. Our HET RVs combined with others establish an upper limit of about one Saturn mass for possible companions in a dynamically stable range of companion-star separations, 0.2 ⩽ a ⩽ 1.2 AU. RV residuals to a model incorporating components b and c contain a signal with an amplitude equal to the rms variation with a period, P ∼194  days and an inferred a ∼ 0.75 AU. While dynamical simulations do not rule out interpretation as a planetary mass companion, the low S/N of the signal argues for confirmation.