RECONNAISSANCE OF THE HR 8799 EXOSOLAR SYSTEM. I. NEAR-INFRARED SPECTROSCOPY

The classification of HR 8799 as a γ Dor variable (Kaye et al. 1999) was secure by 1997. Periodicity of the variations from Hipparcos data (Aerts et al. 1998) and a long photometric campaign revealed multiple periodic signals with frequencies ranging from 0.2479 to 1.9791 days with amplitudes up to 1.5% and with a range of ±4% in the longest mode (Zerbi et al. 1999). The derived stellar properties from this study are L = 5.0 L_☉, T_eff = 7230 K, and R = 1.44 R_☉ (Baines et al. 2012).

Large time-series data sets of photometry and spectra have been taken on HR 8799 and a number of asteroseismology analyses have been conducted, first to constrain the age, unfortunately with only weak constraints, (Moya et al. 2010) and later to constrain the inclination angle of the star's rotation axis (i > 40°; Reidemeister et al. 2009; Wright et al. 2011). The rotational velocity of the star is vsin i = 37.5 to 49 km s⁻¹, based on a variety of different measurements (Zorec & Royer 2012; Kaye & Strassmeier 1998).

Though its age is still under debate, several indicators suggest that the star is younger than 100 Myr (Marois et al. 2010b; Doyon et al. 2010; Zuckerman et al. 2011; Malo et al. 2013; Moya et al. 2010). One estimate of this system's age comes from two independent studies (Doyon et al. 2010; Zuckerman et al. 2011), each of which used galactic space motions and spectroscopic age indicators to conclude that HR 8799 has a very high probability (∼98%) of being a member of the 30 Myr Columba Association. In contrast, Hinz et al. (2010) convincingly demonstrate that this association is poorly supported and that the age of HR 8799 cannot reliably be associated with Columba. Earlier works argued for a system age younger than 50–60 Myr (Marois et al. 2008; Reidemeister et al. 2009). While significantly different from the 1 Gyr estimate put forth by Moya et al. (2010) using asteroseismology data, the Moya et al. (2010) work does not explicitly exclude such a young age.

The age of the system is important because it constrains the masses of the planets when used in conjunction with theoretical cooling models and the measured luminosity of the objects. It is interesting to note, as discussed in Skemer et al. (2012), that some of the strongest constraints on the age of the HR 8799 system come from the dynamical simulations of the four planets, combined with their broadband photometry. As Goździewski & Migaszewski (2009) discuss, only a small fraction of three-body systems resembling HR 8799 can remain stable longer than ∼100 Myr. Furthermore, the stability of a four planet system for ∼10–100 Myr is consistent with lower masses for the companions (Goździewski & Migaszewski 2009; Reidemeister et al. 2009; Fabrycky & Murray-Clay 2010; Moro-Martín et al. 2010). Specifically, the study presented in Goździewski & Migaszewski (2009) finds only very limited dynamical configurations that foster long term stability: such stability will occur if the planets exhibit low-order, two- or three-body mean motion resonances.

These conclusions are consistent with the work presented by Fabrycky & Murray-Clay (2010), which find a longer system lifetime if the two innermost planets maintain a 2:1 commensurability. Furthermore, the persistence of this 2:1 resonance to the current epoch would put a specific constraint on the masses to be less than 10 M_Jup. Taking this one step further, Fabrycky & Murray-Clay (2010) point out that these dynamical constraints serve as early tests of the so called "hot start" models (Baraffe et al. 2003). Specifically, if all three planets exhibit low-order commensurabilities with each other, then the hot start models are easily consistent with the measured luminosities and inferred masses. On the other hand, if only the two innermost planets have a 2:1 commensurability, then the hot start models are only barely applicable. These studies point to the promise of using high contrast imaging to place empirical constraints on evolutionary models for the brightness of substellar objects (e.g., Crepp et al. 2012). Nonetheless, the published absolute H-band magnitudes, M_H ∼ 14 to 15 (Marois et al. 2008), are consistent with standard cooling evolutionary models for 30 Myr (Burrows et al. 1997).

Furthermore, the degree of excess infrared emission measured at 24 μm (Su et al. 2009) is consistent with an age significantly younger than 1 Gyr. Observations of A-stars with circumstellar material plus well-constrained ages show a marked decline with age in the strength of the 24 μm infrared excess (e.g., Rieke et al. 2005). This empirical result can be well-modeled by dynamical simulations (e.g., Wyatt et al. 2007), illustrating the clear decay of this hot dust emission due to the shorter dynamical time scales associated with the innermost portions of these systems.

Based on the ages of 30, 60, and 430 Myr, the four companions have masses in the following ranges: 7, 13, and 34 M_Jup for c, d, and e; and 5, 14, and 23 M_Jup for b using the cooling models of Baraffe et al. (2003). An important caveat exists regarding the age of the system, as discussed in Moya et al. (2010). If the inclination angle of the star is near 50° and the mass of the star is 1.45 M_☉, then the age range from asteroseismology analysis is 26 to 430 Myr, with the larger ages pushing the masses of all components into the "brown dwarf" mass range. As Moro-Martín et al. (2010) point out, the inclination of the star is probably smaller than the value initially assumed. Nonetheless, spectroscopy can also be used to constrain mass, or more accurately, log (g), as described in, for example, Marley et al. (2012) and Madhusudhan et al. (2011).

A number of studies have accomplished spectroscopy of the b component in the L, K and, H bands (Janson et al. 2010; Bowler et al. 2010; Barman et al. 2011). In addition, preliminary observations of the c component in the K band were presented by Konopacky et al. (2013). In Janson et al. (2010), a noisy L-band spectrum of b, ultimately sub-sampled to very low resolution, suggested that either improvements in the treatment of dust are required to model the object properly, or significant non-equilibrium chemistry is present in b's atmosphere. Bowler et al. (2010) obtained a K-band spectrum of b, and, combined with the Barman et al. (2011) data in H and K, a number of modeling papers, notably Marley et al. (2012), suggest that complex cloud structure, high metallicity, and/or non-equilibrium chemistry are needed to understand this object. In any case, b shows little or no methane absorption, possibly some CO, and the shape of the H-band peak suggested low-gravity Barman et al. (2011). Unfortunately, at the time of writing, the c component spectrum (Konopacky et al. 2013) and its analysis were not available in any publication.

These previously published measurements are all of the outer most planet in the system. However, this currently unique system has four directly detected planets. Thus, it is ideal for the study of multiple planets around a single star to investigate the diversity of planets. This is the purpose of this paper.

HR 8799 has become a benchmark target of observation in the field of direct exoplanet imaging to demonstrate whether a given project or technique is capable of finding even fainter objects, relative to the star, than the four planets orbiting this star. As a result, numerous images of the system have been obtained mainly in the near-IR through L and M bands (Hinz et al. 2010; Esposito et al. 2013; Currie et al. 2011, 2012; Hinkley et al. 2011a; Galicher et al. 2011; Janson et al. 2010; Lafrenière et al. 2009; Marois et al. 2008, 2010b; Serabyn et al. 2010; Skemer et al. 2012; Soummer et al. 2011). The conclusion of this work from various different instruments and telescopes indicates that the three outer planets are roughly coplanar, orbiting at an inclination angle of i ∼ 27.3 to 314, with eccentricities below 0.1 and semi-major axes of 68, 42, and 27 AU for b, c, and d, respectively (see especially Soummer et al. 2011; Currie et al. 2012). Orbital motion of e is poorly constrained at present, and only two groups attempted dynamical studies including all four planets (Esposito et al. 2013; Marois et al. 2010b).

The astrometry of this system is complicated by the background speckle field and, due to the extensive nature of the data reduction, we present that analysis in a second paper (Paper II; Pueyo et al. 2013, in preparation). However, HR 8799 will likely be the first exoplanetary system with a complete characterization including dynamical masses from astrometry through atmospheric analysis via spectroscopy, part of the purpose of combining imaging with spectroscopy in a single instrument such as Project 1640. We defer further discussion of the orbital characteristics to Paper II.

The HR 8799 system was first identified by Zuckerman & Song (2004) as a source with a prominent IRAS infrared excess, indicating the presence of cool dust. Shortly thereafter, the analysis of Spitzer spectroscopic data presented in Chen et al. (2006) indicated that the debris structure in this system was organized into two spatially distinct components. After the initial discovery of the planets in the system, Su et al. (2009) presented the analysis of much deeper Spitzer spectroscopy as well as photometry at 24 and 70 μm. In addition to constraining the spatial extent of the dust belts, this work identified an extensive outer halo of dust particles at hundreds of AUs, and these larger regions were later probed through observations presented in Hughes et al. (2011) and Patience et al. (2011) at 350 and 880 μm, respectively. These observations were used to resolve the outer debris structure surrounding the planets, solidifying the framework presented in Su et al. (2009). A coherent picture of the overall architecture of the system near the time of the planet's discovery, including the third outer component to the debris disk, was synthesized in Reidemeister et al. (2009)

The inner structure of the HR 8799 debris disk is shaped through dynamical interactions with the four planets. Specifically, the radial gap in the debris structure from 15–90 AU is clearly consistent with clearing caused by the four known planets. However, the source of clearing interior to 6 AU in the inner warm belt remains a mystery, with only Hinkley et al. (2011a) placing firm constraints in this region. If a planetary mass companion is responsible for the clearing, then Hinkley et al. (2011a) put a firm upper limit of ∼11 M_Jup on its mass, ruling out the possibility of a companion star or brown dwarf between 0.8 and 10 AU. Spectroscopy in the published literature has not revealed any evidence of binarity either.

The only companion object to HR 8799 with published spectral information in this study's wavelength range is component b (Barman et al. 2011). Our spectrum in the same wavelength region agrees with the Barman et al. (2011) spectrum, except for two points which deviate by roughly 2σ. Figure 5 shows both spectra overlaid to demonstrate the consistency of the two results. The current spectrum reveals two distinct features however. These are discussed in detail in Section 8.3.

To evaluate the S4 spectral extraction fidelity for each of the four components, we also extracted spectra at 6 different, randomly selected locations at the same angular radius from the central star. These spectra are essentially flat and do not reveal the detection of any object. Rather, they indicate the background against which the planets themselves are being detected. These are shown in Figure 4 by the dashed lines for each component. Note that in the case of d and e, these background spectra show that the components are only barely detected in the Y and J bands. However, integration of these channels permits the photometry for b and c and the upper limits for d and e.

Along these lines, we also derived spectra of 14 other similarly bright peaks in the detection map shown in Figure 3, to see whether they were real sources. All of them, which are either not point sources or are of lower significance than the bona-fide companions, either showed peaks in the insensitive water band between 1350 and 1430 nm—indicating spurious starlight contamination—or flat, featureless spectra at levels similar to the nearby background, indicating that they are not real sources.

As another test of the ability of our data analysis to retrieve spectra with high fidelity, we inserted multiple fake companion sources into the same data cubes we analyzed for this paper. The sources were placed at the same overall brightness as each component and at the same radii from the star. The fake sources were given the T4.5 spectrum of the standard 2MASS J0559−1404, chosen for its strong molecular features (Burgasser & Kirkpatrick 2006), which are useful to understand whether the features are correctly reproduced. The spectrum was resampled at Project 1640's resolution and wavelength range and applied to a cube of a fiducial PSF derived from an unocculted observation of the internal calibration white-light source. The instrumental response of the system was also applied to this PSF and then it was reduced in intensity by a factor sufficient to make it of the same intensity as each of the b, c, d, and e components. We then ran this composite data cube through the S4 algorithm, detected each object with the same significance as the real components, and extracted their spectra as described in Appendix B. This was repeated five times with different, randomly selected azimuthal locations for the fake objects at the same radii as the planets.

The spectra of fake sources extracted matched the input spectrum, with an rms error of less than 1σ for all four cases. All four spectra of the fake sources are shown in Figure 6. In general, the extractions are identical in shape to the input spectrum and in all cases reproduce small kinks in the input T4.5 object. For reference, each fake component was set to the same 2%–3% of the local speckle brightness as the real planets are (Section 2; i.e., the companions are roughly 40 times fainter than the local speckle background). The average deviation from the input spectrum over all wavelengths is 1.6%, 8.8%, 5.6%, and 15.0% for b through e, respectively, consistent with the photometry error values given in Table 2. The largest discrepancies appear in the fake source placed at the radial distance of e, as expected, due to its close proximity to the star. In particular, the J-band peak is over-estimated, although the broadband shape and slope of the spectra are retrieved to better than 20% accuracy. We note that all molecular features are reproduced, except for the very reddest part of the H band, which is noisy in any case (see Figure 4). The critical point with respect to interpreting the real spectra of the companions is that every minor feature (including, for example, the small dent in the J-band peak at about 1250 nm) is reproduced with fidelity in the simulated source spectra. Generally Y- and J-band data are worse because of poorer AO performance at those wavelengths.

Zoom In Zoom Out Reset image size

Because of the possibility of apparent systematic, but broad-band over or under estimates of the peaks, especially for the radial location of e, as described in the last section, we investigated the behavior of injected fake sources with perfectly white spectra positioned at various locations in the field of view. This test was designed to reveal any field-dependent biases in the spectral extraction by S4 and to ascertain whether a further minor correction to the spectra extracted is necessary. No such correction was deemed necessary, as we explain below.

We injected fake white-spectrum companions at a range of intensities (expressed as a fraction of the local speckle intensity) and at a range of radial distances from the star (from 15 to 50 pixels, in steps of 5 pixels). For each radial location, we injected fake sources with a range of 0.5%–5% intensity relative to the local speckle brightness, with steps of 0.5%. For each of these 140 fake sources, we retrieved the spectrum as we did with the tests in Section 7.3, and computed the rms error between the injected and retrieved spectra as well as the difference in slope between the two. None of the extracted spectra had a linear slope that deviates from white by more than 0.2 rms. To visualize the results of this test, we present a surface plot of the rms error versus radial location and brightness Figure 7. The rms error at the locations of each of the four planets is less than 0.4 in all cases, consistent with the tests utilizing a T4.5 spectrum (Section 7.3). No location in the explored parameter space exceeds an rms error of 0.75.

Zoom In Zoom Out Reset image size

Since no systematic error was detected and no additional correction of spectral shape is required, we conclude that the spectra shown in Figure 4 are real and as accurate as the data permit. Thus, it is appropriate to interpret them.

As a very rudimentary initial measurement, we calculated the statistical covariance among the six possible pairs of four spectra and the two spectra of c, to establish how similar or dissimilar they are from each other. Except for the comparison of the two epochs of the observation of c, no two have a covariance above 0.8 (a common measure of a real correlation) and all are below 0.2, a general benchmark for a significant statistical difference. Components b and c have a covariance of 0.152, which also indicates a significant statistical difference between the two objects, though their covariance is higher than in the other cases. The calculated values are tabulated in Table 4. For the calculations in which d or e were involved, only H-band data were used, since the Y and J spectra are only marginal detections.

The covariance values suggest that a very weak similarity exists between b and c (essentially the main hump in the H band). These values also suggest that the others are all significantly different from each other, confirming an inspection by eye of Figure 4.

For comparison purposes, we attempted to find spectra of other celestial point sources that exhibit at least some similarity to those of the components of HR 8799. We primarily used the library of spectra in Rayner et al. (2009). Initially, we plotted the four spectra overlaid with those of the outer planets of the solar system. However, the planets do not provide reasonable matches, with the sole exception of the e component as described below. We expanded the pool of comparison objects to stars and brown dwarfs. In Figure 8, we show comparisons with objects of best match presented for each of the four companions of HR 8799. In all cases, there are discrepancies, which clearly shows that these four objects are currently unique. Sources b and c seem to exhibit the usual water features between the Y, J and H bands, while for d and e we only detect the red edge of the water band between the J and H bands. The main methane feature at 1.65 μm seems present only in the d and e components, and strongest in e. See Section 8.3.

Zoom In Zoom Out Reset image size

In Figure 8 (top-left), we overlay the spectrum of the peculiar L5 brown dwarf 2MASS J035523.51+113337.4 (2M0355; Faherty et al. 2013) with the b component of HR 8799. This brown dwarf, sometimes referred to with a L5γ spectral type, seems best explained with low gravity and a possible reddening due to dust (Faherty et al. 2013). Although the general shape of the spectrum of 2M0355 may be a match, none of the small features are reproduced. We believe these to be molecular features and they are discussed in Section 8.3.

In the top-right panel of Figure 8, we compare HR 8799c to the T2 dwarf SDSS J125453.90−012247.4. Clearly, this object is not completely similar to a field T dwarf, as significant discrepancies appear at many wavelength channels. The comparison T dwarf is somewhat bluer, with HR 8799c showing a lower brightness in the Y band and a lack of methane in the H band, while the J-band peak appears sharper or "pointier." We note that the covariance of these two spectra is −0.41 which is dominated by the broad slope difference.

In the bottom-left panel of Figure 8, we overlay the spectrum of 2M0355 (Faherty et al. 2013) with the d component of HR 8799, as we did in the top-left panel for the b component. Important discrepancies exist in several features around 1500 nm and 1660 nm in the d spectrum. These are discussed further in Section 8.3.

The bottom-right panel of Figure 8 shows perhaps the best match between another object and one of the HR 8799 planets, the e component. Overlaid is the spectrum of Saturn, but reddened by a function of λ⁴ and normalized to match the overall flux level. This reddening function is really an attempt to approximate the spectrum of Saturn with the Sun's Rayleigh–Jeans black body spectral tail removed, as though it were a night-side spectrum of Saturn. We are unaware of any night-side spectra of the giant planets in this wavelength range, and so employed this simplistic approach. We note that there is a strong indication of methane in this object.

Recently, Marley et al. (2012) and Morley et al. (2012) have produced new models that incorporate far more complex and realistic cloud and condensate behavior. These have been applied to the previously published photometry of the HR 8799 system and indicate that clouds at varying levels in the atmosphere and with incomplete covering fraction of the atmosphere (as is seen in the giant planets of our solar system) may explain the red behavior of the spectral energy distributions. The new data presented here should provide additional constraints to these models.

Using the database of Freedman et al. (2008) and in private communications with M. Marley and D. Saumon (2012), we have indicated in Figure 4 the four main sources of opacity that most closely match the features we see in these four objects (aside from the water opacities between the astronomical bandpasses). These include ammonia, NH₃, in the 1450–1550 nm range, acetylene, C₂H₂, in the 1500–1550 nm range, methane, CH₄, redward of 1650 nm, and possibly CO₂ in the 1560–1630 nm range, a feature which is commonly used in terrestrial atmospheric science (e.g., O'Brien & Rayner 2002). Given the temperatures and masses of these objects as discussed in Section 1, these seemed the simplest molecules to consider. In Figure 4, we indicate the CO₂ feature with a question mark because it is unclear that it should be visible at the high temperatures inferred for these objects.

Based on the spectra and assuming that we have identified the sources of opacity correctly, we propose the following.

• 1.  
  b: contains ammonia and/or acetylene as well as CO₂ but little methane.
• 2.  
  c: contains ammonia, perhaps some acetylene but neither CO₂ nor substantial methane.
• 3.  
  d: contains acetylene, methane and CO₂ but ammonia is not definitively detected.
• 4.  
  e: contains methane and acetylene but no ammonia or CO₂.

We emphasize here that these are tentative identifications, although the lack of methane in the b component is consistent with the work in Skemer et al. (2012), Barman et al. (2011), and Bowler et al. (2010). In contrast, Barman et al. (2011) identify CO, not CO₂ as we indicate here. Detailed modeling by theorists will be necessary to definitively test these identifications.

Acetylene has never been convincingly identified in a substellar object outside the solar system, and unfortunately the fundamental band is in the thermal infrared, for which we have no spectral coverage. Therefore, it is difficult to confirm with other observations at this point.

As mentioned in Section 6, there is a weak 2σ difference in the c spectrum between 2012 June and 2012 October that is suggestive of the presence of a weak CO₂ feature at ∼1610 nm in 2012 October that is not apparent in 2012 June. Given the new models with patchy clouds which imply variability in these objects, it is not a leap of faith to consider variability of certain molecular features in these objects. It is probably unlikely that broadband variability in the primary star (as discussed in Section 1) is the cause of such a feature, mainly because the emergent flux from each of these planets is on the order of a few times 10⁷ ergs cm⁻² s⁻¹ (assuming the ∼800 K temperatures others have derived) while the incident starlight ranges from three to four orders of magnitude smaller.

However, the UV flux from the primary star is considerable because it is an A5V star—with more than 1000 times the UV flux of the Sun. If the Lyα line is also prominent in HR 8799, then there may be larger UV radiation incident upon even the most distant b component compared to that upon Jupiter. As Zahnle et al. (2009) demonstrate, UV flux combined with strong mixing in the atmosphere can greatly complicate the chemistry represented in the emergent spectra of planets (such as Jupiter, in their paper). Acetylene and even hydrogen cyanide (HCN) can become abundant. We mention HCN because it has an absorption feature at 1580 nm, comparable to the feature we have tentatively assigned to CO₂. A known FeH feature exists there as well, but has a broader extent.

We will continue to monitor the spectra of these planets into the future to assess the reality of any variability. We note that for years, photometric variability of L dwarfs has been known and detected in some T dwarfs as well (e.g., Khandrika et al. 2013 and references therein). From this point of view, it would not be surprising to detect variability in any or all of these planets.

As discussed in Section 1 and shown in Table 2, c, d, and e differ by less than a magnitude in the H band, but their spectra are not the same. Furthermore, b has some similarity in spectral shape when compared with c, which is nearly a magnitude brighter. This includes the peak near 1350 nm and the general shape of the large feature in H band. In addition, the e component has similarity with a reddened spectrum of Saturn. All these facts point to a diversity in the salient properties of objects for which mass, age, and metallicity are not the only factors determining observables.

Our results also highlight the power of discerning spectral differences between planets in multi-planet systems, which can provide constraints on the formation locations for each object. Specifically, Öberg et al. (2011) point out that the measured carbon-to-oxygen ratios in exoplanet atmospheres may constrain the formation locations of such objects, since various snow lines of carbon and oxygen-rich ices form at various radial locations from the star.

The results presented here consist of the first comparative spectroscopic study of multiple planets around a star other than our Sun.

Finally, we note that it would be of tremendous value to the exoplanet community for any space missions passing behind any of the planets of the solar system to take night-side thermal spectra of the planets in the near-IR, to enable direct comparisons with exoplanets.

The extracted data-cubes are first cleaned for remaining spurious bad pixels using median filtering and high-pass filtering to remove the contribution of the residual atmospheric halo. Their radial scaling and relative registration are then calculated using the procedure discussed in Crepp et al. (2011) and Pueyo et al. (2012a). For each wavelength from λ₀ = 995 nm, all of the subsequent slices in each cube of the observing sequence are then compressed or stretched and registered. In these pre-processed cubes, all speckle patterns appear at the same spatial scale and the PSF at λ₀ is the true image of the sky. This yields a series of N_λ pre-processed cubes of dimensions N_λ × N_Exposures × N_Spaxels × N_Spaxels. We then proceed to speckle removal using the KLIP algorithm described in Soummer et al. (2012). The target images consist of the N_Exposures PSFs at λ₀ in the pre-processed cube $n_{\lambda _0}$. We then proceed to partition the images into search zones, $\mathcal {S}$, that are equivalent to the optimization zones in LOCI as presented by Lafrenière et al. (2007; radial location r, radial width Δr, azimuthal location ϕ, azimuthal width Δϕ). Note that in contrast to LOCI, for KLIP these are not subtraction zones. For each target image, the reference ensemble is chosen as the subset of rescaled slices in the same exposure whose wavelength λ is such that the image of a putative planet in the compressed/stretched images is sufficiently far from its location at λ₀ (parameter N_δ in Lafrenière et al. 2007; Crepp et al. 2011; Pueyo et al. 2012a). Finally the N_Exposures speckle-reduced images at λ₀ are co-added and we proceed to the next wavelength.

We conduct a Monte Carlo search over the ensemble of parameters (N_δ, Δr, Δϕ, K_KLIP) with N_δ = 0.5, 0.75, 1, 1.25, 1.5 W (W is the un-occulted PSFs FWHM), Δr = 10, 20, 30 pixels, (Δϕ/2π) = (1/2), (1/4), (1/6), (1/8), (1/12), (1/16), (1/20), (1/24), and K_KLIP = 1.25. This approach is reminiscent of the parameter search in Soummer et al. (2011). However, when compared to LOCI, it is greatly facilitated by the absence of subtraction in KLIP and by our moderate number of reference PSFs. HR 8799 b and c are visible in H-band reduced images (S/N >8) and detected using matched filtering with un-occulted PSFs, at S/N  ∼  3, in J-band reduced images for (Δϕ/2π) < (1/6) and K_KLIP > 10. Since HR 8799 d and e lie at closer angular separations, they are detected in the H band by matched filtering at S/N  ∼  3 for (Δϕ/2π) > (1/8) (values of N_A, see Lafrenière et al. 2007, similar to the ones at the location of HR 8799 b, c) and K_KLIP > 10. These findings are summarized in the integrated H-band image shown in Figure 2.

Large biases of the spectral information can arise when estimating the spectro-photometry of faint companions unraveled by an aggressive PSF subtraction routine. In Pueyo et al. (2012a), we identified two sources of potential biases in IFS data: self-subtraction from over fitting the companion's signal, and wavelength-to-wavelength cross talk. When using KLIP, the former can be calibrated using the forward modeling methodology discussed in Soummer et al. (2012), provided that the self-subtraction is not too severe (e.g., provided that it does not radically change the morphology of the PSF). However, since our detection pipeline relies on search zones of large radial extent, some companion flux is present in the basis-set of principal components. Choosing a small K_KLIP can in principle mitigate the wavelength-to-wavelength cross talk, because the contribution of the companion's flux to the reference PSF is most likely in the small eigenvalue components. However, this does not fully alleviate this phenomenon. The full details of this characterization pipeline will be reported in a subsequent paper focusing on the astrometric characterization of the HR 8799 using Project 1640. Below, we only outline the main steps of our method.

Once the rough location (r_c, ϕ_c) of the companions is known from the detection pipeline, our spectral extraction pipeline is composed of the following steps.

• 1.  
  Minimization of wavelength-to-wavelength cross talk. We partition the images in such a way that ensures that there is no companion signal in the search zones. For the wavelengths λ₀ < λ_Mid, where λ_Mid = 1380 nm is the central wavelength of the Project 1640 bandpass, we choose search zones azimuthally centered at ϕ_c over the radial interval [r_c − N_δW, r_c − N_δW + Δr] and only use reference images with λ > λ₀. For the wavelengths λ₀ > λ_Mid, we choose search zones azimuthally centered at ϕ_c over the radial interval [r_c + N_δW − Δr, r_c + N_δW] and only use reference images with λ < λ₀. This ensures that there is no companion signal in the search zones but dramatically reduces the number of PSFs in the reference ensemble when only using slices from the same exposure. We alleviate this issue by aggregating slices at all the relevant wavelengths from the full observing sequence in a larger reference ensemble. We then proceed through KLIP reduction for these search zones over the range of parameters for which the companion has been previously detected.
• 2.  
  Calibration of the self-subtraction using forward modeling. Using un-occulted PSFs, we apply the forward modeling methodology described in Soummer et al. (2012). We use an un-occulted PSF at each wavelength as a template. For each reduction parameter, this process produces single wavelength detection maps of width nine pixels around (r_c, ϕ_c). We then derive the spectro-photometry in each channel as a value at the peak of an area of the detection map contained in a small circle of diameter thee pixels centered at (r_c, ϕ_c). The random spectro-photometric uncertainty is estimated as the scatter within the 63% confidence interval of the detection map around this peak. The rather large 9 × 9 spatial extent of the detection maps is chosen as a sanity check to rule out the cases of residual speckles contaminating the companion's signal. This is necessary to rule out the pathological case for which the tail of a speckle located ∼2 pixels away from the companion is considered to be companion flux by the matched filter.
• 3.  
  Estimation of the systematic uncertainty due to the reduction. The previous steps yield a series of spectra over a large collection of reduction parameters (N_δ, Δr, Δϕ, K_KLIP). Direct inspection of the detection maps determine whether or not the companions have been detected for a given set of parameters. For a given sub-set of (N_δ, Δr, Δϕ) where detection occurs, there are three regimes of PCA truncation: K_KLIP < K_KLIP, Min, for which residual speckles in the neighborhood of the companion are still visible; K_KLIP > K_KLIP, Max, for which detection occurs but the morphology of the PSF has been significantly altered by the self-subtraction and thus cannot be calibrated by forward modeling; and K_KLIP, Min < K_KLIP < K_KLIP, Max, which are of scientific interest. We analyze our extracted spectra in conjunction with single wavelength detection maps in order to establish the subset of parameters over which each planet is detected. The spectra are then derived as the mean of this subset of spectra and their uncertainty due to the reduction is derived as the scatter of this ensemble—for HR 8799 b and c, this subset is rather large (∼100 spectra) while it is smaller for HR 8799 d and e (∼10 spectra).
• 4.  
  Spectral calibration. The resultant spectra from the previous steps are expressed in fractional units of un-occulted core intensity of the primary star (e.g., units of contrast). We derive the spectrum of HR 8799 using the A5V template standard from the Pickles library (Pickles 1998) normalized at m_H = 5.28, or M_H = 2.30 (Cutri et al. 2003) and binned at the Project 1640 resolution. The final spectra are obtained by multiplying the results from the previous step with this stellar spectrum. Note that this procedure for spectral calibration is slightly different than in earlier Project 1640 results presented in Hinkley et al. (2010) and Roberts et al. (2012). Those relied on the ratio of aperture photometry estimates of both the PSF cores and companion. As a sanity check, we conducted an analysis analogous to Hinkley et al. (2010) and Roberts et al. (2012): we used the relevant reduced images selected in step 2 and conducted step 3. By projecting un-occulted PSFs onto the relevant principal components and using aperture photometry to calibrate the self-subtraction, we computed the spectral calibration using aperture photometry of un-occulted images. The two yielded results that are consistent within our uncertainties.

Figure 9 compares the KLIP spectral extraction with the S4 method.

Zoom In Zoom Out Reset image size

The S4, algorithm for post-processing speckle suppression (Fergus et al. 2013), based on principal component analysis (PCA), has not yet been published. Thus, we provide some details of it here.

S4 takes as input a 4D data block of dimension $N_{\rm pixels_y} \times N_{\rm pixels_x} \times N_{\lambda } \times N_{\rm Exposures}$ which has been pre-processed as follows: (1) application of a 2 × 2 median filter to each band/exposure to remove dead pixels; (2) spatial alignment of all bands and exposures relative to one another; and (3) spatial centering so that the star lies precisely (within 0.1 pixel) at $(N_{\rm pixels_y}/2,N_{\rm pixels_x}/2)$.

In S4, the data at each spatial position is decomposed into a speckle component and a companion component, whose sum reconstructs the original data to the limit of Gaussian noise. To model the speckle component effectively, which evolves radially with wavelength, this decomposition is performed in a polar reference frame, as shown in Figure 10 (left). To examine a location at radius d and angle θ from the center, an annular region of width R at a radius d is transformed into a region of size Θ × R × N_λ × N_Exposures, where Θ = 2πd to ensure the region is well-sampled. This region is divided into a $test$ zone around the location θ (of size δθ) and a $training$ zone of all other angles, as shown in red and green, respectively, in Figure 10 (left). We use the training zone to build a model of the speckles that is then applied to the test region, decomposing it into speckle and companion components.

Zoom In Zoom Out Reset image size

We treat the speckles in both the training and test zones as being independent over angle θ and exposure n, justified by the limited angular extent of the speckles and their variation due to atmospheric turbulence and instrument flexure or temperature variations between exposures. We thus assume the structure of the speckles to be confined to a joint radius–wavelength space (of dimension λR), as illustrated in Figure 10 (right). This joint space is modeled using PCA, which approximates each λR-dimensional slice of the data at angle θ cube n as a linear combination of K orthonormal basis vectors W = [w₁, ..., w_K] (K being a user-defined parameter). These basis vectors capture the majority of the variance in the radius–wavelength space and are computed by performing an eigendecomposition of the covariance matrix built from the training region, reshaped into a matrix of λR dimensions by (Θ − δθ)N_Exposures samples. We then use the PCA basis W to infer the speckle component of the test zone. W is visualized in Figure 10 (upper right).

In detection, we first use the PCA basis W to fit the data in the test zone and then correlate the residual error with a fixed companion model P (previously obtained by calibration with a point white light source). The response for position (θ, d) is stored in a correlation map. The intuition is that the PCA basis W will effectively model the speckles, but not the companion, thus the residual should only contain the companion signal which will respond strongly to the matched filter P (see Figure 10, lower right). The location of the test zone is systematically moved across all angles θ at a given radius, with the PCA basis being recomputed at each location, because the training zone also changes. The process is then repeated for a new radius d, until all spatial locations in the visual field have been examined.

Finally, a normalization is performed on the correlation map that compensates for the variation in flux with radius in the original data. The resulting normalized map is then converted back to the Cartesian coordinate frame and is the output of the detection algorithm. (See Figures 3 and 10(c) for examples.)

Promising peaks in the normalized detection map are then selected for spectrum extraction. This follows the same overall modeling approach as detection, except that the spectrum of the companion model P is no longer fixed to be white. For a peak at location (d, θ), let the observed data for exposure n be y_n (represented as λR dimensional vector). We now must estimate both the spectrum of the planet m (an N_λ dimensional vector) and the K dimensional PCA coefficients z_n for exposure n that reconstruct y_n. Using a Gaussian noise model, this is equivalent to minimizing the following convex objective:

$$\begin{equation} \sum _{n=1}^{N_{\rm Exposures}} \Vert y_n - (Wz_n + Pm)\Vert ^2_2, \end{equation} \tag{ B1 }$$

where Wz_n is the speckle component for exposure n and Pm is the companion component (constant across all exposures). We impose a non-negativity constraint on m, since negative spectra are not physically plausible. Minimization is performed using standard optimization software in Matlab. The spectra shown in Figure 4 are the resulting m vectors for the four different planet locations. An estimate of the uncertainty is obtained by measuring the variance in m over different settings of the model parameters, K. In the case of the data presented here, the number of principal components was varied from 50 to 500 with steps of 25, and the angular fitting size, Θ, was 3, 5, and 7 pixels. The error bars were calculated using principal components of 150, 175, and 200 and Θ of 5 and 7 pixels, settings which yielded the best results. For October, due to significantly worse observing conditions, more principal components were needed for detection and spectrum extraction: 250, 300, and 350 with Θ of 5 and 7 pixels.