MagAO IMAGING OF LONG-PERIOD OBJECTS (MILO). I. A BENCHMARK M DWARF COMPANION EXCITING A MASSIVE PLANET AROUND THE SUN-LIKE STAR HD 7449^*

We observed HD 7449 using the Magellan Clay Telescope at the Las Campanas Observatory in Chile on the nights of UT 2014 November 5 and 22. We used MagAO paired with VisAO and Clio-2, for which we used the narrow camera (plate scale = 001585; Morzinski et al. 2015). On the first night, the observing conditions were fair, with seeing varying around 1${}^{\prime\prime }$, therefore only 200 modes of AO correction were employed. We observed the star with VisAO at Ys (0.99 μm) and with Clio-2 at H (1.65 μm) and Ks (2.15 μm). Unsaturated photometric images were also acquired in each filter. On the second night, the seeing was much better, with stable seeing under 1${}^{\prime\prime }$, therefore the maximum 300 modes of AO correction were employed. We observed the star with VisAO at ${r}^{\prime }$ (0.63 μm), ${i}^{\prime }$ (0.77 μm), ${z}^{\prime }$ (0.91 μm), and with Clio-2 at J (1.1 μm). Unsaturated photometric images were acquired in each filter. All observations were acquired with the instrument rotator off to enable angular differential imaging (ADI, Marois et al. 2006).

A bright object was identifiable in the raw images at each wavelength, separated by ∼05 from the star. Therefore ADI PSF subtraction was not needed to enhance contrast, and little integration was required in each filter. We obtained total integrations of 2.3 minutes at ${r}^{\prime }$, 1.2 minutes at ${i}^{\prime }$, 1.17 minutes at ${z}^{\prime }$, 1.9 minutes at Ys, 0.5 minutes at H, 18.67 minutes at J, and 4.33 minutes at Ks.

All data reduction was performed with custom scripts in Matlab. The Clio-2 images were divided by the number of coadds, corrected for nonlinearity (Morzinski et al. 2015), divided by the integration times, sky-subtracted, and then registered and cropped. The VisAO images were dark-subtracted, divided by the integration times, and then registered and cropped. All images were rotated to north-up, east-left and then median-combined into final images at each wavelength. Finally, 2D radial profiles were subtracted from each image to remove the majority of the stellar flux (see Figure 1). The object at ∼054 is unambiguously detected in each filter.

Zoom In Zoom Out Reset image size

RV data on HD 7449 were first acquired as part of the Magellan Planet Search Program, which originally made use of the MIKE echelle spectrometer (Bernstein et al. 2003) on the Magellan Clay telescope until 2009 September. The reported precision achieved by MIKE was 5 m s⁻¹ on solar-type stars (Minniti et al. 2009). Observations with MIKE were made using a 035 slit, which results in a spectral resolution of R ∼ 70,000 in the blue and ∼ 50,000 in the red. The wavelength coverage ranges from 3900 to 6200 Å, capturing the iodine region (5000–6300 Å), and is divided into two CCDs covering the red and blue wavelength regions.

HD 7449 was subsequently observed using the Carnegie PFS (Crane et al. 2010), a temperature-controlled high resolution spectrograph, which now carries out all observations for the Magellan Planet Search Program. PFS covers 3880–6680 Å and the 05 slit is used, which results in a spectral resolution of ∼80,000 in the iodine region. Continuous monitoring of stable stars reveals that the Magellan/PFS system achieves an average measurement precision of 1.5 m s⁻¹ (Arriagada et al. 2013).

The RVs for both instruments were obtained using the iodine technique (Butler et al. 1996). Briefly, an iodine absorption cell provides the wavelength scale and instrumental PSF for each stellar observation, which are computed in 2 Å chunks. A forward modeling procedure of each observation is carried out for each chunk, thus providing an individual measurement of the wavelength, PSF, and Doppler shift. The final measured RV is the weighted average of all the chunks for a given observation. Internal uncertainties are computed as the standard deviation of the velocities derived from each chunk. The new RVs for HD 7449 obtained from MIKE and PFS are listed in Table 1.

Table 1.  RVs for HD 7449

Julian Date	RV (m s−1)	${\sigma }_{\mathrm{RV}}$ (m s−1)	Instrument
2451459.55882	78.57	10.00	1
2451480.14706	86.90	10.00	1
2451490.44118	76.19	10.00	1
2451541.91176	76.90	10.00	1
2451747.79412	52.62	10.00	1

Note. Instrument 1 corresponds to CORALIE (Dumusque et al. 2011), 2 corresponds to HARPS, 3 corresponds to Magellan/MIKE, and 4 corresponds to Magellan/PFS.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We also included in our analysis RVs measured with HARPS and CORALIE. These RVs were originally reported in Dumusque et al. (2011). However, HD 7449 has been observed by HARPS since that publication, so we downloaded all available HARPS data on HD 7449 from the ESO archive. Starting from the ESO extracted and calibrated spectra, we obtained new Doppler measurements using the HARPS-TERRA software (Anglada-Escudé & Butler 2012). The CORALIE data were not explicitly reported by Dumusque et al. (2011), nor are they available in any archive, so we used DataThief (http://datathief.org) to retrieve the RVs. To account for possible errors in the extraction, we assumed 10 m s⁻¹ errors for the CORALIE data in our subsequent RV analysis. The entire RV data set is shown in Figure 2(a), revealing the clear long-term, parabolic trend, and the individual RVs are reported in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Photometry was measured as follows. First, a circular aperture of radius = 1 FWHM, corresponding to the size of a diffraction-limited PSF at each wavelength, was placed at the detected object's photocenter in each image. The same aperture was placed at the stellar photocenter in each unsaturated image, and then the fluxes within all the apertures were summed. Uncertainties were calculated as the standard deviations of the fluxes in apertures placed around the star at the same radius.

Astrometry was measured by calculating the photocenters in the same apertures, and astrometric uncertainties were assumed to be 5 mas at each wavelength based on previous imaging with MagAO (e.g., Rodigas et al. 2015). Table 2 lists the object's photometry and astrometry. The object has a separation of ∼054 and P.A. ∼340°. Because the star has high proper motion (van Leeuwen 2007), the two epochs of direct detections separated by only 17 days are enough to show that the object is inconsistent with being background at 2σ confidence. In Section 3.2, we will show that the object's SED confirms that it is unlikely to be background.

Table 2.  HD 7449B Photometry and Astrometry

Parameter	Value
${\rm{\Delta }}{r}^{\prime }$ (0.63 μm)	${8.82}_{-0.11}^{+0.13}$
${\rm{\Delta }}{i}^{\prime }$ (0.77 μm)	${7.32}_{-0.11}^{+0.13}$
${\rm{\Delta }}{z}^{\prime }$ (0.91 μm)	${6.53}_{-0.13}^{+0.15}$
${\rm{\Delta }}{Ys}$ (0.99 μm)	${5.87}_{-0.23}^{+0.29}$
${\rm{\Delta }}{J}_{\mathrm{MKO}}$ (1.1 μm)	${5.81}_{-0.10}^{+0.11}$
${\rm{\Delta }}{H}_{\mathrm{MKO}}$ (1.65 μm)	${5.11}_{-0.10}^{+0.11}$
${\rm{\Delta }}{{Ks}}_{\mathrm{Barr}}$ (2.15 μm)	${4.85}_{-0.03}^{+0.03}$
${M}_{{r}^{\prime }}$	13.39 ± 0.17
${M}_{{i}^{\prime }}$	11.51 ± 0.17
${M}_{{z}^{\prime }}$	10.75 ± 0.20
${M}_{J}$	9.26 ± 0.16
${M}_{H}$	8.33 ± 0.16
${M}_{{K}_{s}}$	7.97 ± 0.09
${\rm{\Delta }}{{\rm{R.A.}}}_{{t}_{1}}$ (${}^{\prime\prime }$)	−0.19 ± 0.003
${\rm{\Delta }}\mathrm{Decl}{.}_{{t}_{1}}$ (${}^{\prime\prime }$)	0.52 ± 0.003
${\rm{\Delta }}{{\rm{R.A.}}}_{{t}_{2}}$ (${}^{\prime\prime }$)	−0.19 ± 0.006
${\rm{\Delta }}\mathrm{Decl}{.}_{{t}_{2}}$ (${}^{\prime\prime }$)	0.51 ± 0.005
${\rho }_{{t}_{1}}$ (${}^{\prime\prime }$)	0.55 ± 0.007
${{\rm{P.A.}}}_{{t}_{1}}$ (°)	339.99 ± 1.84
${\rho }_{{t}_{2}}$ (${}^{\prime\prime }$)	0.54 ± 0.007
${{\rm{P.A.}}}_{{t}_{2}}$ (°)	339.99 ± 1.88

Note. ${t}_{1}\quad =$ UT 2014 November 5; ${t}_{2}\quad =\quad$ UT 2014 November 22. ${M}_{{Ys}}$ is not reported (or used in any photometric analysis) because the primary star has no reported measurements near 1 μm.

Download table as:  ASCIITypeset image

Henceforth, we will refer to the outer object as HD 7449B. Note that Roell et al. (2012) suggest that HD 7449 has a common proper motion companion at >2000 AU. The candidate companion was identified using the PPMXL proper motion catalog (Roeser et al. 2010). Examining the relevant images from PPMXL reveals that the object is actually one of the diffraction spikes and is therefore not a real astrophysical source. Therefore HD 7449 does not have any stellar companions at $\gt 2000\;{\rm{AU}}$.

Because we have detections of HD 7449B in both the visible and the near-infrared (NIR), we can use its colors and absolute magnitudes to constrain its spectral type, effective temperature (T_eff), and mass. To accomplish this, we compared its photometry to both known objects and to the low-mass stellar models of Kraus & Hillenbrand (2007) and Baraffe (Baraffe et al. 1998, 2002, 2015).

To create a comparative SED for HD 7449B, we began with the MagAO photometry for the primary and the outer companion. We used catalog 2MASS and SDSS photometry for HD 7449A and then used the color transformation relations in Carpenter (2001) to put the NIR photometry on the MKO system, which is comparable to the MagAO filters. Photometry for HD 7449B was then obtained by computing the magnitude differences relative to the primary. We used the Hipparcos parallax of 25.69 ± 0.48 mas (van Leeuwen 2007) to compute the absolute magnitudes and then converted each to $\lambda {F}_{\lambda }$ (e.g., see Faherty et al. 2013). Figure 4 shows the resulting SED for HD 7449B as well as the similarly computed SEDs of comparative M dwarfs from the 8 pc sample (Reid & Gizis 1997). The best matching SED corresponds to an M4.5, which also confirms that HD 7449B is at the distance to the primary (38.9 pc).

Zoom In Zoom Out Reset image size

To demonstrate that HD 7449A and B fall along the main sequence together (hence verifying that they are likely coeval), we constructed color–magnitude diagrams (CMDs) at several wavelengths from the visible to the NIR. We used the low-mass star Hipparcos sample, the NSTARS parallax sample, and the brown dwarf parallax sample from Dupuy & Liu (2012) and Faherty et al. (2012). Because these report photometry in the 2MASS system, we converted our MagAO photometry to 2MASS (assuming MKO comparable) using the Carpenter (2001) relations. All CMDs generally showed that the A and B components fall on the main sequence together, indicating that they are coeval and that the companion is not a background or foreground object. Figure 5 shows an example CMD.

Zoom In Zoom Out Reset image size

Using the Kraus & Hillenbrand (2007) models, comparing HD 7449B's colors and absolute magnitudes yielded a best-matching spectral type of M5, ${T}_{\mathrm{eff}}=3010\;{\rm{K}}$, and $M=0.15\;{M}_{\odot }$. Using the pre-2015 Baraffe models (Baraffe et al. 1998, 2002), and assuming the stellar age is between 1 and 3 Gyr, the colors were best matched by a 0.10 M_⊙, ${T}_{\mathrm{eff}}=2824\;{\rm{K}}$, 1 Gyr old star. The absolute magnitudes were best matched by a 0.15 M_⊙, ${T}_{\mathrm{eff}}=3161\;{\rm{K}}$, 1 Gyr old star. Using the 2015 Baraffe models, for stellar ages between 1 and 3 Gyr, the colors were best matched by a star with M = 0.20 ${M}_{\odot }$, ${T}_{\mathrm{eff}}=3261\;{\rm{K}}$, and the absolute magnitudes were best matched by a star with $M\quad =$ 0.09 ${M}_{\odot }$ and ${T}_{\mathrm{eff}}=2643\;{\rm{K}}$. Based on all of the above analysis, we classify HD 7449B (from photometry alone) as an M4.5 ± 0.5 with mass = $0.15\pm 0.05\;{M}_{\odot }$.

RVs have been obtained on HD 7449 for the past ∼15 years by HARPS and CORALIE (Dumusque et al. 2011), and by Magellan/MIKE and PFS (this work; see Table 1). To explain the periodic RVs (and the clear long-term trend), previous works (Dumusque et al. 2011; Wittenmyer et al. 2013) searched for solutions explained by one or more planets. We have the advantage that we know from direct imaging that the system contains an ∼M4.5 companion whose current projected separation is ≳21 AU. Can this companion explain the long-term trend and in doing so help revise the parameters of the inner planet(s)?

To test this, we first analyzed the RVs using log-likelihood periodograms (Baluev 2009; Anglada-Escudé et al. 2014) for preliminary period detection and confidence evaluation. Then we used a Bayesian Markov Chain Monte Carlo (MCMC) approach to produce posterior distributions of the allowed parameter values (Ford 2005). The likelihood function L contains the Keplerian model and a handful of nuisance parameters to account for the arbitrary zero-points of each RV instrument and the different levels of instrumental excess noise (also called jitter, which typically contains the contribution from stellar activity). The likelihood function is given by

$$\begin{eqnarray}&&L=\displaystyle \prod _{I}\displaystyle \prod _{i}^{{N}_{\mathrm{obs}}}{l}_{i,I}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{l}_{i,I}=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{1}{\sqrt{{\epsilon }_{i}^{2}+{s}_{I}^{2}}}\mathrm{exp}\left[-\displaystyle \frac{1}{2}\displaystyle \frac{{({v}_{i,I}-v(t,I))}^{2}}{{\epsilon }_{i,I}^{2}+{s}_{I}^{2}}\right]\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{v}_{i,I} & = & {\gamma }_{I}+\displaystyle \sum _{p}u({\hat{\kappa }}_{p};t)\\ & & +{\dot{v}}_{r}(t-{t}_{0})+\displaystyle \frac{1}{2}{\ddot{v}}_{r}{(t-{t}_{0})}^{2},\end{eqnarray} \tag{ 3 }$$

where i indexes the observations acquired with the Ith instrument, ${\epsilon }_{i,I}$ is the nominal uncertainty of each RV measurement, ${\gamma }_{I}$ and s_I are the zero-point and extra noise parameters (also called jitter) of each instrument, and the Doppler signal from a companion on the star is encoded in the model $u({\hat{\kappa }}_{p};t)$, which is a function of time t and the Keplerian parameters $\hat{{\kappa }_{p}}$. The Keplerian parameters of the pth companion in the system are: the orbital period P_p (in days), the semi-amplitude K_p (in m s⁻¹), the mean anomaly ${\mu }_{0,p}$ at the reference epoch t₀ (in degrees), the eccentricity e_p, and the argument of periastron ${\omega }_{p}$. The second and third terms in Equation (3) account for the possible presence of a long-period candidate whose orbit is only detected as a trend (acceleration, ${\dot{v}}_{r}$) plus some curvature (jerk, ${\ddot{v}}_{r}$). These two terms are especially important for the analysis that follows.

When performing the Bayesian MCMC analysis, one needs to specify some prior distributions for these parameters. In this paper, we use uniform distributions for the angles ${\mu }_{0}$ and ω, as any value would be equally likely a priori. Given that the objects involved are rather massive and the signals are large, we also allow for a uniform eccentricity distribution between [0,0.95). For K_p, ${\gamma }_{I}$, s_I, ${\dot{v}}_{r}$, and ${\ddot{v}}_{r}$, we assume unbound non-normalized uniform priors. While this can cause issues when normalizing the posterior, we are only using the MCMC analysis to sample the shape of the posterior, so precise values of the bounds and the normalization factors are unnecessary. Furthermore, because we will later try to constrain the signal with a period much longer than the span of the observations, the possible values of these three quantities will be correlated, so bound priors might eliminate many long period solutions that would otherwise be highly probable. For example, a large K in general requires a large γ unless the companion is precisely crossing the plane of the sky at t₀ (which is highly unlikely).

Regarding the prior on the period, in this work we assume that the prior is uniform in $1/P$ (equivalent to uniform in frequency). This is motivated by the following. When analyzing time series, the local solutions to periodic signals are approximately equally spaced in frequency. For example, if one produces a Lomb–Scargle periodogram of a time series and plots period versus power, one will quickly appreciate that the peaks become much broader with increasing period (Scargle 1982). However, when making the same plot in frequency, the peaks appear uniformly distributed over the possible frequencies. As discussed in Tuomi and Anglada-Escudé (2013), in Bayesian statistics the choice of the parameter automatically imposes implicit priors on all other alternative parameterizations. In the case of the period, a uniform prior at very long periods can outweigh the information content on the likelihood, producing a biased result. While this issue is not very severe for periods shorter than the time-span of the observations (typical RV planet search domains), the disrupting effects of the uniform prior become serious if one attempts to constrain very long orbits and becomes strongly dominated by the chosen period cut-off. On the other hand, the frequency parameter does not suffer from such singularities (all very long periods become packed in a single likelihood maxima close to 0) and preserves the role of the likelihood function as the most informative element in the posterior distribution.

Our MCMC algorithm is based on the one described in Ford (2005), which uses a Gibbs sampler with independent Gaussian jump functions for each parameter. For each parameter, the proposal function of our Gibbs sampler depends on a scale parameter that needs to be tuned to ensure acceptance rates between 10% and 30%. This is automatically done by tuning all the scale parameters until they reach the aforementioned acceptance rates (burn-in period). These samples typically amount for 10⁶ iterations and they are not used for the final MCMC analysis. In this paper we only focus on the detection and characterization of the two most significant signals in the RV data (HD 7994Ab and the long-period trend). While there have been other claims of possible candidates in the system (Dumusque et al. 2011; Wittenmyer et al. 2013), we suspect these were artifacts caused by sampling issues with the rather eccentric orbit of HD 7994Ab and the presence of the long-period parabolic trend.

Our first analysis consisted of a likelihood function with two Keplerian signals: one initialized at P ∼ 1200 days (which roughly corresponds to the preferred period for the most significant planet in Dumusque et al. 2011 and Wittenmyer et al. 2013), and the other one at P ∼ 8000 days, as suggested by the maximum likelihood periodograms in Figure 6. While a maximum likelihood orbit could be obtained with a second planet at ∼10,000 days (30 years), long MCMC runs indicated that the possible parameters of this object were heavily correlated. As a result, the parameter space was broadly unconstrained, making it difficult for the chains to achieve convergence even after 10⁸–10⁹ steps. Such strong degeneracy indicates that only a subset of the 5 Keplerian parameters can be constrained by the current data. As we will see later, the trend in the RV data can be well-described by two terms: an acceleration (linear trend) + jerk (curvature).

Zoom In Zoom Out Reset image size

As an attempt at better constraining the orbital elements of the outer companion, in our second analysis we included a Keplerian model for the outer companion's predicted orbital separation in order to use our direct imaging constraints (e.g., see Lucy 2014). Unfortunately, the orbital motion of the imaged companion was not very large between the direct imaging runs, and the imaging provides just two observables (projected separation in R.A. and decl.) while introducing three more free parameters (orbital inclination i, longitude of ascending node Ω, and the mass ratio between the companion and the primary star). As a result, these MCMC chains had even more difficulty converging to a meaningful equilibrium distribution (e.g., companion masses up to 100 ${M}_{\odot }$ and periods up to millions of years were consistent with the data).

Given that the entire RV data set (including HARPS, CORALIE, MIKE, and PFS RVs) can be well-fit by a simple parabola, for our third analysis we implemented a Doppler model containing a single inner planet plus a linear and quadratic term. In this case, the model contains a single Keplerian initialized at 1200 days (inner planet), plus the last two terms in Equation (3). Since ${\dot{v}}_{r}$ and ${\ddot{v}}_{r}$ are linear parameters, the MCMC quickly converged to the best-fit solution, which had an almost identical value of the likelihood function to the full two Keplerian solution attempted in Section 3.4. This means that the entire RV data set is best (and most simply) described by a single inner planet along with a long-term trend consisting of linear plus quadratic terms. We therefore use this final MCMC's results to constrain the parameters of the companions around HD 7449.

The posterior distributions of the inner planet's parameters are shown in Figure 7. Clearly, HD 7449Ab is eccentric (median e_b = 0.8), in agreement with Dumusque et al. (2011). To constrain the planet's mass (m), we used the distributions of K_b, P_b, and e_b, drew random Gaussian-distributed values for the stellar mass M_* having mean = 1.05 ${M}_{\odot }$ and standard deviation = 0.09 ${M}_{\odot }$ ¹⁵ , assumed ${M}_{*}\gg m$, and then solved for $m\mathrm{sin}{i}_{b}$ using the well-known relation

$$\begin{eqnarray}&&m\mathrm{sin}{i}_{b}={K}_{b}\sqrt{1-{e}_{b}^{2}}{M}_{*}^{2/3}{\left(\displaystyle \frac{{P}_{b}}{2\pi G}\right)}^{1/3}.\end{eqnarray} \tag{ 4 }$$

We find that the median $m\mathrm{sin}{i}_{b}$ = 1.09 M_J, which is in excellent agreement with preferred mass found by Dumusque et al. (2011). HD 7449Ab's parameters and their possible ranges are listed in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  HD 7449A Companion Parameters

	Mass	Period	a (AU)	e	ω (°)	i (°)
HD 7449Ab	$\gt {1.09}_{-0.19}^{+0.52}$ MJ	${1270.5}_{-12.1}^{+5.92}$ days	${2.33}_{-0.02}^{+0.01}$	${0.80}_{-0.06}^{+0.08}$	$-{25.2}_{-5.22}^{+6.87}$	unconstrained
HD 7449B	${0.23}_{-0.05}^{+0.22}$ ${M}_{\odot }$	${65.7}_{-56}^{+227}$ years	${17.9}_{-12.9}^{+32}$	unconstrained	unconstrained	${59.7}_{-25.8}^{+20.1}$

Note. All uncertainties correspond to symmetric 68% confidence intervals around the median values.

Download table as:  ASCIITypeset image

Here we develop and apply a new statistical procedure, expanding on the one developed in Torres (1999), that uses the slope and quadratic terms discussed in Section 3.6 to tightly constrain the outer companion's properties.

For the case of an imaged companion producing a long-period RV trend, Torres (1999) formulated a numerical Monte Carlo approach to marginalize over unknown parameters under some uninformative priors. It is based on using the fact that the linear trend observed in a Doppler curve of the primary star can be written in terms of the mass of the long-period companion only (M_B), and that the observed separation at a given epoch t₀ can be written as a function of M_B, the direct imaging observables, and a function that can be easily marginalized over the unknown orbital parameters (P_B, e_B, ${\omega }_{B}$, ${\mu }_{0,B}$, i_B, ${{\rm{\Omega }}}_{B}$). The method described in Torres (1999) uses only the measured linear part of the trend and produces a distribution of possible masses. We now develop a method that exploits the second derivative of the RV (the jerk), which allows us to obtain a probability distribution for the companion's orbital period as well. In this section, all the quantities refer to the secondary companion, so we will avoid using sub-indices for clarity, except for the period P_B and mass M_B of the secondary.

We begin with Equations (3) and (5) from Torres (1999) to write the derivative of the RV of the primary component ${\dot{v}}_{r}$ as

$$\begin{eqnarray}&&{\dot{v}}_{r}=\displaystyle \frac{{{GM}}_{B}}{{{\ell }}^{2}}{\rm{\Psi }},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{\rm{\Psi }} & = & {[(1-e)(1+\mathrm{cos}E)]}^{-1}(1-e\mathrm{cos}E)\mathrm{sin}i\\ & & \times (1-{\mathrm{sin}}^{2}(\nu +\omega ){\mathrm{sin}}^{2}i)(1+\mathrm{cos}\nu )\mathrm{sin}(\nu +\omega ),\end{eqnarray} \tag{ 6 }$$

where ℓ is the projected separation between the primary star and the companion in physical units (e.g., mks), and Ψ is a rather intricate function that encapsulates all the orbital elements to be marginalized (time-dependencies included). Equation (5) can be evaluated by solving Kepler's equation

$$\begin{eqnarray}&&E-e\mathrm{sin}E=\mu \end{eqnarray} \tag{ 7 }$$

to obtain the eccentric anomaly E, and then using

$$\begin{eqnarray}&&\mathrm{tan}\displaystyle \frac{\nu }{2}=\sqrt{\displaystyle \frac{1+e}{1-e}}\mathrm{tan}\displaystyle \frac{E}{2}\end{eqnarray} \tag{ 8 }$$

to derive the true anomaly ν. To account for all possible combinations of periods and orbital phases, Torres (1999) realized that the mean anomaly $\mu =\frac{2\pi }{{P}_{B}}t+{\mu }_{0}$ could be assumed to be uniformly distributed in $(0,2\pi ]$. That is, irrespective of the values of the observation time t and P_B, μ still can be assumed to have any orbital phase because ${\mu }_{0}$ can also have any value between 0 and $2\pi$.

To use the information on the quadratic term in the RV curve, we first need to compute the second derivative of the RV. To do this efficiently, it is enough to realize that all the time dependence is included in μ. Therefore, we can apply the chain rule and the fact that $d\mu /{dt}=2\pi /{P}_{B}$ to obtain

$$\begin{eqnarray}&&{\ddot{v}}_{r}=2\pi \;\displaystyle \frac{G}{{{\ell }}^{2}}\displaystyle \frac{{M}_{B}}{{P}_{B}}{{\rm{\Psi }}}^{\prime }\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}^{\prime }=\displaystyle \frac{d{\rm{\Psi }}}{d\mu }.\end{eqnarray} \tag{ 10 }$$

This is an important result because we have found that ${\ddot{v}}_{r}$ is proportional to ${M}_{B}/{P}_{B}$, and all the dependencies can again be marginalized by evaluating ${{\rm{\Psi }}}^{\prime }$ (which is a function of time because it depends on E). While an analytic expression for ${{\rm{\Psi }}}^{\prime }$ could be derived, it is far simpler (and requires fewer operations) to compute this numerically. We found that a simple two point formula with an infinitesimal increment of 10⁻⁴ radians for μ works to sufficient precision. By rearranging terms in Equation (5) and combining Equation (5) with Equation (9), we find

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{B}}{{M}_{\odot }}=5.341\times {10}^{-6}\;{\dot{v}}_{r}{\left(\displaystyle \frac{\rho }{{\rm{\Pi }}}\right)}^{2}\displaystyle \frac{1}{{\rm{\Psi }}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{B}}{{yr}}=2\pi \;\displaystyle \frac{{\dot{v}}_{r}}{{\ddot{v}}_{r}}\displaystyle \frac{{{\rm{\Psi }}}^{\prime }}{{\rm{\Psi }}},\end{eqnarray} \tag{ 12 }$$

where ρ and Π are the projected separation and parallax in arcseconds, respectively. Equation (11) was already derived in Torres (1999). The additional relation that we present here (Equation (12)) can be used to constrain a companion's period using the same observables and marginalization method outlined in Torres (1999). The numerical factors in the equations come from the numerical substitution of the gravitational constant G, and the choice of units in Torres (1999), which assumes that ${\dot{v}}_{r}$ is in m s⁻¹ yr⁻¹, and ${\ddot{v}}_{r}$ is in m s⁻¹ yr⁻². Note that Equations (5) and (12) assume that we can produce a Taylor expansion of the RV near the epoch of the direct image(s). The most straightforward way to impose this condition is to set the reference epoch in Equation (3) to ${t}_{0}={t}_{{\rm{image}}}$, thus deriving consistent MCMC samples for ${\dot{v}}_{r}$. This consideration is unnecessary when no curvature is detectable in the RV curve, as the first derivative of the RV is then independent of time.

With Equations (5) and (12) in hand, we set out to constrain the mass and period of the outer companion. All the quantities have uncertainties, including ${\dot{v}}_{r}$, ${\ddot{v}}_{r}$, ρ, and Π. To account for these, we applied an additional refinement to the marginalization procedure of Torres (1999). That is, in addition to drawing random values for e, μ, and $\mathrm{sin}{i}_{B}$, we also drew randomly generated values of the observables: ${\dot{v}}_{r}$ and ${\ddot{v}}_{r}$ pairs were drawn from randomly selected states of the MCMC, and Gaussian distributions consistent with the error on parallax and the error on projected separation (Table 2) were used to generate plausible pairs of Π and ρ, respectively. We know that the outer companion must be less massive than the primary (or it would be a known visual binary), so we excluded all parameters that corresponded to mass >1 ${M}_{\odot }$. Subsamples of the resulting distributions for the outer companion's mass, period, and inclination are shown in Figure 8. Based on these distributions, the median mass of HD 7449B is 0.23 ${M}_{\odot }$, and the mass of HD 7449B is larger than 0.17 ${M}_{\odot }$, with a 99% probability, consistent with our constraints from photometry. The median period is 65.7 years, corresponding to a semimajor axis of 17.9 AU. The inclination distribution is broad and has a median value of i = 597  because it is mostly inherited from the uniform distribution in the Monte Carlo generated test values. However, because we exclude masses larger than 1 ${M}_{\odot }$, values of i_B < 84  are also excluded with 99% probability, thus ruling out strictly face-on orbits. We do not show the distributions for e and ω because they were completely unconstrained. Table 3 lists the outer companion's constrained parameters and ranges.

Zoom In Zoom Out Reset image size

Based on the above analysis, we were able to constrain the outer companion's mass, period, and inclination, but not its eccentricity. To get a sense of the allowed ranges, we explored the dynamical stability of the system. We used the MERCURY integration package (Chambers 1999) with a Bulirsch-Stoer integrator and simulated 100 different realizations of the system for 1 Gyr. The initial semimajor axis and eccentricity of the inner planet (HD 7449Ab) were held fixed at 2.32 AU, and 0.78, respectively, while its mass was set to 1 M_J.¹⁶ The outer companion's semimajor axis was fixed at 18 AU and we assumed near-coplanarity such that its initial inclination relative to the planet was randomly drawn from values between 0 and 1°.¹⁷ In each of the 100 simulations, the companion's eccentricity was varied between 0 and 1 in increments of 0.01. Its mass was set to 0.17 ${M}_{\odot }$. For both the outer companion and the planet, the arguments of pericenter, longitudes of ascending node, and mean anomalies were all drawn randomly from a uniform distribution in each simulation.

In this case, the outer companion's critical eccentricity e_crit (the eccentricity above which the planet's orbit becomes unstable) was 0.45. Based on this initial result, we performed additional simulations in which the planet and outer companion had mutual inclinations ${\rm{\Delta }}i$ = 30°, 60°, 90°, 120°, 150°, and 180°. For the 90° case, the planet was never stable regardless of the outer companion's eccentricity. The critical eccentricities for the other inclinations were (in ascending order of mutual inclination) 0.4, 0.4, 0.4, 0.42, and 0.5, respectively. Based on these results, the outer companion's eccentricity is constrained to be $\lesssim 0.5$. Figure 9 summarizes the results of our stability analysis.

Zoom In Zoom Out Reset image size