Where Is the Water? Jupiter-like C/H Ratio but Strong H₂O Depletion Found on τ Boötis b Using SPIRou

The high temperature that enables water to exist in vapor form in the upper atmosphere of hot Jupiters is in sharp contrast to the solar system gas giants. Jupiter and Saturn are comparatively much colder and therefore have some of their atmospheric constituents (such as H₂O) condensed out to deep layers beneath the cloud deck, where they cannot easily be probed. For example, Jupiter was found to have a factor of ∼3 enhancement of C, N, S, P, Ar, Kr, and Xe compared to solar by the Galileo mission, but its H₂O content, and hence the O/H ratio, remains poorly constrained (Niemann et al. 1998; Owen et al. 1999; Atreya et al. 2018). The O abundance was only recorded as a lower limit of 0.3× solar because the Galileo probe entered an anomalously dry region of Jupiter's atmosphere and still recorded rising abundances before it stopped transmitting (Orton et al. 1998; Atreya & Wong 2005). More recently, the Juno satellite was used to probe deeper layers and infer an oxygen abundance of ${2.7}_{-1.7}^{+2.4}$ relative to solar in the 0° to +4° latitude equatorial regions of the Jovian giant, although a value of zero cannot be ruled out at the 2σ level (Li et al. 2020). Moreover, the authors caution that Jupiter exhibits large latitudinal variations, and so its equatorial regions may not be representative of the average composition. Obtaining precise constraints on the bulk oxygen content is crucial for comparing Jupiter's C/O ratio relative to the solar value of 0.54 (Asplund et al. 2009).

Outside of our solar system, H₂O absorption features have been detected on a large sample of hot Jupiters, although often at weaker levels than would be expected from cloud-free solar compositions (e.g., Madhusudhan et al. 2014b; Sing et al. 2016). While clouds can play a significant role in damping the observed signals, current retrieval studies suggest that subsolar H₂O abundances remain necessary to fully explain the suite of observed water features in the majority of cases (Barstow et al. 2017; Pinhas et al. 2019; Welbanks et al. 2019). It is, however, harder to place these H₂O measurements in context from a chemical standpoint without accompanying abundances of other molecules, such as CO or CH₄. Unfortunately, these are more difficult to constrain with the limited wavelength coverage of the Hubble Space Telescope Wide Field Camera 3 (HST/WFC3) data sets, from which almost all current hot Jupiter H₂O measurements are made. This then raises the following question: is the trend of subsolar water abundances observed on hot Jupiters inherently due to low overall metallicities, or rather the results of high C/O ratios, in which case oxygen would preferentially be bound in CO rather than H₂O (Madhusudhan 2012)?

The possibility of high C/O (carbon-rich) planets was first explored by Madhusudhan et al. (2011), who reported a C/O ≥ 1 by analyzing the dayside spectrum of the hot Jupiter WASP-12b. However, the high C/O claim rests primarily on a single data point, the photometric Spitzer Space Telescope secondary eclipse depth at 4.5 μm, and has widely been both contested (Crossfield et al. 2012; Line et al. 2014) and reaffirmed (Föhring et al. 2013; Stevenson et al. 2014). Meanwhile, the analysis of follow-up transit observations using HST/WFC3 favors C/O values closer to solar (Swain et al. 2013; Benneke 2015; Kreidberg et al. 2015), although carbon-rich scenarios are still plausible under different model parameterizations in some cases (Kreidberg et al. 2015).

Comparing instead planetary measurements to that of their host stars, Brewer et al. (2017) found that hot Jupiters generally have superstellar C/O ratios, although most cases remain consistent with the C/O value of their host star to 1σ owing to large measurement uncertainties. In a separate study, Welbanks et al. (2019) obtained constraints on the water and alkali metal abundances of 19 exoplanets and found prevalent supersolar Na and K abundances in conjunction with low H₂O concentrations on hot Jupiters. This combination would suggest that superstellar C/O, Na/O, and K/O ratios exist, and that oxygen is depleted relative to other species. While large error bars prevent definitive conclusions, the general trend found in these studies favors high C/O ratios as being common in hot Jupiter atmospheres.

The difficulty in obtaining precise C/O measurements for exoplanets with current space telescopes lies primarily in their limited wavelength coverage, which currently does not enable us to obtain spectroscopy past 1.7 μm to directly probe the prevalent bands of the main carbon-bearing molecules. This limits us to broadband photometric measurements beyond 1.7 μm, which are inherently difficult to interpret, rendering conclusions often sensitive to model assumptions (e.g., Kreidberg et al. 2015; Spake et al. 2020).

Circumventing the need for space-based low-resolution instruments and partially overcoming these limitations has proven to be possible in some cases for directly imaged planets, providing C/O measurements ranging from being subsolar or near solar (Konopacky et al. 2013; Barman et al. 2015; Todorov et al. 2016; Lavie et al. 2017; Mollière et al. 2020; Nowak et al. 2020; Wang et al. 2020) to being supersolar close to unity in the case of HR 8799b (Lee et al. 2013; Lavie et al. 2017). However, direct imaging is currently only applicable to a small sample of nearby systems with young giant planets that orbit at very wide separations.

In the past decade, high-dispersion cross-correlation spectroscopy (HDCCS) has emerged as an extremely powerful alternative method for identifying molecular species in hot Jupiter atmospheres. Key to the technique is the fact that short-period planets experience large line-of-sight velocity changes that induce Doppler shifts corresponding to multiple resolution elements per hour on a high-resolution spectrograph. This can allow for atmospheric signatures of exoplanets to be disentangled from contributions of both their host star and Earth's transmittance, which, in contrast, are essentially stationary in wavelength (Birkby 2018). HDCCS has successfully been used to produce a wealth of unambiguous molecular detections in both transiting and nontransiting exoplanet atmospheres using CRIRES (Snellen et al. 2010; Brogi et al. 2012, 2013, 2014, 2016, 2017; Rodler et al. 2012; Birkby et al. 2013, 2017; de Kok et al. 2013; Schwarz et al. 2015; Hawker et al. 2018; Cabot et al. 2019; Webb et al. 2020), Keck NIRSPEC (Rodler et al. 2013; Lockwood et al. 2014; Piskorz et al. 2016, 2017, 2018; Buzard et al. 2020), Subaru HDS (Nugroho et al. 2017), GIANO (Brogi et al. 2018; Guilluy et al. 2019; Giacobbe et al. 2021), CARMENES (Alonso-Floriano et al. 2019; Sánchez-López et al. 2019, 2020), IGRINS (Flagg et al. 2019), and Subaru IRD (Nugroho et al. 2021).

A cross-correlation analysis can be insightful in determining whether an exoplanet atmosphere contains a certain species or a thermal inversion. However, one area in which HDCCS has shown more limited success is in retrieving molecular abundances and temperature structures in a robust Bayesian fashion. This is a consequence of observing through Earth's atmosphere and the detrending procedure necessary for analyzing such data sets, which removes all continuum information, framing all time-varying effects of the planet's atmosphere in relative terms with respect to the stellar flux. Such a treatment of the data effectively prevents the use of standard "data minus model" likelihood calculations used in conventional atmospheric retrieval frameworks (Brogi & Line 2019). As a chi-squared-based approach is also nonideal for model comparison to data that may have leftover broadband variations or correlated noise from poorly removed telluric lines, the challenge is then to obtain a robust goodness-of-fit estimator from the more outlier-resistant cross-correlation function (CCF). Frameworks to perform retrievals on high-resolution data sets have since been proposed and implemented, making use of CCF-to-likelihood mappings (Brogi & Line 2019; Gandhi et al. 2019; Gibson et al. 2020). Applied on CRIRES observations, this has led to more robust atmospheric constraints on H₂O and CO abundances, even with only the ≤0.1 μm span of wavelength covered by such data sets. Results from these earlier Bayesian inference analyses have shown evidence for water abundances that are underabundant relative to carbon monoxide (Brogi & Line 2019; Gandhi et al. 2019). For time series observations taken with wider spectral range instruments with hundreds of thousands of pixels per exposure, performing full-scale high-resolution atmospheric retrievals remains a computationally challenging feat.

In this work, we present the detection and characterization of τ Boo b's dayside atmosphere using the new SPIRou high-resolution spectrograph and the SCARLET atmospheric retrieval framework adapted for HDCCS data sets. In Section 2 we describe the observations and the data reduction and detrending procedures. We discuss the modeling and adopted likelihood prescription in Section 3 and the cross-correlation and retrieval analyses in Section 4. In Section 5 we explore possible nonastrophysical causes for the inferred water abundance upper limit. Finally, we discuss the implications of our results for planet formation in Section 6 and conclude in Section 7.

The hot Jupiter τ Boo b was one of the first exoplanets ever discovered (Butler et al. 1997) and remains one of the brightest and nearest exoplanet-bearing systems known to date. It was the first exoplanet to be successfully observed at high resolution in thermal emission, with CO absorption detected on its dayside hemisphere (Brogi et al. 2012; Rodler et al. 2012). With evidence of water vapor also having been reported (Lockwood et al. 2014), τ Boo b presents one of the best opportunities for atmospheric characterization and obtaining a precise measure of its C/O ratio. The τ Boo planetary system is nontransiting (Baliunas et al. 1997) and also hosts a widely separated M-dwarf companion (Hale 1994), with all three components orbiting in a coplanar configuration (Justesen & Albrecht 2019). Systemic parameters of τ Boo b and its host star are given in Table 1.

Table 1.  τ Boo System Properties

Stellar Parameters	Value
Stellar type (1)	F7V
K magnitude (2)	3.36
Stellar mass (M*) (3)	1.38 ± 0.05 M⊙
Stellar radius (R*) (3)	1.42 ± 0.08 R⊙
Luminosity (L*) (3)	3.06 ± 0.16 L⊙
Temperature (T*) (3)	6399 ± 45 K
Metallicity (Fe/H) (3)	0.26 ± 0.03
Age (3)	0.9 ± 0.5 Gyr
Distance (4)	15.66 ± 0.08 pc
Systemic velocity (Vsys) (5)	−16.4 ± 0.1 km s−1
Velocity semiamplitude (K*) (3)	471.73 ± 2.97 m s−1
Planet parameters	Value
Orbital period (P) (3)	3.3124568(69) days
Semimajor axis (a) (3)	0.049 ± 0.003 au
Eccentricity (e) (3)	0.011 ± 0.006
Longitude of periastron (ω) (3)	1134 ± 322
Epoch of periastron (T0) (3)	2,456,400.94 ± 0.30 (BJD)
Inclination (i) (6)	435 ± 07
Planet mass (Mp ) (6)	6.24 ± 0.23 MJup

References. (1) Gray et al. 2001; (2) Belle & Braun 2009; (3) Borsa et al. 2015; (4) Gaia Collaboration et al. 2018; (5) Brogi et al. 2012; (6) this work (Section 4.2).

Download table as:  ASCIITypeset image

We observed the extremely bright naked-eye τ Boötis system with SPIRou on five different nights at orbital phases targeting the dayside of τ Boo b (Figure 1). We obtained continuous 5 hr observation sequences on the nights of 2019 April 15 and 21 (Program 19AC29/19AC12, PI Pelletier/Herman), with additional 3.3 hr sequences obtained on 2020 May 10, June 2, and June 5 (Program 20AC32, PI Pelletier). The epochs were chosen to sample a wide portion of τ Boo b's orbit while also avoiding regions where a larger fraction of the colder nightside is in view. During the observations, the line-of-sight component of the planet's velocity shifts by approximately 6–8 km s⁻¹ per hour, corresponding to roughly 3 full pixels on the SPIRou detector. The observations were scheduled such that the radial velocity shift of τ Boo b's atmospheric trail does not overlap with the telluric rest frame, which could otherwise have resulted in leftover telluric residuals overlapping and interfering with the planetary signal of interest (e.g., Brogi et al. 2018; Alonso-Floriano et al. 2019; Sánchez-López et al. 2020).

Zoom In Zoom Out Reset image size

For each night we obtained 90 s exposures, averaging a signal-to-noise ratio (S/N) per pixel of ∼265, with values typically ranging between 100 and 400 depending on the spectral order, with peak values being obtained in the H and K bands. For these observations we opted not to use the polarimetric mode of SPIRou that measures the Stokes parameters. While in theory this option provides polarimetric information for "free," it requires the position of the rhombs to be moved after each exposure. To avoid having elements in the optical path move unnecessarily and to maximize stability, the rhombs were maintained at a fixed position for all observations. Sky conditions were dry and photometric throughout all the observations, with only the night of 2020 June 2 having elevated humidity levels.

We reduced the raw SPIRou data using A PipelinE to Reduce Observations (APERO) v0.5 (N.J. Cook et al. 2021, in preparation), which outputs wavelength-calibrated and science-ready spectra for each order of each exposure. APERO also provides the data in an already telluric-corrected format, including the reconstructed transmission spectrum used to correct each exposure. This reconstruction of Earth's transmittance is done via principal component analysis (PCA) using a library of SPIRou observations of telluric standards collected at various air masses and water vapor concentrations over many nights (E. Artigau et al. 2014a, 2021, in preparation). For our analysis, we use the non-telluric-corrected version of the data and instead remove tellurics from the data themselves as is commonly done in the literature (e.g., Snellen et al. 2010). This is done out of precaution, to avoid the automated removal of tellurics that may unexpectedly affect the underlying planetary signal in the data. We note, however, that starting from the pre-telluric-corrected data product and performing the same detrending produces near-identical results and does not change the conclusions of this paper. It may be beneficial to use this data product for SPIRou time series with planetary trails that overlap in radial velocity shift with the telluric rest frame. For our analysis we use the data in their extracted order-by-order flat-fielded format as provided by APERO. An example time series for one order of the 2019 April 15 night of observations is shown in the top panel of Figure 2.

Zoom In Zoom Out Reset image size

2.3.1. Bad-pixel Correction

2.3.2. Blaze Removal and Continuum Alignment

Before attempting to uncover the planetary signal, it is crucial to remove unwanted contributions from Earth's atmosphere and the host star. Many techniques exist to achieve such a "cleaning" of the data, all relying on the key concept that, contrary to Earth's and the star's atmospheres, which are stationary or quasi-stationary in wavelength, the exoplanet atmospheric features are rapidly being Doppler-shifted in time (Birkby 2018). With an orbital period of just 3.31 days, τ Boo b's orbital line-of-sight velocity accelerates by tens of kilometers per second over the course of the 3–5 hr long observations, moving across multiple pixels on the SPIRou detector. Therefore, the removal of all time-varying effects that are constant in wavelength/pixel can get rid of nearly all traces of both Earth's atmosphere and the star, while leaving the planetary signature mostly unaffected.

2.4.1. Median Spectrum and Time Polynomial Fits

2.4.2. Principal Component Analysis

2.4.3. Masking

We generate synthetic spectra of τ Boo b's atmosphere using the SCARLET framework (Benneke & Seager 2012, 2013; Knutson et al. 2014; Kreidberg et al. 2014; Benneke 2015; Benneke et al. 2019a, 2019b). SCARLET calculates high-resolution line-by-line thermal emission models for a given composition, thermal structure, and cloud properties. Molecular opacities and associated line lists used in this work include H₂O (Barber et al. 2006; Polyansky et al. 2018), CO (Rothman et al. 2010), CH₄ (Yurchenko & Tennyson 2014), CO₂ (Rothman et al. 2010), HCN (Harris et al. 2006; Barber et al. 2014), TiO (McKemmish et al. 2019), C₂H₂ (Chubb et al. 2020), and NH₃ (Yurchenko et al. 2011). Cross sections for these molecules are computed using the HELIOS-K opacity generator (Grimm & Heng 2015; Grimm et al. 2021) and are accessible from the open-access DACE database. ¹¹ All included molecules have prominent cross sections in the SPIRou bandpass and could be detectable if they are present in high-enough volume mixing ratios (VMRs) in τ Boo b's atmosphere (Figure 3). We do not add opacities from additional molecules or alkali metals, as they are not expected to have noticeable effects in our data. H₂–H₂ and H₂–He collision-induced absorption is computed following Borysow (2002). Models also include an optically thick gray cloud deck at a given cloud top pressure P_c in bars. Assuming that no other major carbon- or oxygen-bearing gaseous species are present in τ Boo b's atmosphere, the gas-phase C/O ratio can then be calculated as

$$\begin{eqnarray}&&{({\rm{C}}/{\rm{O}})}_{\mathrm{gas}}=\displaystyle \frac{\mathrm{CO}+{\mathrm{CH}}_{4}+{\mathrm{CO}}_{2}+\mathrm{HCN}+2{{\rm{C}}}_{2}{{\rm{H}}}_{2}}{{{\rm{H}}}_{2}{\rm{O}}+\mathrm{CO}+2{\mathrm{CO}}_{2}+\mathrm{TiO}}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

Given molecular abundances, a temperature–pressure (TP) profile, and a cloud top pressure, all spectra are computed via line-by-line radiative transfer at a spectral resolution of R = λ/Δ λ = 250,000 and subsequently convolved with a Gaussian profile to match the R = 70,000 instrumental resolution of SPIRou. The planetary thermal emission models F_p (λ) are then scaled to the blackbody continuum level of the star and planet-to-star area ratio

$$\begin{eqnarray}&&{F}_{\mathrm{scaled}}(\lambda )=\displaystyle \frac{{F}_{p}(\lambda )}{\pi B(\lambda ,{T}_{* })}{\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where B(λ, T_*) is the Planck function at the stellar effective temperature T_*, while R_p and R_* are the planetary and stellar radii, respectively (Table 1). As τ Boo b is a nontransiting planet, its radius is not well constrained, and thus we follow previous work and assume that R_p = 1.15 R_Jup, the average radius of all hot Jupiters in τ Boo b's mass regime measured from transiting planets (Brogi et al. 2012; Hoeijmakers et al. 2018).

While a generated forward model of τ Boo b's atmosphere may be a good description of the true astrophysical signal, it is not necessarily representative of the signal contained in the data after the removal of telluric and stellar contamination. This is because the detrending procedure that is applied to the data in Section 2.4 to remove telluric and stellar features also stretches and scales the planetary signature (Brogi & Line 2019; see also their Figure 2, step 7—noiseless). This warping must be accounted for to produce an accurate atmospheric analysis and avoid introducing any biases when retrieving molecular abundances.

To account for the distortions experienced by the underlying atmospheric signal during the telluric removal process, we replicate the same detrending procedure that was applied to the data onto each model. To achieve this, we first shift the scaled planetary model F_scaled(λ) in velocity space. For a given combination of K_p and V_sys, the radial velocity shift V_p (t) that τ Boo b would have at the observed epoch of each exposure is given by

$$\begin{eqnarray}&&{V}_{p}(t)={K}_{p}\left(\cos (f(t)+\omega )+e\cos \omega \right)+{V}_{\mathrm{sys}}+{V}_{\mathrm{bary}}(t),\end{eqnarray} \tag{ 3 }$$

where ω is the longitude of periastron, e is the eccentricity, V_bary(t) is the barycentric velocity of Earth, and f(t) is the true anomaly at the time of each observation. This produces the expected planet-to-star contrast level F_scaled(λ, t) Doppler-shifted as a function of time for a given atmospheric model. We store the detrending steps outlined in Section 2.4 that were divided out of the data to remove contributions from both the star F_s (λ) and Earth's transmittance T_E(λ, t) and assume that the observed flux is of the form

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}(\lambda ,t)={F}_{s}(\lambda ){T}_{{\rm{E}}}(\lambda ,t)\left(1+{F}_{\mathrm{scaled}}(\lambda ,t)\right).\end{eqnarray} \tag{ 4 }$$

The reconstruction of the processing steps applied to the data, consisting of the fitted median spectrum, time polynomial, and principal components, then acts as an approximation of F_s (λ)T_E(λ, t). We can therefore inject each generated forward model into this reconstructed F_s (λ)T_E(λ, t), as per Equation (4), and reapply the same detrending steps as was done to the data in Section 2.4. The result is then effectively 1 + F_scaled(λ, t), but as the planet spectra would appear after having gone through the same detrending process as the real data. Equation (4), but passed through the detrending procedure, is then what serves as model input for the likelihood evaluation.

With the planetary model (m) now representative of what the atmospheric signal would be in the processed data (d), we calculate the likelihood (L) via the CCF-to-lnL mapping

$$\begin{eqnarray}&&\mathrm{ln}(L)=-\displaystyle \frac{N}{2}\mathrm{ln}\left({s}_{d}^{2}-2{{as}}_{d}{s}_{m}(v)\mathrm{CCF}(v)+{a}^{2}{s}_{m}^{2}(v)\right)\end{eqnarray} \tag{ 5 }$$

derived in Brogi & Line (2019). Here a is a scaling factor and the CCF is given by

$$\begin{eqnarray}&&\mathrm{CCF}(v)=\displaystyle \frac{R(v)}{{s}_{d}{s}_{m}(v)},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&R(v)=\displaystyle \frac{1}{N}\sum _{n}d(n)m(n-v),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{s}_{d}^{2}=\displaystyle \frac{1}{N}\sum _{n}{d}^{2}(n),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{s}_{m}^{2}(v)=\displaystyle \frac{1}{N}\sum _{n}{m}^{2}(n-v),\end{eqnarray} \tag{ 9 }$$

for N spectral channels (Zucker 2003). The CCF is then calculated at a given velocity shift (v) for each order of each exposure by summing over each wavelength bin (n). Both d and m are mean normalized immediately before the CCF is computed. Keeping a as a free parameter relaxes the presumption that the imposed planet-to-star line contrast is representative of the true signal. In particular for a nontransiting planet such as τ Boo b, for which the precise planetary radius is unknown, fitting a simultaneously is crucial to encompass modeling uncertainties.

With the CCF-to-lnL mapping and an appropriate treatment of the generated models to replicate the detrending procedure added to the SCARLET framework, we can now perform atmospheric retrievals on high-resolution spectroscopy data sets. We opt for classical "free" retrievals that fit for molecular abundances directly to avoid making any assumptions on the chemistry on τ Boo b. Parameters explored in this analysis include molecular abundances, temperature structure points, cloud properties, orbital parameters, and the scaling factor a. We adopt log-uniform priors on molecular VMRs between 10⁻¹⁰ and 1, on the cloud top pressure P_c between 0.1 mbar and 100 bars, and for the scaling parameter a between 0.1 and 10. Linear-uniform priors are used for K_p and V_sys between +50 and −50 km s⁻¹ from their respective known values.

The precise shape of spectral lines probed in exoplanetary thermal emission spectra depends on the temperature profile of the planet. For example, if a planet has a thermal inversion, with temperature increasing with altitude, spectral features can be observed in emission rather than absorption. Similarly, the steepness of the temperature gradient dictates the narrowness of spectral lines. Hot Jupiters in particular can have different temperature structures depending on the presence of ultraviolet and visible opacity sources driving at what altitude the stellar irradiation is absorbed (Arcangeli et al. 2018; Merritt et al. 2020; Nugroho et al. 2020). It is thus crucial that, when performing an atmospheric retrieval on high-resolution hot Jupiter thermal emission observations, an appropriate TP profile parameterization able to capture all possible spectral line shapes is adopted.

We parameterize the TP profile with N free temperature points T uniformly distributed in log pressure between p_top = 10⁻⁷ and p_bot = 10² bars, the respective pressures at the top and bottom of the atmosphere. We impose a uniform prior between 0 and 5000 K on each temperature point and a Gaussian prior of mean zero and variance ${\sigma }_{s}^{2}$ on the rms of the second derivative of the profile. Dropping constant terms, this smoothing prior on the TP profile (π) is given by

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}\pi (T(\mathrm{log}p))=-\displaystyle \frac{1}{2{\sigma }_{s}^{2}}{\left(\mathrm{rms}\left(\displaystyle \frac{{d}^{2}T}{d{\left(\mathrm{log}p\right)}^{2}}\right)\right)}^{2}\\ =\displaystyle \frac{-1}{2{\sigma }_{s}^{2}\mathrm{log}({p}_{\mathrm{bot}}/{p}_{\mathrm{top}})}{\displaystyle \int }_{\mathrm{log}{p}_{\mathrm{top}}}^{\mathrm{log}{p}_{\mathrm{bot}}}{\left(\displaystyle \frac{{d}^{2}T}{d{\left(\mathrm{log}p\right)}^{2}}\right)}^{2}d(\mathrm{log}p),\end{array}\end{eqnarray} \tag{ 10 }$$

where rms is the root mean square of the continuous second derivative of the TP profile T(log p). This smoothing prior effectively punishes TP profiles that have high second derivatives and prevents nonphysical zig-zaggy temperature oscillations in pressure regions that are less well constrained by the data. The standard deviation σ_s is physically motivated and has units of kelvin per pressure decade squared (K dex⁻²), i.e., the unit of the second derivative. It controls how strongly the prior punishes the curvature in the temperature structure. Large values of σ_s allow for more curvature and deliver profiles with high vertical resolution. Assuming grid points uniformly distributed in log pressure, Equation (10) can be written in discretized form as

$$\begin{eqnarray}&&\mathrm{ln}\,\pi ({\boldsymbol{T}})=\displaystyle \frac{-1}{2{\sigma }_{s}^{2}\mathrm{log}({p}_{\mathrm{bot}}/{p}_{\mathrm{top}})}\sum _{i=2}^{N-1}\displaystyle \frac{{\left({T}_{i+1}-2{T}_{i}+{T}_{i-1}\right)}^{2}}{{\left({\rm{\Delta }}\mathrm{log}p\right)}_{i}^{3}},\end{eqnarray} \tag{ 11 }$$

with the summation over each layer i in the pressure interval between p_top and p_bot.

A curvature-punishing smoothing prior was previously used for brown dwarf retrievals (Line et al. 2015); however, their parameterization using a γ parameter was variant across different pressure grids, and moreover, we argue that σ_s should not be simultaneously fit for most exoplanet applications. In particular, fitting σ_s often leads to overly smooth profiles where extrapolation effects may falsely deliver tight constraints in pressure regions not probed by the data. This is particularly true for highly irradiated hot Jupiters whose dayside TP profiles are not as smooth as those of brown dwarfs. Instead, we argue that the σ_s parameter is effectively a user choice that sets the trade-off between the vertical spatial resolution and the noisiness and stability of the temperature structure. For τ Boo b, we explore different σ_s values between 50 K dex⁻² (strong smoothing) and 1000 K dex⁻² (weak smoothing). We find the smoothness prior to be too stringent for σ_s ≲ 80 K dex⁻², resulting in overly smooth temperature profiles that propagate constraints from the well-constrained mid-atmosphere into the deeper layers of the atmosphere. On the other extreme, we find nonphysical oscillations in pressure regions not readily probed by the data for σ_s ≳ 500 K dex⁻². For our final analysis we opt for a value of σ_s = 240 K dex⁻², which we find for this data set to be the best compromise between vertical resolution of the temperature structure and the noisiness and stability of the profile. We note that the retrieved abundances do not vary significantly under different adopted smoothing factor values.

The advantages of our temperature retrieval prescription for exoplanets are that it is independent of the pressure grid and that it does not make any prior assumption on the overall shape of the temperature structure, which is generally a disadvantage of parametric forms for the TP profile (e.g., Madhusudhan & Seager 2009; Guillot 2010; Benneke & Seager 2012; Line et al. 2012, 2013, 2016). We set the number of layers to N = 15, enough to properly sample the TP profile for our τ Boo b data set, while keeping the number of retrieval parameters small. We note that our physically motivated formulation of σ_s in K dex⁻² is independent of the number of pressure levels, and we obtain the same results by introducing more pressure levels. Too few temperature points, however, can result in an overly coarse TP profile being retrieved.

Performing full Bayesian atmospheric retrievals on HDS data sets is a computationally challenging task. Relative to low- or medium-resolution spectroscopic data sets from space-based instruments, the added computational cost of running HDS retrievals comes primarily from both the need to generate each forward model at very high spectral resolution (R ≥ 250,000) and the need to reapply the detrending steps to each reconstructed data set. Calculated over the full 0.95–2.50 μm wavelength range, a SCARLET model takes 2–3 s to compute when using only a single CPU core. Processing such a model as per Section 3.2 for an observation sequence of 100–150 exposures further takes 2–5 s. The code is fully parallelizable to run simultaneously on many CPU cores on a cluster but can still take up to multiple weeks to run to full convergence. A typical SCARLET high-resolution retrieval has 27 parameters simultaneously being fitted, namely, 15 temperature points, two orbital parameters, a cloud top pressure, a scaling parameter, and the VMR of eight molecules.

As τ Boo b moves in time, the maximum of the CCF of each exposure should follow an orbital path as determined by the combination of the motion of Earth around the Sun, of the τ Boo system relative to the Sun, and of τ Boo b relative to its host star (Equation (3)). As such, a phase folding of the CCF for different combinations of V_sys and K_p should reveal the strongest signal at the true values. A comparison of the sum of the CCF in the known rest frame of the planet relative to erroneous velocity combinations can then give an estimate of the S/N of the atmospheric detection.

We begin by performing a blind search for τ Boo b's atmospheric signature in velocity space. The search is performed by setting up a velocity grid ranging from −300 to +300 km s⁻¹ in steps of 2 km s⁻¹ roughly corresponding to the velocity sampling of the SPIRou spectra on the detector pixel grid. For each of these velocity shifts, a planetary model is Doppler-shifted and then interpolated to the same wavelength points as the data. With both data and model on a common wavelength grid, the CCFs (Equation (6)) of each spectral order of each exposure are calculated and combined. We convert summed CCFs into S/N estimates by dividing them by the standard deviation of all values excluding those near the expected location.

The resulting K_p − V_sys S/N maps derived from CO and H₂O templates for the night of 2019 April 15 are shown in Figure 4. CO is clearly detected at K_p and V_sys values consistent with the known radial velocity of the star and the K_p inferred by previous CRIRES observations of τ Boo b (Brogi et al. 2012). On the contrary, even though the SPIRou data are more sensitive to H₂O than CO, no evidence for any water detection is seen. Here we see SPIRou's capability for detecting atmospheric signatures in even a single high-quality 5 hr observation sequence.

Zoom In Zoom Out Reset image size

A similar search performed on the other four nights shows CO being added constructively at the expected location, but, again, no sign of H₂O absorption is evident. The strength of the CO signal varies from being strong to marginal depending on the night, although this is expected owing to the varying observing conditions of each night. The data quality obtained from each night and how sensitive each is to an underlying planetary signal due to the varying integration time, water column density in Earth's atmosphere, and presence of systematics is investigated in more detail in Appendix B.

Having prominent features in the wavelength range covered by SPIRou (Figure 3), we also searched for CH₄, CO₂, HCN, TiO, C₂H₂, and NH₃ in all nights but found no statistically significant signal for any of these. This is, however, not surprising, as these molecules are not necessarily expected to be present in detectable levels in τ Boo b's atmosphere. We note that all conclusions regarding the detection or nondetection, as well as any abundance constraints, of molecular species is dependent on the accuracy of the line lists used (Section 3.1).

We calculate posterior distributions for the orbital and atmospheric parameters of τ Boo b using the SCARLET retrieval framework adapted to analyze high-resolution spectroscopic data. The retrieved abundances, orbital velocity components, line shape, and temperature structure for all spectral orders of all five SPIRou nights combined are shown in Figure 5 and tabulated in Table 2. We recover a slightly supersolar carbon monoxide abundance of log(CO) = $-{2.46}_{-0.29}^{+0.25}$ and find that models with very little water (3σ upper limit of log(H₂O) = −5.66, corresponding to 0.0072× solar) best match the data. We find no evidence for a significant abundance of CH₄, CO₂, HCN, TiO, C₂H₂, or NH₃ in τ Boo b's atmosphere, with upper limits for each of these molecules reported in Table 2. The dayside atmosphere of τ Boo b is consistent with being relatively cloud-free, with a cloud deck disfavored above the 250 mbar–1 bar pressure level. The cloud top pressure is partially degenerate with the CO abundance, as the muting of spectral lines by a higher-altitude cloud deck can be counteracted by an increase in line depth as a result of a higher CO abundance. The recovered temperature structure is noninverted, consistent with previous findings, and has a subadiabatic lapse rate of roughly dT/dlog₁₀(P) ∼ 300 K per pressure decade in the 10 mbar–1 bar regions most readily probed by the data. A typical CO line drawn from a random sampling of the scaled model thermal spectra in the Markov Chain Monte Carlo (MCMC) chain (all normalized to the same line continuum) shows a median depth of ∼70% relative to the continuum and no evidence of partial emission in the line core (Figure 5, top middle panel).

Zoom In Zoom Out Reset image size

Table 2.  τ Boo b Atmospheric Retrieval Results

Parameter	Value
VMR log(CO)	$-{2.46}_{-0.29}^{+0.25}$
VMR log(H2O)	<−5.66 (3σ upper limit)
VMR log(CH4)	<−3.78 (3σ upper limit)
VMR log(CO2)	<−3.99 (3σ upper limit)
VMR log(HCN)	<−5.37 (3σ upper limit)
VMR log(TiO)	<−7.54 (3σ upper limit)
VMR log (C2H2)	<−4.88 (3σ upper limit)
VMR log(NH3)	<−6.10 (3σ upper limit)
Cloud top pressure (Pc )	>0.26 bar (3σ lower limit)
Scaling parameter (a)	1.04 ± 0.03
Keplerian velocity (Kp )	109.2 ± 0.4 km s−1
Systemic velocity (Vsys)	−15.4 ± 0.2 km s−1
Derived parameters:
(C/H)gas	${5.85}_{-2.82}^{+4.44}$ × solar
(O/H)gas	${3.21}_{-1.56}^{+2.43}$ × solar
(C/O)gas	${1.00}_{-0.00}^{+0.01}$
(C/O)all	>0.60 (extreme limit)

Download table as:  ASCIITypeset image

We recover a planet velocity semiamplitude K_p = 109.2 ± 0.4 km s⁻¹ that is consistent with previous results using CRIRES (Brogi et al. 2012; Cabot et al. 2019; Watson et al. 2019). When combined with the stellar mass and stellar velocity semiamplitude (Table 1), the inferred K_p of τ Boo b corresponds to a mass of 6.24 ± 0.23 M_Jup and an orbital inclination of 435 ± 07. The velocity of the τ Boo system is constrained directly from the planetary signal to be V_sys = −15.4 ± 0.2 km s⁻¹. This is slightly offset from previous values between −16 and −16.4 km s⁻¹ derived from radial velocity measurements (Nidever et al. 2002; Donati et al. 2008; Brogi et al. 2012; Borsa et al. 2015) and from the predicted present-day systemic shift of roughly −17 km s⁻¹ expected owing to the motion of the wide-orbit M-dwarf companion τ Boo B (Justesen & Albrecht 2019). The global scaling parameter does not show significant degeneracy with recovered molecular abundances and is retrieved to be a = 1.04 ± 0.03, a value roughly consistent with unity, perhaps indicative of a slight underestimation of τ Boo b's atmospheric scale height. An insufficiently puffy modeled atmosphere may be caused by the assumed planet radius R_p = 1.15 R_Jup (Section 3.1) being slightly too small.

While a near-solar slightly enhanced abundance of CO on τ Boo b is in line with expectations from both theory and previous results, the very low upper limit determined for the water abundance is somewhat surprising. To better understand the retrieved results, we investigate how the best-fit model fares when injected in the raw data before any processing has been applied.

First, we verify the planetary nature of the CO signal by injecting the negative of the retrieved best-fit planetary model into the raw data at the inferred K_p and V_sys. If the signal is indeed from the planet and the model is a good fit, then the CO detection should be perfectly canceled and disappear. We find that this is indeed what happens. After injecting the negative planet signal and reapplying the detrending procedure (Section 2.4) and repeating the cross-correlation analysis (Section 4.1), there is indeed no CO signature left in the data (Figure 6). The near-perfect cancellation of the signal indicates that the atmospheric model recovered is at a spectroscopic contrast that well matches τ Boo b's underlying signature. We argue that this negative injection test is a verification procedure that should be widely applied in the literature to validate how well an atmospheric model truly represents the underlying planetary signal.

Zoom In Zoom Out Reset image size

We also investigate the lack of a clear water signature seen in the data. As H₂O has significant absorption throughout most of the SPIRou bandpass (Figure 3), any underlying signal should be readily found in the data with even more sensitivity than CO. To better understand the retrieved H₂O abundance upper limit, we explore how sensitive the data are to various injected water VMRs. Using the best-fit temperature structure, we inject simulated atmospheric models with different H₂O concentrations into the raw data and reperform the full data reduction and cross-correlation analysis for each case. The results of this water sensitivity test for the night of 2019 April 15 are shown in Figure 7, where the observed water K_p − V_sys S/N map is compared to different injected subsolar H₂O signals. Here we conclude that the absence of a water detection in the data is not due to a lack of sensitivity as even an injected VMR${}_{{{\rm{H}}}_{2}{\rm{O}}}=0.01\times$ solar can robustly be recovered from a single 5 hr SPIRou observation sequence.

Zoom In Zoom Out Reset image size

With the retrieved VMRs of all major carbon- and oxygen-bearing gaseous species in τ Boo b's atmosphere (Section 4.2), we can calculate corresponding elemental abundance ratios. Relative to protosolar abundances (C/H = 2.95 × 10⁻⁴, O/H = 5.37 × 10⁻⁴; Asplund et al. 2009), we find a carbon abundance of ${5.85}_{-2.82}^{+4.44}$ × solar, an oxygen abundance of ${3.21}_{-1.56}^{+2.43}$ × solar, and a C/O ratio of ${1.00}_{-0.00}^{+0.01}$ (Figure 8). The inferred C/H and O/H ratios of τ Boo b are supersolar and consistent with the corresponding values for Jupiter (Atreya et al. 2018; Li et al. 2020). With CH₄ as the next molecule after CO with the least strict upper limit abundance constraint, the distribution of inferred gas-phase C/O ratio values (Equation (1)) is slightly skewed toward carbon-rich (C/O > 1) scenarios. However, as refractory species can condense out and preferentially remove oxygen from the gas phase, the directly measured O/H ratio is realistically only a lower limit. We must therefore account for this to infer the true atmospheric C/O ratio of τ Boo b's envelope.

Zoom In Zoom Out Reset image size

4.4.1. Oxygen Sequestration

One plausible explanation as to why a true underlying water signal is not seen in the SPIRou data would be if the analysis is performed using an erroneous line list. While CO absorption produces a fewer number of deep, well-defined spectral lines, the H₂O spectrum is much more convoluted, composed of millions of weaker lines blended together. Variations across different line lists can cause significant biases in retrieved abundances (Brogi & Line 2019), or even cause strong detections to completely disappear such as in the case of HD 189733b (Brogi et al. 2016; Flowers et al. 2019). Most HDCCS water detections in the literature have utilized the HITEMP database based off of the BT2 line list (Barber et al. 2006). Also having shown some success is the more recent ExoMol updated opacities based off of the POKAZATEL line list (Polyansky et al. 2018). Tested on high-resolution data sets, these line lists have shown mostly consistent results in the wavelength ranges of 3.459–3.543 μm for HD 179949b (Webb et al. 2020), 0.92–2.45 μm for HD 209458b (Giacobbe et al. 2021), and 2.27–2.35 μm and 3.18–3.27 μm for HD 179949b and HD 189733b, respectively (Gandhi et al. 2020). We perform our analysis using both the POKAZATEL/ExoMol and BT2/HITEMP opacities but find near-identical results in both cases.

Another possible explanation could be that the retrieved TP profile, here constrained from the shape and contrast of CO lines (Figure 5), does not well represent the water lines. We investigate this by performing a free retrieval at the location of the planet and only including opacity contributions from water. Excluding other molecules allows the temperature structure to be fit entirely from H₂O lines, as opposed to being driven by the CO signal. Such a search, however, did not reveal any hint of appreciable water absorption or emission.

Alternatively, as the cores of CO and H₂O lines form for the most part at different pressure levels of the atmosphere, differences in low- and high-altitude winds may cause the H₂O signal to be blurred below a detectable level (Brogi & Line 2019). Accounting for convoluted effects from contrasts in low- and high-altitude winds would require the use of more sophisticated high-resolution 3D models (e.g., Flowers et al. 2019; Beltz et al. 2021; Harada et al. 2021). Recent work has shown that such 3D atmospheric models that take into account broadening effects from winds and rotation fare significantly better than 1D models at recovering molecular signatures in high-resolution dayside emission data sets (Beltz et al. 2021).

In comparison to CO molecules that only have appreciable opacity near the 1.6 and 2.3 μm spectral band heads, H₂O vapor has transition lines throughout the entire SPIRou wavelength range (Figure 3). However, the accuracy of the line list, the level of contamination from residual tellurics, and the presence of instrumental systematics can all vary greatly across different spectral orders. It is therefore possible that positive correlations from water in some orders are drowned out by negative contributions from more problematic noisier orders. We investigate this by analyzing each of SPIRou's photometric bands (Y, J, H, K) individually. For this we fix the temperature structure points at their median retrieved values and perform retrievals for each data subset at the location of the signal. H₂O posterior distributions for each individual band of all five SPIRou nights combined are shown in Figure 9. No photometric band shows evidence of a water signal, with each providing an independent upper limit on the water abundance. The H and K bands, where the SPIRou data have the highest S/N and the planet-to-star flux contrast is more favorable, provide the most strict upper limits on the H₂O VMR.

Zoom In Zoom Out Reset image size

It is also possible that the detrending algorithm used on the data fares well for uncovering CO lines but is not optimal for extracting a water signal. As some parts of the SPIRou bandpass are more plagued by tellurics than others, spectral orders falling in these heavily contaminated wavelength regions would potentially require a higher level of detrending or masking to uncover weaker absorption lines from H₂O on τ Boo b. While it is common practice to optimize the number of PCA iterations or the telluric mask cutoff percentage of each order individually to maximize the retrieved significance of either the real signal or an injected signal, doing so can result in noise patterns being optimized into false positives (Cabot et al. 2019; Zhang et al. 2020b). This is particularly true given that SPIRou has many spectral orders, which would produce a high number of free parameters to be optimized. As the data already achieve a high degree of sensitivity to the water abundance under the general reduction of Section 2.4, we prefer to avoid this practice and treat all orders equally to avoid introducing any biases.

The potential partial inaccuracies of the line lists used, the neglect of 3D wind and rotation broadening effects, and the suboptimal but unbiased generalized practice of treating all orders identically during the detrending process may all have some effect as to the exact value of the upper limit of the H₂O abundance. However, even considering these caveats, it is clear from Figure 7 that even a modestly subsolar abundance of water is hard to reconcile with the SPIRou data.

Finally, at high temperatures it is possible for H₂O molecules to thermally dissociate. This would act to reduce the strength of observed water spectral features. Meanwhile, CO would remain mostly unaffected given its significantly stronger molecular bond (Lodders & Fegley 2002). If the temperature is high enough on τ Boo b's dayside, it is therefore plausible that H₂O would be partially dissociated while CO molecules remain intact, thus naturally explaining the retrieved supersolar CO and significantly subsolar water abundances from the SPIRou data.

Water dissociation is efficient in high-temperature, low-pressure environments such as the upper atmospheres of ultrahot Jupiters with strong thermal inversions. Specifically, H₂O starts to break down above the ∼100 mbar level at temperatures greater than about 2200 K (Parmentier et al. 2018). However, τ Boo b shows no sign of a strong temperature inversion and most likely does not reach high-enough temperatures in its upper atmosphere for water dissociation to play a major role. It is thus improbable that H₂O dissociation is prominent enough on τ Boo b to explain the reported low water abundance.

In a static protoplanetary disk, the different condensation temperatures of the main carbon- and oxygen-bearing ices introduce distinct gradients in the metallicity and C/O of the gas and solids (Öberg et al. 2011). Primarily, the disk is composed (∼99%) of molecular hydrogen and helium that remains mostly in a gaseous state. The heavier elements (C, O, N, etc.) composing the remaining ∼1%, however, can be in either gas or solid form, depending on the local temperature. Near the star it is warm enough that nearly all of the material is in vapor form for part of the gaseous phase. At larger distances from the star, where it is colder, different molecules condense out of the gas and add to the dust grains in the disk. This causes the metallicity of the gas to drop with orbital separation as more and more metals freeze out. Since oxygen-rich molecules such as H₂O and CO₂ condense out at higher temperatures than CO, the C/O ratio decreases in solids and increases in the gas farther out in the disk. The result is then a protoplanetary disk with oxygen-rich, low-C/O solids and metal-poor, superstellar C/O gas beyond the water ice line (Madhusudhan et al. 2014a).

Within this simple picture of a static disk model, different regimes of gas- or planetesimal-dominated accretion occurring at distinct orbital separations can be inferred from the present-day atmospheric compositions of giant planets. In a gas-dominated accretion regime (no substantial planetesimal accretion), the planet's present-day envelope may reflect the gas-phase C/O of the disk at the location where the runaway growth occurred (Öberg et al. 2011). Such a scenario would most likely occur if the planet migrated to its current orbit without accreting a significant amount of solids (disk-free migration). Assuming an isolated core, gas-dominated accretion beyond the water ice line preferably enables the formation of low-metallicity, high-C/O planets.

Alternatively, in a regime dominated by planetesimal accretion, the results will be an atmosphere that is heavily polluted by metals from captured solids. Planetesimal accretion can occur during the envelope contraction phase of runaway growth, during the clearing of the feeding zone, and while the planet migrates to becoming a hot Jupiter while interacting with the still-present protoplanetary disk (Type II migration). Solid-dominated accretion typically results in planets that are metal-rich but have low C/O ratios (Madhusudhan et al. 2014a, 2017; Cridland et al. 2016, 2019, 2020; Mordasini et al. 2016; Ali-Dib 2017; Espinoza et al. 2017; Nowak et al. 2020). Core erosion can also contribute to increase the metallic content of the envelope with stellar or substellar C/O material, depending on where in the disk the core formed (e.g., Öberg & Wordsworth 2019). However, by building a planet from metal-poor, superstellar C/O gas with metal enrichment from substellar to stellar C/O solids either from the core of from planetesimals, there is no outcome that allows for the formation of high-metallicity, high-C/O planets like τ Boo b regardless of formation location and gas/solid accretion fractions (Madhusudhan et al. 2014a, 2017).

One possible mechanism that naturally explains a carbon-rich, high-C/O gas giant envelope is by enrichment of the gas in the disk due to pebble drift across condensation fronts (Öberg & Bergin 2016; Booth et al. 2017). Indeed, observations show that protoplanetary disks contain large reservoirs of millimeter- to centimeter-sized pebbles (e.g., Rodmann et al. 2006). Such grains are large enough to decouple from the gas and feel a headwind that causes them to drift radially inward. The migration of icy grains to the warmer inner regions of the disk can then cause the transported frozen volatiles to evaporate back to the gas phase and enhance the local chemical abundance near the condensation front (Ali-Dib et al. 2014). For example, pebbles drifting inward from the outer disk will begin to sublimate once they cross the CO snowline, which will cause a local metal enrichment of the gas phase. The deposited excess CO molecules can then diffuse and advect throughout the disk, increasing the metallicity of the gas phase inward of the CO snowline while maintaining a C/O ∼ 1 (Öberg & Bergin 2016). Evidence of pebble drift in the form of excess C/H within the CO snowline was recently observed in the HD 163296 disk (Zhang et al. 2020a). Gas-dominated accretion in a pebble-drift-induced, CO-enriched environment provides a natural mechanism to form metal-rich, high-C/O planets like τ Boo b, without the need for substantial enrichment from planetesimal accretion or core dredging. In the scenario where the C/O ratio is overestimated (but still supersolar) owing to significant oxygen condensed out of the gas phase, pebble drift past the CO or CO₂ snowlines would still be favored, although a more significant enrichment from low-C/O, high-O/H planetesimals may have occurred.

We note that planet formation is likely to be much more complex than the simple picture depicted here, and a number of processes can act to alter our conclusions. For one, protoplanetary disks can undergo significant evolution, changing the chemical composition and position of snowlines over time (Cridland et al. 2016, 2019, 2020; Eistrup et al. 2016, 2018). It is also possible that the contribution of accreted planetesimals to the net envelope metallic enrichment is overestimated if captured planetesimals reach the core instead of being dissolved and mixed into the atmosphere (Pinhas et al. 2016). Similarly, the degree to which the core interacts and exchanges material with the envelope is uncertain. The assumed oxygen-rich, low-C/O composition of the solids past the water ice line may also be erroneous if more carbon than expected is held in refractory organic material (e.g., Lodders 2004). All interpretations also rest on the fundamental assumption that the probed atmospheric abundances in hot Jupiters are representative of the bulk composition of the planet, which may not necessarily be the case (Leconte & Chabrier 2012; Vazan et al. 2018). Finally, τ Boo b may have formed in a truncated disk owing to the system's highly eccentric wide-orbit companion, which would likely further complicate this simple picture (Justesen & Albrecht 2019).

Another possible pathway to form giant planets is via gravitational instability (e.g., Boss 1997, 2001). This is a relatively short timescale process that results in an envelope initially representative of the local bulk composition of the disk. With subsequent enrichment from solids, planets formed via gravitational instability are expected to have stellar or superstellar metallicities and C/O ratios that are stellar or substellar (Madhusudhan et al. 2014a). It is unclear how a gravitational instability formation scenario could produce a τ Boo b−like high-C/H, high-C/O planet.

So how does τ Boo b and its Jupiter-like C/H and elevated C/O fit in relative to other hot Jupiters? While the lack of water observed in the SPIRou data is unexpected, subsolar H₂O abundances on exoplanets are not uncommon. Systematic analyses of a sample of hot Jupiters observed in transit spectroscopy have shown a favorability for water VMRs that are depleted relative to expectations from solar (Barstow et al. 2017; Pinhas et al. 2019; Welbanks et al. 2019). A comparison of the retrieved H₂O VMR upper limit on τ Boo b to water measurements on hot Jupiters obtained from transmission spectroscopy with HST/WFC3 is shown in Figure 10. While most reported H₂O abundances are consistent with being subsolar, the stringent upper limit derived from the SPIRou data is lower than any of these. However, the low water abundance is not necessarily evidence of a low O/H if, instead of in H₂O, the oxygen budget of τ Boo b is held in other molecules such as CO. This is why constraining the abundance of multiple species is crucial for understanding the chemical budget of exoplanetary atmospheres.

Zoom In Zoom Out Reset image size

With a measure of its CO abundance, τ Boo b joins HD 189733b and HD 209458b in having precise constraints on both H₂O and CO obtained from high-resolution spectroscopy observations (Brogi & Line 2019; Gandhi et al. 2019). With only these very few robust and precise measurements that currently exist, current trends of CO abundances on hot Jupiters as determined from high-resolution spectroscopy favor slightly enhanced Jupiter-like C/H ratios (Figure 11). As the respective main carriers of O and C under equilibrium chemistry, this allows for the C/O ratios of these planets to be estimated. In all three cases, the CO abundances are significantly higher than the H₂O, and the CO/(H₂O+CO) ratio, a proxy for the C/O in the absence of condensates or other carbon- and oxygen-bearing species, can constrained to be ∼1. In the case of τ Boo b, gas-phase C/O ratios significantly deviating from unity are disfavored owing to the upper limit constraints derived on other potential carbon- or oxygen-bearing molecules CH₄, CO₂, HCN, C₂H₂, or TiO.

Zoom In Zoom Out Reset image size