Evidence That the Directly Imaged Planet HD 131399 Ab Is a Background Star

HD 131399 A was observed with GPI at two epochs, 2017 February and 2017 April, as part of the GPI Exoplanet Survey (GS-2015B-Q-501). Three data sets were obtained on consecutive nights in 2017 February, with a total on-source integration time of 1.87 hr at $K{1}_{\mathrm{GPI}}$ (${\lambda }_{\mathrm{eff}}=2.06$ μm), 1.38 hr at H (${\lambda }_{\mathrm{eff}}=1.64$ μm), and 1.60 hr at J (${\lambda }_{\mathrm{eff}}=1.23$ μm). An additional data set was obtained on 2017 April 20 at H with an on-source integration time of 1.03 hr. Each data set was obtained in the spectral coronagraphic mode of the instrument.

To create spectral data cubes, the raw data were reduced with the GPI Data Reduction Pipeline v1.4.0 (DRP; Perrin et al. 2014, 2016), which subtracts the dark current, removes the microphonic noise (Chilcote et al. 2012; Ingraham et al. 2014b), and identifies and removes bad pixels. Instrument flexure is compensated for using observations of an argon arc lamp taken immediately prior to each sequence at the target elevation (Wolff et al. 2014). Microspectra are then extracted to create 37-channel data cubes (Maire et al. 2014), which are corrected for any remaining bad pixels and finally for distortion (Konopacky et al. 2014). The last step consists of measuring in each image the location of the four satellite spots—attenuated replicas of the central point-spread function (PSF) created by a diffraction grating in the pupil plane—to accurately measure the position and flux of the central star during the sequence (Wang et al. 2014). The position of each satellite spot flux is written in the header so that it can be used for calibration.

Further processing to remove the stellar PSF and extract the astrometry and spectrophotometry of HD 131399 Ab was performed using two different pipelines to mitigate biases and systematics introduced by the data processing.

In the first pipeline, a Fourier high-pass filter with a smooth cutoff frequency of four spatial cycles was applied to each image. The speckle field was then estimated and subtracted using the classical ADI algorithm (cADI; Marois et al. 2006; following the definition of Lagrange et al. 2010 as a median combination) for each sequence in each wavelength slice, which was then rotated to align north with the vertical axis and averaged over the sequence. Broadband images were further created from the stack of the individual slices, examples of which are shown in Figure 1. The astrometry and broadband contrasts of HD 131399 Ab were extracted in each data set from the broadband images using the negative simulated planet technique (Lagrange et al. 2010; Marois et al. 2010). A template PSF of HD 131399 Ab was created from the temporal and spectral average of the four satellite spots. The template was injected in the raw data cubes at a trial position but opposite flux of HD 131399 Ab, and the same reduction as for the original set was executed. The process was iterated over these three parameters (separation, position angle, flux) to minimize the integrated squared pixel noise in a wedge of 3 × 3 FWHM centered at the trial position. The minimization was performed with the amoeba-simplex optimization algorithm (Nelder & Mead 1965) and provided the best-fit broadband contrast and position. Uncertainties on HD 131399 Ab location and contrast were calculated by injecting independently 20 positive templates at the same separation and contrast as HD 131399 Ab but different position angles. The fitting procedure was repeated for each simulated source and the measurement errors obtained from the statistical dispersion on the three parameters. Finally, the contrasts—and associated measurement errors—in individual slices in each set were then extracted following the same procedure at the best-fit position, which is fixed, and varying only the flux of the template, which is built for each wavelength from the corresponding satellite spots.

Zoom In Zoom Out Reset image size

The second pipeline used pyKLIP (Wang et al. 2015), an open-source Python implementation of the Karhunen-Loève Image Projection algorithm (KLIP; Soummer et al. 2012). Before PSF subtraction, the images were high-pass filtered using a seven-pixel FWHM Gaussian filter in Fourier space to remove the smooth background. KLIP was run on a 22-pixel-wide annulus centered on the location of the source. To build the model of the stellar PSF, we used the 150 most-correlated reference images in which HD 131399 Ab moved at least a certain number of pixels due to ADI and SDI observing methods (the exclusion criteria). Since we will forward-model the PSF of the planet, including the effects of self-subtraction, we use an aggressive exclusion criteria of 1.5 pixels for all wavelengths except the J-band, where we found that using images very close in time most accurately modeled the speckles and thus a 0.2 pixel exclusion criteria worked best. As the source is far from the star and thus from the majority of the speckle noise, we used only the first five KL basis vectors to reconstruct the stellar PSF. All images were then rotated to align north up and collapsed in time and wavelength, resulting in one 2D image per epoch. The astrometry and broadband photometry were measured from these images using the Bayesian KLIP-FM Astrometry (BKA) technique (Wang et al. 2016) that is implemented in pyKLIP. In BKA, we concurrently forward-model the PSF of HD 131399 Ab during KLIP. To do this, we used the average of the satellite spots to model the instrumental PSF at each wavelength, and we assumed HD 131399 Ab had a spectral shape that was the same as HD 131399 A. As noted in Wang et al. (2016), spectra differing by even 20% did not affect the astrometry, so we did not require a precise input spectral template for our forward model. After generating the forward model, we used the affine-invariant sampler implemented in emcee (Foreman-Mackey et al. 2013) to compute the posterior distribution of the location and flux of HD 131399 Ab. Our MCMC sampler used 100 walkers, each iterating for 800 steps after 300 steps were discarded as the "burn in." To obtain accurate uncertainties, the residual speckle noise in the image was modeled as a Gaussian process with a spatial correlation described by the Matérn covariance function. We adopt the 50th percentile values as the position of HD 131399 Ab and the 16th and 84th percentile values as the 1σ uncertainty range. To obtain the spectrum of HD 131399 Ab in each filter, we performed a PSF subtraction with KLIP that only used ADI to model the stellar PSF, allowing us to forward-model the PSF of HD 131399 Ab without any spectral dependencies. Then, we modified BKA to run independently on each spectral channel to obtain the flux and the uncertainty on the flux at each wavelength. As the planet has significantly lower signal-to-noise ratio (S/N) in each spectral channel than in a collapsed broadband image, we restricted the position of the planet to be within 0.1 pixel of the position we measured in the broadband data.

Astrometric calibration was obtained with observations of the ${\theta }^{1}$ Ori B field and other calibration binaries following the procedure described in Konopacky et al. (2014) and used to convert the detector positions into on-sky astrometry. The astrometric error budget consists of the following added in quadrature: the measurement errors described previously; a star registration error of 0.7 mas from Wang et al. (2014); a plate scale error of 0.007 mas lenslet⁻¹; and position angle offset error of 0.13 deg, the last two from Konopacky et al. (2014).

The raw astrometric and photometric measurements from the two pipelines (i = 1, 2) agreed very well to better than 1σ at each epoch. The pairs (${x}_{i},{\sigma }_{i}$) from the two pipelines for each data set were combined with a weighted average ${x}_{\mathrm{tot}}={\sum }_{i}{w}_{i}{x}_{i}/{\sum }_{i}{w}_{i}$, where ${w}_{i}=1/{\sigma }_{i}^{2}$. The measurement errors were computed as ${\sigma }_{\mathrm{tot}}=\sqrt{{\sum }_{i}{\sigma }_{i}^{2}{w}_{i}/{\sum }_{i}{w}_{i}}$ since they are not independent. The systematic errors (registration, calibration) were then added in quadrature to calculate the final astrometric uncertainties.

Photometric measurements from different epochs (j = 1, 2) were also combined with the same weighted mean, but the errors were computed as ${\sigma }_{\mathrm{tot}}=\sqrt{1/{\sum }_{j}{w}_{j}}$ since they are independent. Finally, the systematic uncertainties of the star-to-satellite-spot ratios (0.03 mag in the J band, 0.06 mag in the H band, and 0.07 mag in the K1 band, Maire et al. 2014) were added in quadrature to the final contrast errors. The spectrum was then obtained by multiplying the contrasts with the spectrum of the central star (see Section 3.1).

The S/Ns for each data set were computed using the pyKLIP implementation of the Forward Model Matched Filter (FMMF) algorithm (Ruffio et al. 2017), using the stellar spectrum of HD 131399 A as the spectral template in the matched filter. Like the two pipelines to extract astrometric and photometric data, FMMF similarly utilizes forward modeling of point sources through the PSF subtraction process for the data analysis, but is better optimized for planet detection. Thus, FMMF produces S/Ns that are comparable or slightly better than the S/Ns inferred from the astrometric for photometric errors.

All measurements are discussed in Sections 3.2 and 4.

2.2.1. Reanalysis of Public Data

Four epochs of observations were obtained with SPHERE by W16 between 2015 June and 2016 May, all of which are publicly available on the ESO archive.³⁷ We downloaded the data as well as the associated raw calibration files. Briefly, the HD 131399 system was observed with the IRDIFS_EXT mode using simultaneously the Integral Field Spectrograph (IFS; Claudi et al. 2008) instrument in spectroscopic mode over 0.95–1.65 μm (YJH) and the Infra-Red Dual-beam Imaging and Spectroscopy (IRDIS; Dohlen et al. 2008) instrument in dual-band imaging mode (DBI; Vigan et al. 2010) at $K{1}_{\mathrm{SPH}}$ (${\lambda }_{\mathrm{eff}}=2.10$ μm) and $K{2}_{\mathrm{SPH}}$ (${\lambda }_{\mathrm{eff}}=2.25$ μm), with all SPHERE filter profiles being different from those of GPI (see Section 3.1 and Figure 8). The IRDIS detector was dithered on a 4 × 4 pattern. A total of 0.44 hr, 0.75 hr, 0.50 hr, and 0.50 hr were obtained on the IFS on 2015 June 12, 2016 March 06, 2016 March 17, and 2016 May 07, respectively, and 0.43 hr, 0.56 hr, 0.50 hr, and 0.50 hr on IRDIS, the difference between the two detectors being due to readout overheads. Each observing sequence started and finished with a brief "star-center" coronagraphic sequence in which four satellite spots are created from a periodic modulation introduced on the deformable mirror, the barycenter of these spots being used to measure the position of the star behind the focal plane mask during the sequence. In practice, the star position is very stable (Zurlo et al. 2016; Vigan et al. 2015). A brief off-axis ($\simeq 0\buildrel{\prime\prime}\over{.} 4$) "flux" sequence with the neutral density filter ND3.5 (attenuation factor from 4 × 10² to 2 × 10⁴ dependent on wavelength) was then executed to obtain a template and the flux of the target PSF. The on-axis coronagraphic sequence was then carried out. Calibration data were obtained during the following days: darks, detector flat fields, integral field unit flat (broadband lamp image to register the IFS microspectra), and a wavelength calibration frame.

IFS data processing. The raw data and calibration files were reduced using the SPHERE IFS preprocessing tools v1.2³⁸ (Vigan et al. 2015), which make use of custom IDL routines and the ESO Data Reduction and Handling (DRH) package v22.0 (Pavlov et al. 2008). These tools were updated with the latest calibration values provided by Maire et al. (2016) and the ESO SPHERE user manual 7th edition³⁹ : instrument angle updates (pupil offset of 135.99 deg, and IFS angle offset of −100.48 deg), the IFS anamorphism correction (1.0059 along the horizontal direction, 1.0011 along the vertical direction), and the parallactic angle correction , a small factor to correct the parallactic angle calculation for a missynchronization between the VLT and SPHERE internal clock that affects data taken before 2016 July 13. Additionally, the tools were updated to process the entire field of view (it was originally cropped by five pixels on the edges). The preprocessing tools used the DRH package to create the master darks, bad pixel maps, the microspectra position map, the IFU flat field, and the wavelength calibration file. Detector flats were created with a custom IDL routine. The data preprocessing was then executed by a custom IDL routine, which subtracts the dark current, removes the bad pixels, and corrects for cross-talk. This was followed by processing through the DRH, which corrects for flat-fielding and extracts the microspectra to create 39-channel data cubes. The 3D data cubes were then digested by a custom IDL routine to remove the remaining bad pixels, to correct for the anamorphism, to register the spot locations in the star-center frames, to align the coronagraphic and the off-axis PSF frame at the center, and to recalibrate the wavelengths.

IRDIS data processing. A custom set of tools to reduce IRDIS DBI data was developed following the IFS philosophy, combining both DRH and IDL routines. IRDIS DBI raw data are made from images in two side-by-side quadrants, being associated with the K1 (left) and K2 (right) filters. The DRH first created the master darks, flat fields, and associated bad-pixel maps. Our IDL routine then performed the dark current subtraction, flat-field division, bad-pixel removal, vertical anamorphism correction by a factor of 1.006 (Maire et al. 2016), and parallactic angle calculation and correction by the factor. For each image, the two quadrants were separated at the end to create a master data cube for each filter. The locations of the satellite spots and frame registration, taking into account the dithering offset from the header keywords, were performed as for the IFS data as final processing steps.

Similarly to the GPI data, the speckle field in both IRDIS and IFS data cubes was removed using the two postprocessing pipelines as described in Section 2.1. Final broadband images at K1, K2, and YJH, created from the stack of the 39-channel IFS data cubes, are shown for each epoch in Figure 2. The position, contrast, and measurement uncertainties of HD 131399 Ab were also obtained using the same techniques as for GPI; the PSF templates for the IRDIS and IFS data were built from the unsaturated off-axis images of the star. The astrometric calibrations of the plate scale and position angle for both instruments are given by Maire et al. (2016) and the ESO SPHERE user manual 7th edition to convert the on-chip measurements into on-sky positions. These calibration values have been stable since the commissioning of the instrument, when taking into account the missynchronization correction between the SPHERE and VLT clocks. The final astrometric error budget consists of the following added in quadrature: the measurement errors described in Section 2.1; a star registration error of 0.1 px (Vigan et al. 2015; Zurlo et al. 2016); a plate scale error of 0.02 mas lenslet⁻¹ (IFS) and 0.021 mas px⁻¹ (IRDIS); a pupil angle offset error of 0.11 deg; a position angle offset error of 0.08 deg; and an IFS angle offset error of 0.13 deg.

Zoom In Zoom Out Reset image size

The spectrophotometric and astrometric measurements from the two pipelines agreed very well to better than $1\sigma$ at each epoch and were combined following the procedure used for the GPI data (see Section 2.1). The S/Ns for all of the data sets, except the 2016 May 7 IFS data, were also computed using the same FMMF algorithm as the GPI data. Due to the short amount of time HD 131399 Ab stays on the chip and some artifacts on the edge of the images, the 2016 May 7 IFS data seemed to be a pathological data set for the FMMF algorithm. Instead, for this data set, we computed the S/Ns by cross-correlating each broadband-collapsed image with a Gaussian PSF and comparing the peak of the cross-correlation of Ab with the standard deviation of the cross-correlation of the noise at the same separation.

Figure 3 shows the spectrum of HD 131399 Ab extracted from each epoch of IFS data. The spectra are very noisy because HD 131399 Ab is barely detected in individual slices, especially in the Y and J bands. As reported in Table 1, the source lies on the detector a small fraction of the total time in three data sets (as low as 30%), lies very close to the edge of the detector another significant portion, particularly in the 2016 May data set, and falls off the chip up to 39% of the time in our reduced IFS images. Ultimately, this reduced effective observing time strongly affects the data quality. The continuum and flux are nevertheless consistent between the different epochs, except between the J and H bands, where the atmospheric transmission is low. However, the third epoch strongly differs from the other three in the H band, exhibiting a steep slope with a peak at $1.61-1.63$ μm. To assess this feature, the 2016 March 17 data were reduced using LOCI, and the spectrum was extracted following the same procedure as for the cADI/pyKLIP analysis. In both cases, the slope and peak were both recovered. A visual inspection of the reduced data cubes reveals the presence of a speckle very close to HD 131399 Ab, which becomes more prominent in the 1.61 and 1.63 μm channels (see Figure 4). More aggressive high-pass filters and algorithm parameters are not able to suppress this speckle. We therefore propose that the peak of the spectrum of HD 131399 Ab in the 2016 March 17 may be biased by this speckle, particularly in less aggressive reductions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To mitigate this effect and also to improve the S/N of the spectrum, we followed the strategy of W16 and combined the four data sets using both pipelines.

In the first pipeline, the cADI flux loss ($\simeq 5 \%$) was compensated for in each data set by injecting and reducing simulated sources at the same separation as HD 131399 Ab, but at 20 other position angles. A stamp of 30 × 30 pixels centered at the measured position of HD 131399 Ab was then extracted in the cADI-reduced image at each epoch. The stamps of the four epochs were averaged for each wavelength slice. To extract the flux at each wavelength from the combined data, we created a forward model of the PSF of HD 131399 Ab. At each epoch and wavelength slice, the off-axis PSF was injected in a noise-free data cube at the separation and position angle of HD 131399 Ab and reduced using the parallactic angle exploration of each epoch with cADI. Stamps of the model were then extracted and combined similarly. The combined model was used to fit the flux of HD 131399 Ab using the amoeba-simplex minimization procedure. To estimate the uncertainties, the exercise (injection of simulated sources in the raw data and forward model computation) was repeated at the same separation but at 20 different position angles. The statistical dispersion of the extracted fluxes was used as the uncertainty in the spectrum at each wavelength.

In the second pipeline, we extracted from the pyKLIP-reduced data and forward-modeled PSF a 11 × 11 pixel stamp centered at the location of HD 131399 Ab at each epoch. The stamps of both the data and forward model were averaged over the four epochs at each wavelength slice, resulting in one stamp of both the data and forward model at each wavelength. Then, we follow the same BKA technique as before to measure the flux and quantify the uncertainties in each wavelength channel.

The spectra were then combined in the same way as discussed previously. The results are discussed in Section 3.2.

Astrometry and photometry of HD 131399 B in the IRDIS K1 and K2 unsaturated off-axis images were obtained using the same technique used for the NIRC2 data and are discussed in Sections 3 and 4.

2.2.2. New Observations

HD 131399 was observed on 2017 March 15 (098.C-0864(A), PI: Hinkley) with SPHERE IRDIS in dual polarimetric imaging (DPI) mode at J (${\lambda }_{\mathrm{eff}}=1.23$ μm) as part of a program to measure the polarization of directly imaged planets. The same "star-center," "flux," and "coronagraphic" sequences, as were executed for the public DBI observations described in Section 2.2.1, were carried out in this program, for a total on-source integration time of 0.36 hr. Calibration data were obtained on subsequent days, following the standard calibration plan for the instrument.

The raw data were reduced following the same procedure as the public IRDIS DBI data. However, since the data were taken in DPI mode, images in the two quadrants, corresponding to two orthogonal polarization states, were summed to create total-intensity images. PSF-subtracted (see Figure 5) photometric and astrometric measurements were also obtained using the two postprocessing pipelines and the same parameters as described previously. Finally, the measurements from these two pipelines were combined with a weighted mean as for the other data sets and are reported in Table 5 and in Table 6.

Zoom In Zoom Out Reset image size

HD 131399 A was observed with the narrow camera of Keck/NIRC2 in the L' filter (${\lambda }_{\mathrm{eff}}=3.72$ μm) serving as its own natural guide star on consecutive nights 2017 February 7 and 8. We used only the February 8 data in our final analysis because high winds and poor seeing degraded the quality of the February 7 data. This resulted in 166 exposures of 0.9 s and 30 coadds each for a total integration time of 1.25 hr. The 400 mas diameter coronagraph mask occulted the star in all exposures, and the instrument was in vertical angle mode to enable ADI. The raw data were reduced with a custom set of tools that subtracts dark current and thermal background and then aligns all frames to a common star position.

To recover HD 131399 Ab, we subtracted the stellar halo and speckle pattern using a customized LOCI algorithm ("locally optimized combination of images"; Lafrenière et al. 2007). We tested various levels of algorithm aggressiveness and present here a compromise between noise suppression and astrophysical source throughput, with LOCI parameter values of ${N}_{\delta }=0.3$, W = 10 px, dr = 10 px, g = 0.9, and N_a = 10 following the conventional definitions in Lafrenière et al. (2007). Speckle suppression in this data set particularly benefited from temporal proximity of reference images (i.e., small ${N}_{\delta }$), possibly due to high air mass and varying seeing conditions diminishing PSF stability. The PSF-subtracted frames were rotated to place north up and collapsed into a final median image (see Figure 6, left). We also performed a separate reduction using pyKLIP on the same aligned frames. The algorithm divided images into annuli that were 20 pixels wide radially and further divided into 10 azimuthal subsections each. To build the model of the stellar PSF, we used the first 50 KL basis vectors of the 200 most-correlated reference images where HD 131399 Ab moved at least three pixels due to the ADI observing method (Figure 6, right).

Zoom In Zoom Out Reset image size

In neither reduction was a source detected at the location of Ab with greater than 3σ confidence over the background noise levels (see Figure 6). Therefore, we report only a lower limit of 11.10 mag for its L' contrast.

HD 131399 B and C are detected in individual images in which HD 131399 A is unocculted and unsaturated, so we performed astrometry on brighter component B as an independent confirmation of our SPHERE astrometry. To locate A, we fitted it with a bivariate Gaussian function using a least-squares minimization. We then jointly fitted B and C using the PSF of A as a template for a least-squares minimization. We repeated this process for six images divided between two dither positions, and we report in Section 4 the mean separation and PA of B from those fits. The measurement errors were estimated as the standard deviation of the separation and PA across the six images. The final astrometric uncertainties were calculated as the quadrature sum of these measurement errors, the star registration error estimated at 5 mas, and the plate scale error of 0.004 mas pixel⁻¹ and position angle offset error of 0.02 deg (Service et al. 2016).

A flux-calibrated spectrum of the primary was required to convert the measured contrast between HD 131399 A and Ab within the SPHERE and GPI data sets. As no near-IR spectrum of HD 131399 A was available within the literature, we used a stellar evolutionary model and a grid of synthetic stellar spectra to fit the observed SED of HD 131399 A. From this fit we estimated both the spectrum of the star and synthetic photometry within the GPI and SPHERE passbands. The properties of HD 131399 A and Ab are given in Table 2.

Optical and near-infrared photometry were found in the literature for a number of systems: Tycho (${B}_{T}{V}_{T};$ Høg et al. 2000), Hipparcos (H_p; ESA 1997), and 2MASS (${{JHK}}_{s};$ Skrutskie et al. 2006). Optical color indices in the Strömgren uvby (Hauck 1986) and Geneva⁴⁰ systems were also found (Mermilliod et al. 1997). An uncertainty of 0.1 mag was assumed for these color indices as none were presented within the literature. As the angular separation between HD 131399 A and BC is comparable to the angular resolution of the telescopes used to obtain these photometric measurements, the measures reported within these catalogs are of the blended system rather than of HD 131399 A. At shorter wavelengths, the contrast between HD 131399 A and the BC pair is large enough that the faint pair has a negligible impact on the optical photometry of the system. At longer wavelengths, this effect becomes significant, approximately 10% at K. To account for this, we simultaneously fit the combined flux of the three stars using the photometric measurements of the system described previously and apparent magnitudes of the BC pair obtained from the literature.

We used the emcee parallel-tempered affine-invariant MCMC sampler to fully explore the parameter space and estimate uncertainties on the near-IR spectrum of HD 131399 A. At each step within a chain, an age t, parallax π, mass for each component ${M}_{{\rm{A}}}$, ${M}_{{\rm{B}}}$, ${M}_{{\rm{C}}}$, and extinction A_V were selected. We used a Gaussian prior for age (16 ± 7 Myr) and a Gaussian (10.20 ± 0.70 mas) multiplied by a ${\pi }^{-4}$ power law—to account for a uniform space density of stars as expected at the distance to HD 131399 Ab—as the prior for parallax. The prior on the three masses was based on the Kroupa (2001) initial mass function. Age and mass were converted into an effective temperature (${T}_{\mathrm{eff}}$) and surface gravity ($\mathrm{log}g$) using the MIST evolutionary models (Choi et al. 2016; Dotter 2016). Given the rapid rotation seen for young early-type stars (e.g., Strom et al. 2005), we used the evolutionary models that incorporated stellar rotation ($v/{v}_{\mathrm{crit}}=0.4$). A solar metallicity was assumed, consistent with the observed metallicity of other stars within the ScoCen association (Bubar et al. 2011).

Synthetic photometry and color indices were computed from a BT-NextGen model atmosphere (Allard et al. 2012)⁴¹ of the appropriate ${T}_{\mathrm{eff}}$ and $\mathrm{log}g$, scaled by the ${R}^{2}/{d}^{2}$ dilution factor, where R is the radius of the star computed from M and $\mathrm{log}g$, and $d=1/\pi$ is the distance to the star. Model atmospheres at temperatures and surface gravities between grid points were estimated using a linear interpolation of the logarithm of the flux. These synthetic spectra were first reddened using the selected A_V value and the Cardelli et al. (1989) extinction law, and then convolved with the throughput of each filter to obtain synthetic photometry. Filter transmission profiles and zero points were obtained from Mann & von Braun (2015) for the optical filters, and from Cohen et al. (2003) for the 2MASS filters. A probability ($\mathrm{ln}p=-{\chi }^{2}/2$) was calculated at each step by comparing the synthetic magnitudes and color indices for the blended system to the observed values, the synthetic magnitudes of the B and C components to the $K{1}_{\mathrm{SPH}}$ contrasts (the SPHERE/IRDIS filters are described later in this section) given in W16 (${\rm{\Delta }}K{1}_{\mathrm{SPH}}=1.86\pm 0.10$ mag and 3.86 ± 0.10 mag for B and C, respectively), and the apparent H_p magnitude for the blended BC pair of 11.161 ± 0.187 mag reported in the Catalog of the Components of Double and Multiple Stars (CCDM; Dommanget & Nys 2002).

We initialized 512 walkers at each of 16 different temperatures to ensure the parameter space was fully explored; lower temperatures sample the posterior distribution, while higher temperatures fully explore the prior distributions. Each walker was advanced for 1000 steps as an initial burn-in stage, and then advanced for a further 9000 steps to fully sample the posterior distribution for each parameter. The median and 1σ range calculated from the posterior distribution of the six fitted parameters (t, π, ${M}_{{\rm{A}}}$, ${M}_{{\rm{B}}}$, ${M}_{{\rm{C}}}$, A_V), and that of the derived ${T}_{\mathrm{eff}}$ and $\mathrm{log}g$ for each component, are given in Table 3.

Table 3.  Stellar Parameters Derived from SED Fit

Property	Unit	HD 131399 system
t	Myr	${21.9}_{-3.8}^{+4.1}$
π	mas	${9.27}_{-0.37}^{+0.33}$
d	pc	${107.9}_{-3.7}^{+4.5}$
AV	mag	0.22 ± 0.09
		A	B	C
M	M⊙	${2.08}_{-0.11}^{+0.12}$	0.95 ± 0.04	0.35 ± 0.04
${T}_{\mathrm{eff}}$	K	${9480}_{-410}^{+420}$	${4890}_{-170}^{+190}$	3460 ± 60
$\mathrm{log}g$	[dex]	4.32 ± 0.01	4.40 ± 0.03	4.45 ± 0.05

Download table as:  ASCIITypeset image

We find a mass of ${2.08}_{-0.11}^{+0.12}$ M_⊙, a temperature of ${9480}_{-410}^{+420}$ K, and a surface gravity of $\mathrm{log}g=4.32\pm 0.01$ [dex] for HD 131399 A. These parameters are consistent with an A1V spectral type (Houk 1982) at an age of 16 Myr. The extinction toward HD 131399 of ${A}_{V}=0.22\pm 0.09$ mag estimated from the SED fit is consistent with literature estimates that range from 0.14 to 0.28 mag (de Geus et al. 1989; Sartori et al. 2003; Chen et al. 2012). The photometric distance of ${107.9}_{-3.7}^{+4.5}$ pc is 1.2σ discrepant from the trigonometric distance of ${98.0}_{-6.3}^{+7.2}$ pc from the Hipparcos parallax (van Leeuwen 2007). Repeating the SED fit using only a $p(\pi )\propto {\pi }^{-4}$ prior, corresponding to an assumed uniform space density of stars, results in a similar photometric distance of ${112.2}_{-5.1}^{+5.2}$ pc. The stated uncertainties on the fitted parameters do not incorporate any model uncertainty and are therefore likely underestimated.

The SED of each component and that of the blended system are shown in Figure 7. Uncertainties on the near-IR portion of the SED of HD 131399 A, estimated by sampling randomly from the posterior distributions (t, π, ${M}_{{\rm{A}}}$, and A_V), ranged between 1.5% and 2.0%. The SED of A was degraded to the spectral resolving power of the GPI and SPHERE IFS observations to convert the contrasts between HD 1313199 A and Ab measured in Section 2 into apparent fluxes for Ab.

Zoom In Zoom Out Reset image size

Synthetic photometry of HD 131399 A was also computed for the GPI, SPHERE, and NIRC2 filters to convert the measured broadband contrasts between A and Ab into apparent magnitudes for Ab. Filter transmission profiles for the GPI filters were obtained from the GPI DRP and were combined with a median Cerro Pachón atmosphere (4.3 mm precipitable water vapor) at one air mass (Lord 1992). The SPHERE IRDIS filter curves were obtained from the ESO website,⁴² while the IFS throughput was assumed to be uniform between 0.96 and 1.11 μm at Y, 1.13–1.42 μm at J, and 1.44–1.64 μm at H. These filter curves were combined with a median Paranal atmosphere (2.5 mm precipitable water vapor) at one air mass (Moehler et al. 2014). The NIRC2 ${L}^{{\prime} }$ filter curve was obtained from the Keck website⁴³ and was combined with a Maunakea atmosphere (Lord 1992) with 1.6 mm of precipitable water vapor at two air masses (chosen to match the observing conditions on 2017 February 08). The throughput of the GPI and SPHERE filters are plotted in Figure 8. Zero points and effective wavelengths for all of the filters were estimated using the CALSPEC Vega spectrum⁴⁴ (Bohlin 2014) and are given in Table 4. The properties of the ${H}_{\mathrm{GPI}}$ and $K{1}_{\mathrm{GPI}}$ filters match those derived from observations of the white dwarf HD 8049 B presented in De Rosa et al. (2016).

Zoom In Zoom Out Reset image size

Table 4.  Atmosphere Throughput-corrected Filter Properties

Filter	${\lambda }_{\mathrm{eff}}$	${W}_{\mathrm{eff}}$	Zero point
	(μm)	(μm)	(10−9 Wm−2μm−1)
${J}_{\mathrm{GPI}}$	1.23	0.19	3.12
${H}_{\mathrm{GPI}}$	1.64	0.27	1.15
$K{1}_{\mathrm{GPI}}$	2.06	0.20	0.50
${Y}_{\mathrm{SPH}-\mathrm{IFS}}$	1.03	0.16	5.65
${J}_{\mathrm{SPH}-\mathrm{IFS}}$	1.24	0.24	2.99
${H}_{\mathrm{SPH}-\mathrm{IFS}}$	1.54	0.19	1.41
JSPH	1.23	0.20	3.11
$K{1}_{\mathrm{SPH}}$	2.10	0.09	0.47
$K{2}_{\mathrm{SPH}}$	2.25	0.11	0.36
${L}^{{\prime} }$	3.72	0.59	0.054

Note. The subscript SPH refers to the SPHERE IRDIS filters and SPH-IFS to the derived SPHERE IFS filters, to differentiate between them.

Download table as:  ASCIITypeset image

Photometric measurements, S/N, and spectra obtained from GPI, from NIRC2, from the new SPHERE data, from our reanalysis of the SPHERE data, and from those published by W16 are given in Table 5 and Figure 9. The full SEDs for both HD 131399 A and Ab, from both the SPHERE and GPI data sets, is given in Tables 8 and 9 in the Appendix. The measurements provide YJH contrasts consistent at the 1σ level between the four SPHERE sets, between our average SPHERE contrasts and those published in W16, and between our average SPHERE and average GPI measurements, with the caveat that the GPI and SPHERE filters are different (especially H, see Figure 8). However, the reanalyzed SPHERE contrast at K1 and K2 differs significantly (2σ at K1 and 1σ at K2) from that of W16. The origin of these discrepancies remains unclear since the K1 contrasts of HD 131399 B and C are in agreement between our reanalysis (${\rm{\Delta }}K1=1.95\pm 0.07$ and 3.84 ± 0.10 mag) and that of W16 (${\rm{\Delta }}K1=1.86\pm 0.10$ and 3.86 ± 0.10 mag for B and C, respectively).

Zoom In Zoom Out Reset image size

The GPI spectrum is flat, except for some correlated noise, at a high confidence level, without any indication of the methane absorption beyond 1.6 μm that is expected in the spectra of mid-T dwarfs. The GPI spectrum is also in agreement with that of the combined four SPHERE sets in both the J and H bands. However, the published SPHERE H-band spectrum (W16) peaks at 1.61 μm, a peak that does not appear either in the GPI spectrum or in our reanalysis. The peak flux is also nearly twice the plateau of the other two spectra. These differences might be explained by (1) a different technique used to combine the multiple data sets or (2) the technique used to extract the photometry of Ab, with different techniques being biased by nearby speckles to varying degrees. The presence of a speckle close to HD 131399 Ab in the 2016 March 17 SPHERE data set may be significantly biasing the spectrum at ∼1.6 μm (see Section 2.2), with the spectrum being featureless in the three other sets.

3.2.1. Color–Magnitude and Color–Color Diagrams

The physical nature of HD 131399 Ab can be assessed by placing it on a color–magnitude or color–color diagram (CMD or CCD) and comparing it to the location of other objects of known spectral types. A library of medium-resolution ($R\sim 200$) near-IR spectra of stars and brown dwarfs was compiled from the SpeX Prism library⁴⁵ (Burgasser 2014), the IRTF Spectral Library⁴⁶ (Cushing et al. 2005), and the Montreal Spectral Library⁴⁷ (e.g., Gagné et al. 2015; Robert et al. 2016). The spectra were normalized to literature 2MASS or MKO photometry. Parallax measurements were obtained from Dupuy & Liu (2012), Dupuy & Kraus (2013), Liu et al. (2016) (and references therein) for the brown dwarfs, and from van Leeuwen (2007) for the stars. Synthetic magnitudes in the GPI and SPHERE filters were calculated for each object using the filter curves shown in Figure 8. We generated an M_J versus ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ CMD, and $K{1}_{\mathrm{SPH}}-K{2}_{\mathrm{SPH}}$ versus ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ and $K{1}_{\mathrm{GPI}}-K{1}_{\mathrm{SPH}}$ versus ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ CCDs, all of which are plotted in Figure 10. A $K{1}_{\mathrm{GPI}}-{L}^{{\prime} }$ versus ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ CCD was also created, shown in Figure 11 using literature MKO ${L}^{{\prime} }$ photometry, or estimated from the WISE W1 to MKO ${L}^{{\prime} }$ color transformation given in De Rosa et al. (2016). No extinction correction was applied to the colors, although this is expected to be small (${A}_{V}\sim 0.22$, ${A}_{J}\sim 0.06$, ${A}_{H}\sim 0.04$ mag) at the distance to HD 131399 A, increasing to ${A}_{V}\sim 1$ mag (${A}_{J}\sim 0.29$, ${A}_{H}\sim 0.18$ mag) due to a combination of the extinction within the UCL region and the predicted extinction from galactic dust. The locations of field-gravity standards for spectral types later than M0 are highlighted in each diagram (Burgasser et al. 2006b; Kirkpatrick et al. 2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the contrasts between A and Ab reported in Table 5 and the synthetic magnitudes for A calculated in Section 3.1, we derive colors of ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}=0.54\pm 0.18$ mag, $K{1}_{\mathrm{SPH}}\,-K{2}_{\mathrm{SPH}}=0.22\pm 0.14$ mag, and $K{1}_{\mathrm{GPI}}-K{1}_{\mathrm{SPH}}=-0.09\,\pm 0.18$ mag for HD 131399 Ab. We also derive an upper limit of $K{1}_{\mathrm{GPI}}-{L}^{{\prime} }\gt 1.52$ mag, using the detection limit from the NIRC2 ${L}^{{\prime} }$ observations. On each of the CCDs in Figure 10, HD 131399 Ab is consistent with the colors of M dwarfs and is significantly different from the observed colors of early to mid-T dwarfs, a discrepancy that is most significant for the measured $K{1}_{\mathrm{SPH}}-K{2}_{\mathrm{SPH}}$ color. As a comparison, the upper limit on the color of 51 Eri b of $K{1}_{\mathrm{SPH}}-K{2}_{\mathrm{SPH}}\lt -0.58\,\pm 0.14$ mag (Samland et al. 2017) is more than 3σ discrepant. The position of HD 131399 Ab on the $K{1}_{\mathrm{GPI}}-{L}^{{\prime} }$ versus ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ CCD (Figure 11) only excludes mid to late Ls and late Ts; M dwarfs and mid-Ts are consistent with the measured ${J}_{\mathrm{GPI}}-{H}_{\mathrm{GPI}}$ color and the $K{1}_{\mathrm{GPI}}-{L}^{{\prime} }$ upper limit.

While the absolute M_J magnitude is consistent with an early to mid-T dwarf, the J − H color is far less diagnostic (Figure 10, top panel). If the distance to HD 131399 Ab was not known, the only constraint on the spectral type from the J − H color would be that it is between mid-G and late-M, or between early and mid-T. Excluding the J − H color, the only evidence in support of the bound T-dwarf companion hypothesis from these color–magnitude and color–color diagrams is the absolute J-band magnitude, which relies on the assumption that it is at the same distance as HD 131399 A (${98.0}_{-6.3}^{+7.2}$ pc), and the upper limit on the $K{1}_{\mathrm{GPI}}-{L}^{{\prime} }$ color, which is consistent with either an M dwarf or a mid-T dwarf. The remaining color indices plotted in Figure 10 except for J − H are inconsistent with the observed colors of field T dwarfs. Instead, they are consistent with those of field M dwarfs, which would require HD 131399 Ab to be at a significantly greater distance of between 1 and 10 kpc and not physically associated with HD 131399 A.

3.2.2. Comparison to Spectra of Field Objects

One of the primary reasons why instruments such as GPI and SPHERE use an integral field spectrograph is the ability to immediately distinguish between background stars, which have relatively featureless spectra, and cool substellar companions with strong molecular absorption features. While the J − H of HD 131399 Ab is consistent with both stars between mid-G and late-M and brown dwarfs between early-T and mid-T (Figure 10), the JH spectra of these two groups of objects are significantly different. With a high enough S/N spectrum, it should be possible to confirm or reject the presence of strong molecular absorption features that are seen in the spectra of cool brown dwarfs.

We compared the GPI and our SPHERE spectra of HD 131399 Ab to the library of near-IR spectra described in Section 3.2.1. The spectrum of each object within the library was degraded to the resolution of the GPI/SPHERE spectra by convolving the spectrum with a Gaussian of appropriate width. The scaling factor that minimized ${\chi }^{2}$ was found analytically for the comparison to the SPHERE data and numerically for the comparison to the GPI data, where the separate bands were allowed to float independently to account for uncertainties in the satellite spot ratio (Maire et al. 2014). Wavelengths with a throughput lower than 50% (Figure 8) were excluded from the fit.

The fits of the HD 131399 Ab SPHERE and GPI spectra to objects ranging from a spectral type of G0 to T6 are shown in Figure 12, with the minimum ${\chi }_{\nu }^{2}$ plotted as a function of spectral type in Figure 13. The lower S/N of the SPHERE spectrum is apparent (Figure 12, left panel), with ${\chi }_{\nu }^{2}\lt 1$ for all spectral types except for those between L5–L9 and T5–T9 (Figure 13). The SPHERE spectrum is fit well (${\chi }_{\nu }^{2}\lt 1$) by objects that have significantly different spectral morphologies: from an M5 dwarf (Wolf 47, ${\chi }_{\nu }^{2}=0.50$), with a relatively featureless spectrum, to a T2 (2MASS J12545393–0122474, ${\chi }_{\nu }^{2}=0.63$) or a T4 brown dwarf (2MASSI J2254188+312349, ${\chi }_{\nu }^{2}=0.66$), which exhibit strong molecular absorption features. The YJ portion of the spectrum is consistent within the uncertainties with spectral types earlier than T6, providing little diagnostic power. The H-band spectrum exhibits a rising slope toward longer wavelengths, similar to what is seen in the spectra of brown dwarfs later than L5, although this slope is not measured at a significant level given the low S/N.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The improved S/N and greater wavelength coverage of the GPI JHK1 spectrum provide for better constraints on the spectral type of HD 131399 Ab (Figure 12). The spectrum appears relatively featureless, consistent with the near-IR SED of stars with a spectral type earlier than mid-M. The red end of the H spectrum appears to modulate on a characteristic length scale consistent with the intrinsic resolution of GPI at H. It is likely this is correlated noise due to the presence of speckles at those wavelengths rather than an astrophysical signal. The GPI spectrum is fit well by both an M0 (BD+33 1505, ${\chi }_{\nu }^{2}=0.67$) and an M5 (Wolf 47, ${\chi }_{\nu }^{2}=0.67$) dwarf. Earlier spectral types are also fit well (${\chi }_{\nu }^{2}\lt 1$), although these would require HD 131399 Ab to be at a significantly greater distance, inconsistent with the predictions of Galactic population models described in Section 4.4. The minimum ${\chi }_{\nu }^{2}$ for the fit of the GPI spectrum as a function of spectral type plotted in Figure 13 displays a trend similar to that for the fit of the SPHERE data, with later spectral types being more strongly excluded. Objects earlier than a spectral type of L0 fit the spectrum relatively well (${\chi }_{\nu }^{2}\lt 1$). One limitation of this analysis is the relative dearth of known young/low-surface-gravity T dwarfs. The three within the library—HN Peg B (T2.5, ${\chi }_{\nu }^{2}=3.0$, Luhman et al. 2007), 2MASS J11101001+0116130 (T5.5, ${\chi }_{\nu }^{2}=8.9$, Burgasser et al. 2006a), and CFBDSIR J214947.2-040308.9 (T7, ${\chi }_{\nu }^{2}=16.0$, Delorme et al. 2013)—are all poor fits to the GPI spectrum of HD 1313199 Ab.

Using the color–magnitude and color–color diagrams in Figure 10 and the fit of the SPHERE and GPI spectra to stars and brown dwarfs in Figures 12 and 13, we find no strong evidence to suggest that HD 131399 Ab has a near-IR SED consistent with that of a cool planetary-mass companion of early to mid-T spectral type. We do not detect the characteristic H₂O and CH₄ absorption in the GPI spectrum at either J or, more significantly, at H, nor do we detect it based on the measured $K{1}_{\mathrm{SPH}}-K{2}_{\mathrm{SPH}}$ color, which is sensitive to methane absorption in the spectra of T dwarfs (Figure 10, middle panel). Instead, our analysis of the near-IR SED suggests it has a relatively featureless spectrum and has near-IR colors that are consistent with those of a low-mass star.

Before fitting an orbit, we first consider whether the projected velocity of HD 131399 Ab is less than the escape velocity of the system, as should be true for a bound orbit. The projected velocity (in R.A. and decl.) will in fact be a lower limit on the total velocity, since the total velocity will also include the unmeasured component along the line of sight. We compute projected velocity by fitting straight lines to the astrometry in R.A. and decl. as a function of time. This too represents a lower limit on the velocity, since any curvature not captured by the linear fit would represent a higher velocity. This value is then converted to a physical velocity (km s⁻¹) using the distance to the system.

Escape velocity is given by ${v}_{\mathrm{esc}}=\sqrt{2{GM}/r}$, where G is the gravitational constant, M is the mass of the star, and r is the total separation between star and planet. In the direct imaging case, this corresponds to an upper limit on escape velocity, since we can only measure the projected separation in R.A. and decl. In fact, the presence of the binary BC would lower the effective escape velocity further beyond this upper limit, since the planet would not need sufficient velocity on its own to reach infinity, but only enough velocity to reach the gravitational sphere of influence of BC to eventually escape. Separation is computed using the minimum value over the range of epochs of the astrometry (2015 June 12 through 2017 April 20) from the linear fit, with the minimum value chosen so we continue to define the upper limit of the escape velocity.

In order to compare the projected velocity to the escape velocity limit, we use a Monte Carlo method to draw samples from both velocities given uncertainties in the astrometry, distance to the system, and mass of the star. For each Monte Carlo trial, for both R.A. and decl., we generate values of slope (projected velocity) and intercept (reference position) from the covariance matrix of the linear fit to the astrometry, as well as stellar mass and parallax from Gaussian distributions (2.08 ± 0.11 ${M}_{\odot }$ and 10.20 ± 0.70 mas, respectively). Finally, we compute the ratio of projected velocity to escape velocity, which should be less than unity for a bound orbit. This method accounts for correlations in distance, since the same generated distance is used to calculate projected velocity and escape velocity, as well as between slope and intercept, since the same generated pair is used to compute the minimum separation for the escape velocity as well as the projected velocity.

In Figure 14 we plot the total projected displacement (total distance in both R.A. and decl.) between the first epoch and all subsequent epochs. The reference location is taken as the average of the IRDIS and IFS astrometry at 2015 June 12. We draw 100 values of the escape velocity from our Monte Carlo analysis and plot these as red lines, normalized to pass through the reference location. The astrometry clearly shows a steeper slope (a faster velocity) than the escape velocity. Computing the quotient of projected velocity and escape velocity (bottom panel of Figure 14) shows that the projected velocity is indeed always greater than the escape velocity, with the ratio having a value of 1.89 ± 0.23, which reached a minimum value of 1.07 out of 10⁷ trials. Thus the data are robustly inconsistent with the hypothesis that HD 131399 Ab is a bound planet.

Zoom In Zoom Out Reset image size

The quantity $v/{v}_{\mathrm{esc}}$ is inversely proportional to the square root of stellar mass and directly proportional to ${d}^{1.5}$, with a factor of d coming from the projected velocity and another factor of $\sqrt{d}$ from the escape velocity. The 2σ lower limit on this quantity is 1.56 (95.45% of the samples are larger than this number). In order to bring this 2σ limit to unity, it is therefore necessary to increase the mass of HD 131399 A by a factor of ${1.56}^{2}=2.43$ or decrease the distance by ${1.56}^{1/1.5}=1.35$ (or else have a linear combination of these two changes). Such a change would represent a 27σ deviation in mass or a 5.1σ deviation in distance. Such a change in distance would likely exclude HD 131399 A from the UCL association, and therefore the star and HD 131399 Ab would be much older, which ultimately affects the model-dependent mass estimate of the latter. Of these two, the most susceptible to error is mass, since the mass of the primary comes from an SED fit, and lower-mass stellar companions, too close to be resolved with GPI, could add additional mass that would raise the escape velocity. However, it is difficult to imagine there being 3 ${M}_{\odot }$ of additional stars close to the 2.08 ${M}_{\odot }$ HD 131399 A. The most likely high-mass companion would be an equal-mass binary, which even then is not enough to make the orbital velocity equal escape velocity at the 2σ level, and would be evident in the distance posterior from the SED fit.

This large upper limit on $v/{v}_{\mathrm{esc}}$ is not solely dependent on the astrometric calibration between SPHERE and GPI. When we repeat the same analysis for only the SPHERE astrometry presented here, five epochs from 2015 to 2017, we find this factor has a value of 1.88 ± 0.25, with only four out of 10⁷ trials less than unity. Using the original astrometry reported by W16, this factor becomes consistent with bound orbits, 1.02 ± 0.42, with 48% of generated values below unity. In contrast, when we use our astrometry for these same four epochs from 2015 to 2016, we find a value of 1.80 ± 0.31, with 0.09% of trials less than unity, consistent with our finding of a significantly larger slope in our reduction of the 2015–2016 SPHERE data compared to W16.

In the previous section, we have demonstrated that the projected velocity of HD 131399 Ab is significantly above the escape velocity. We proceed to explore the magnitude of the offset required in the astrometry, mass, and distance in order to fit a bound orbit to the data. We begin by assuming a fixed mass and distance of the star, 2.08 M_⊙ and 98.0 pc. In order to investigate the orbital parameters required to fit the astrometry, we use the rejection sampling algorithm OFTI (De Rosa et al. 2015; Blunt et al. 2017; Rameau et al. 2016). Using OFTI, we generate 100 orbits drawn from the posterior probability distribution and plot them in Figure 15.

Zoom In Zoom Out Reset image size

Unsurprisingly, since the projected motion is faster than the escape velocity, the best-fitting orbit is a poor fit to the data, with a systematically steeper slope in the data than the fit. With ${\chi }_{\nu }^{2}=6.7$, even the best-fitting orbit is clearly a bad fit to the data. This high projected velocity is only possible with very high orbital eccentricity, $e\gt 0.949$ for all generated orbits. This results in a high value of apastron, with 68% confidence between 1017 and 22433 au, and a minimum value of 597 au. The projected separation of the closest of the BC pair, HD 131399 B, with respect to HD 131399 A, is 309 au. Thus a large semimajor axis for BC around A ($\gtrsim 1000$ au) and a highly inclined orbit ($\approx 70\mbox{--}110^\circ$, so that the projected separation is only ∼300 au) would be required for these orbits to not cross each other.

When fitting orbits, it is more correct to incorporate errors on mass and distance by varying these parameters in the orbit fit and imposing priors as Gaussians given the measurements. We noted in Section 4.1 that the escape velocity problem can be ameliorated by increasing the mass of the star or decreasing the distance, so this standard orbit method will have the result of balancing, in a Bayesian sense, the ${\chi }^{2}$ of the orbit fit, the mass, and parallax to find the most likely compromise between the three.

To investigate the effect of allowing the distance and mass of the star to vary within the orbit fit, we use the emcee parallel-tempered affine-invariant MCMC sampler to estimate the orbital elements from the astrometry presented in Table 6. We fit eight parameters: semimajor axis a, inclination i, eccentricity e, position angle of nodes and argument of periastron as ${\rm{\Omega }}+\omega$ and ${\rm{\Omega }}-\omega$, epoch of periastron T₀, parallax π, and M (total mass, the mass of Ab being negligible if bound). Here we define T₀ in units of orbital period from the first epoch of the astrometric record (2015.44). We adopt uniform priors in ${\mathrm{log}}_{10}a$, $\cos i$ (−1 to 1), e ($\lt 1$), ${\rm{\Omega }}+\omega$ and ${\rm{\Omega }}-\omega$ (0–2π), and T₀ (0–1). The prior on π was created by multiplying a Gaussian distribution centered at 10.20 mas with a 1σ width of 0.70 mas, corresponding to the Hipparcos parallax of HD 131399 A, with a ${\pi }^{-4}$ power-law distribution. The prior on the mass, M, was a Gaussian distribution centered at 2.08 M_⊙ with a 1σ width of 0.11 M_⊙ (Table 3). We initialized 512 walkers at each of 32 temperatures. Each walker was advanced for 10⁶ steps, with the first half of each chain discarded as a "burn in" as they converged to their final value.

Allowing the distance and mass to float significantly improved the quality of the fit, reducing the minimum ${\chi }_{\nu }^{2}$ from 6.7 (Figure 15) to 0.92 (Figure 16). This improvement was achieved by the MCMC chains moving to a significantly smaller distance to HD 131399 A (∼73 pc) and a slightly larger total mass (∼2.25 M_⊙). This decreased the measured velocity of HD 131399 Ab and increased the escape velocity of the system so that bound orbits could be fit. The posterior distribution of the parallax ($\pi =13.79\pm 0.46$ mas) is 4.3σ discrepant from the prior distribution (Figure 16) and corresponds to a distance of 72.5 ± 2.4 pc, consistent with the distance required in the escape velocity analysis in Section 4.1. This distance is significantly discrepant from the Hipparcos measurement of 98.0 ± 6.9 pc, the distance obtained from the SED fit of the three stars in the HD 131399 system of 107.9 ± 4.0 pc (Section 3.1), and also the mean distance of UCL members of 140 pc (de Zeeuw et al. 1999). The posterior distribution of the mass ($M=2.27\pm 0.10$ M_⊙) is shifted by 1.3σ relative to the prior distribution ($M=2.08\pm 0.11$ M_⊙). As with the fit using OFTI with a fixed mass and distance, the posterior distribution of e is strongly peaked at very high eccentricities. We find a median and 1σ range of $e={0.980}_{-0.017}^{+0.010}$, and the lowest eccentricity within any of the MCMC chains was e = 0.82.

Zoom In Zoom Out Reset image size

Most stars targeted by direct imaging are relatively nearby (≲100 pc), so they typically have well-measured parallaxes and proper motions (errors ≲1 mas). Thus the motion of the target star across the sky is well determined over time. The analysis then proceeds by comparing the relative astrometry between candidate and target star over multiple epochs and determining whether it follows the background track (which, relative to the target star, moves in the opposite direction of the parallax and proper motion of the star), or is more consistent with common proper motion (e.g. Nielsen et al. 2013).

An added complication is that candidates closer to the star may show significant orbital motion over the time frame of the astrometric observations, as detailed in the case of the planets HR 8799 bcd (Marois et al. 2008), β Pic b (Lagrange et al. 2010), and 51 Eri b (De Rosa et al. 2015). A determination of the status of the candidate can be made when it is shown to follow the background track at a projected velocity (assuming the distance of the target star) inconsistent with a bound orbit, or it does not follow the background track and has motion consistent with a bound orbit.

Figure 17 presents this analysis for HD 131399 Ab and demonstrates that neither the infinitely far background object scenario nor the orbiting planet scenario is fully consistent with all the astrometric data. The background track is tied to the GPI 2017 February 15 H-band point and assumes the Hipparcos proper motion for HD 131399 A of (−26.69 ± 0.59, −31.52 ±0.55) mas yr⁻¹ and parallax of 10.2 ± 0.7 mas (van Leeuwen 2007). The width of the gray track corresponds to the 68% confidence interval, based on a Monte Carlo error analysis that combines the errors on the reference epoch astrometry and the star's proper motion and parallax (Nielsen et al. 2013). The orbital motion cones are calculated using the same 2017 February 15 GPI astrometry, using the Monte Carlo method described in De Rosa et al. (2015). Orbital parameters are drawn from prior distributions: uniform priors in ω and T₀, i following a $\sin i$ distribution, and e following a linear fit to RV planets (Nielsen & Close 2010). Here, a and Ω are chosen to place the reference astrometry on the orbit at the reference epoch, to within a two-dimensional Gaussian centered on the measured astrometry and with standard deviation equal to the measurement errors.

Zoom In Zoom Out Reset image size

The astrometry falls neither on the background track nor within the orbit cone. The 2015 June 12 SPHERE data, in particular, are clearly between the separation values predicted for a zero-proper-motion background object and an orbiting planet. We have demonstrated above that orbital motion cannot explain this offset, so we instead reexamine the assumption that a background object has no proper motion or parallax of its own.

As noted above, the standard assumption when testing for common proper motion is that the background object is at infinite distance with zero proper motion. In reality, the typical distances to background stars are ∼1–10 kpc, with proper motions less than a few mas per year. We investigate this possibility using the MCMC method described in Macintosh et al. (2014) and De Rosa et al. (2015) to find the proper motion and parallax a background object would need to match the relative astrometry of a candidate companion. In these previous works, we assumed no errors on the astrometry at the reference epoch and no errors on the proper motion and parallax of the primary. Here we update this method by incorporating these errors and fitting an eight-dimensional function: proper motion (in R.A. and decl.) and parallax of the candidate, proper motion and parallax of the primary, and separation and PA of the candidate with respect to the primary star at a reference epoch (chosen to be 2017.0 so as to be in the middle of our astrometric record). Priors are taken to be uniform in proper motion in R.A. and decl. of the candidate. We adopt the distance prior described by Bailer-Jones (2015), which combines a uniform space density prior with an exponential drop-off in stellar density, $p(d)\propto {d}^{2}{e}^{-d/L}$, where d is distance and L is a reference length scale, which is set to 1000 pc. Changing variables from distance to parallax ($\pi =1/d$) introduces an additional factor of $1/{\pi }^{2}$, giving us our parallax prior of $p(\pi )\propto {\pi }^{-4}{e}^{-1/(\pi L)}$, which is truncated at 10 kpc. For the primary, Gaussian priors are used for the proper motion and parallax corresponding to the van Leeuwen (2007) Hipparcos measurements and errors. Uniform priors are assumed for the separation and position angle at the reference epoch.

In Figure 18 we present the posteriors on parallax and proper motion of HD 131399 Ab, assuming it is not bound to HD 131399 A. While proper motion in the R.A. direction and parallax are both close to 0 (${\mu }_{\alpha }=-4.7\pm 1.6$ mas yr⁻¹, $\pi \,\lt$ [0.64, 2.01] mas at [68%, 95%] confidence), decl. proper motion is significantly larger (${\mu }_{\delta }=-11.2\pm 1.3$ mas yr⁻¹). This motion in decl. is the departure from the background track seen in Figure 17 and the orbital motion examined by W16 and discussed in Section 4.2. We display the fit proper motions along with the data in Figure 19, where the new proper motion track is the difference between the Hipparcos values of proper motion and parallax of HD 131399 A and our fit values for the proper motion and parallax of HD 131399 Ab. Errors from the Hipparcos measurement and uncertainties in our fit are added in quadrature. While ∼12 mas yr⁻¹ is a relatively large proper motion, about one-quarter the total proper motion of HD 131399 A, is this motion plausible for a star at ∼1 kpc?

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To answer this question, we use the Besançon model of stellar populations (Robin et al. 2003), retrieving a set of simulated stars from the web form at http://model.obs-besancon.fr. We selected stars in the direction of HD 131399 A, with magnitudes of $19.63\lt H\lt 19.87$ to match the 2σ range of apparent magnitudes of HD 131399 Ab, with distances from 0 to 50 kpc, and with a large solid angle of one square degree to give us a large statistical sample (6197 stars were generated). We plot a subset of 1000 stars, along with the constraints on the proper motion and parallax of HD 131399 Ab from the relative astrometry, in Figure 20. Since our parallax posterior is largely set by our choice of prior, and our data cannot distinguish between parallaxes $\lesssim 5$ mas, we extend the contours from 0.1 mas to zero. While the proper motion and parallax constraints do not encompass the majority of the simulated background stars, there is a significant subset that fall within the contours: 0.16% fall within the 1σ contours, 0.89% inside 2σ, and 1.9% within 3σ. We note that these values do not represent the probability that HD 131399 Ab is a background object, but are instead proportional to our constraints on proper motion and parallax. Even if our constraints were in the middle of the cloud of background objects in Figure 20, as better astrometry allowed us to reduce the size of the contours, fewer and fewer simulated background stars would fall inside. Rather, this is a demonstration that the proper motion required to explain the change in relative astrometry seen in Figure 17 is plausible for background stars in the direction of HD 131399 A with the same apparent magnitude as HD 131399 Ab.

Zoom In Zoom Out Reset image size

The Besançon points are consistent not just with the amplitude of the proper motion, but also with the direction we measure for HD 131399 Ab. These points, in the left panel of Figure 20, are not distributed isotropically, but instead preferentially represent stars moving south and west. If we were simply to reverse the direction of HD 131399 Ab, multiplying the measurements of proper motion in R.A. and decl. by −1, the fraction of Besançon stars falling into the [1, 2, 3] σ contours drops to [0%, 0%, 0.048%]. While there is no preferred direction for an orbiting planet, there is clearly a preferred direction for the proper motion of background stars at these coordinates, and HD 131399 Ab is moving in that direction.

We note the constraints at smaller distances (≲1 kpc) are largely driven by our choice of the prior on parallax, but do not strongly affect the overall conclusions. Constraints on proper motion in R.A. and decl. are largely unchanged by the choice of prior. Changing from a $p(d)\propto {d}^{2}{e}^{-d/L}$ prior to a simple $p(d)\propto {d}^{2}$ cutoff at 10 kpc results in the fraction of Besançon stars falling into [1, 2, 3] σ contours dropping from [0.16%, 0.89%, 1.9%] to [0.03%, 0.29%, 1.3%], entirely due to lower values of parallax becoming less probable. Switching to a uniform prior in parallax results in again smaller percentages falling into the contours compared to the exponentially declining prior, [0%, 0.4%, 1.5%], as now more large-parallax points are accepted while more small-parallax points are rejected. Of these three priors, both the uniform and d² prior are poor fits to the Besançon points, while the exponentially declining prior is an excellent fit between 1 and 10 kpc. Since our measurement on parallax is essentially an upper limit, we adopt the values from the more realistic, exponentially declining prior in our analysis, but with the small parallax constraints extended from 0.1 mas to 0 mas.

We are presented with two unlikely scenarios for the nature of HD 131399 Ab: a planet with extreme orbital parameters (most likely currently being ejected from the system), or a background object with an unusually high proper motion. An order-of-magnitude estimate for the likelihood of observing a planet being ejected just as it is detected is the orbital period divided by the system lifetime, which for a circular orbit with semimajor axis equal to the projected separation of 82 au, is 520 yr / 16 Myr, or $9\times {10}^{-6}$. So while a background object whose proper motion is only consistent with ∼1% of Besançon simulated objects (at 95% confidence) is unlikely, the ejected planet hypothesis is several orders of magnitude more unlikely.

4.5.1. Bound Planet or Background Star?

To apply a more rigorous analysis, we construct an odds ratio between the likelihood of the planet and background object scenario. In particular, we consider three elements for each hypothesis: the overall likelihood of each object, the relative likelihood as a function of separation from the star, and relative probability as a function of projected velocity:

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{BG}}}{{P}_{\mathrm{Pl}}}=\displaystyle \frac{P(\mathrm{BG})P(\rho | \mathrm{BG})P(v| {\rho }_{\mathrm{BG}})}{P(\mathrm{Pl})P(\rho | \mathrm{Pl})P(v| {\rho }_{\mathrm{Pl}})},\end{eqnarray} \tag{ 1 }$$

where P is probability, BG and Pl are background and planet, ρ is projected separation, and v is projected velocity.

For the first term, we restrict the separation range to $0\buildrel{\prime\prime}\over{.} 1$ to $1^{\prime\prime}$ (≈10–100 au) and apparent magnitude $19.63\lt H\lt 19.87$. For background objects, this is the total number of Besançon objects generated in Section 4.4, 6197 stars, multiplied by the ratio of solid angles: $6197\times {(\pi {(1^{\prime\prime} )}^{2}-\pi {(0\buildrel{\prime\prime}\over{.} 1)}^{2})/((3600^{\prime\prime} )}^{2})\,=0.15 \%$. For planets, an apparent magnitude constraint of $19.63\lt H\lt 19.87$, at 16 Myr, corresponds to a mass range of 3.64–4.01 ${M}_{\mathrm{Jup}}$, using the COND models (Baraffe et al. 2003). Early analysis of the Gemini Planet Imager Exoplanet Survey puts planet yield at ∼6% of 2 ${M}_{\odot }$ hosting a planet between 10 and 100 au, with mass between 5 and 13 ${M}_{\mathrm{Jup}}$ (E. Nielsen et al. 2018, in preparation). We convert this probability to our mass range, 3.64–4.01 ${M}_{\mathrm{Jup}}$, assuming planets follow the power-law distribution in mass of Cumming et al. (2008), ${dN}/{dM}\propto {M}^{-1.31}$, or 0.76% of stars having such a planet.

Next, we consider the second term, the fraction of background stars and planets at the observed projected separation, taken here to be the 2σ range measured at the discovery epoch with the SPHERE IFS, 836.6–848.2 mas. Over small areas on the sky ($\lesssim 10^{\prime\prime}$), background objects are uniformly distributed, so the probability of lying at this separation is just the ratio of areas between an annulus from 836.6 to 848.2 mas and one from 100 to 1000 mas, or 1.97%. For planets, the semimajor axis distribution (converted from the Cumming et al. 2008 distribution for period for solar-mass stars) is ${dN}/{da}\propto {a}^{-0.61}$. So the second term for planets is then 0.85%. Adopting a uniform log semimajor axis distribution (${dN}/{da}\propto {a}^{-1}$) instead would have a minor effect, changing this term from 0.85% to 0.60%.

Finally, we consider the projected velocity measured between 2015 and 2017. As described in Section 4.4, 0.89% of Besançon objects have a proper motion consistent with our 2σ contours, so we use this as the value for $P(v| {\rho }_{\mathrm{BG}})$. The corresponding term for planets, $P(v| {\rho }_{\mathrm{Pl}})$, is then the fraction of orbits whose projected velocity is within the 2σ range of the measured projected velocity from the relative astrometry. Using OFTI (Blunt et al. 2017), we generate 10⁷ orbits from the first SPHERE IFS astrometric measurement, incorporating errors in separation, position angle, stellar mass, and distance. The maximum projected velocity expected from an orbiting body is then 19.7 mas yr⁻¹. In Section 4.1 we generated 10⁷ measurements of the projected velocity from the full astrometric record of HD 131399 Ab, finding 2σ confidence intervals in the ranges 15.3–21.5 mas yr⁻¹ in R.A. and 18.5–22.6 mas yr⁻¹ in decl. As expected, since we found this projected velocity to be above escape velocity, none of the 10⁷ generated orbits has a projected velocity in this 2σ range, so we can only set an upper limit on $P(v| {\rho }_{\mathrm{Pl}})$ of 10⁻⁷, or 10⁻⁵%.

Entering these values into Equation (1), we find

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{BG}}}{{P}_{\mathrm{Pl}}}\gt \displaystyle \frac{0.15 \% \times 1.97 \% \times 0.94 \% }{0.76 \% \times 0.85 \% \times {10}^{-5} \% }=43,000,\end{eqnarray} \tag{ 2 }$$

so the probability that HD 131399 Ab is a background object is 43,000 times greater than the probability that it is a bound planet. While this analysis depends on multiple assumptions, including the poorly known distribution of giant planets at wide separations, the dominant term in Equation (2) is the velocity term. If we were to consider the stability of these orbits, the likelihood of being a bound planet would drop further, since the best-fitting orbits are likely to cross the orbit of HD 131399 BC. But stability calculations (e.g., Holman & Wiegert 1999) show that the semimajor axis of the planet must actually be at least a factor of two smaller than that of HD 131399 BC, so even noncrossing orbits are not sufficient. As we would expect from our analysis of the escape velocity, the large projected velocity of HD 131399 Ab strongly precludes the possibility that it is a bound planet.

4.5.2. Ejected Planet?

4.5.3. Alternative Explanations