1 Introduction↩︎

Direct imaging is a powerful technique for detection and characterization of exoplanetary systems, providing access to wide orbital separations and enabling astrometric and spectroscopic measurements to constrain planetary bulk, orbital, and compositional properties. This ensemble of information is useful for illuminating planetary demographics [1] as well as distinguishing between possible formation pathways [2]. Early simulations of JWST operating in a standard coronagraphic mode with NIRCam [3] set expectations for increases in direct imaging sensitivity to very low masses (\(\sim\) Saturn mass) at very wide separations (\(>\)​50 AU). This increase in sensitivity has since been confirmed by cycle 1 survey programs [4], although candidate objects remain unconfirmed. Alternatively, NIRSpec’s design was flexible for use in fixed-slit or IFU-based modes with a low-resolution prism (R = 30-330) or higher-resolution gratings (R = 500-1340 or R = 1320-3600) from 0.6-5.3 \(\mu\)m [5], [6]. While NIRSpec was originally intended to image faint galaxies, others quickly realized that NIRSpec was capable of characterizing exoplanets with transmission and eclipse spectroscopy [7].

In this paper, we use JWST NIRSpec in an unconventional observing mode, using fixed-slit spectroscopy to image and characterize 51 Eridani b at high contrast. We build on the methods laid out in [8], [9] which demonstrated spectroscopic detection and atmospheric characterization of HD 19467 B at moderate contrast with the NIRSpec IFU. This method relies on using moderate spectral resolution observations combined with joint modeling techniques to separate the stellar and planetary contributions to the total observed signal. When compared to low resolution (R\(\sim\)​100) spectroscopy, which can detect broadband molecular absorption features from many groups of transitions at proximal wavelengths/energies, moderate resolution (R\(\gtrsim\)​1000) spectroscopy can resolve individual molecular bands from specific rotational and vibrational transitions. When combined with high-pass filtering to remove the stellar and planetary continuum, these molecular features remain in the continuum-subtracted data as oscillations in the residuals with unique spectral fingerprints. Cross-correlation analysis between model templates with known molecular features and these residuals can reveal the existence of the planetary contribution to the total flux. This technique is capable of detecting planets at high contrast even in the absence of a coronagraph to block the starlight. This is possible because spectral correlations can be used to draw planetary signal out of the spatially variable “speckle-noise" from scattered starlight, which is the dominant noise source for standard broadband coronagraphic imaging. The significant heritage of applying this technique on ground based observatories, including with KPIC [10], [11], OSIRIS [12], and CRIRES [13], led to the rapid application of this powerful approach to a space-based platform.

The observational strategy using the fixed slits was originally chosen to address concerns that the very bright (V=5.2 mag) [14] host star would saturate the detector in full frame exposures. This would lead to unrecoverable pixels at the position of the planet as “charge bleeding" effects in the detector cause photoelectrons to spill over into neighboring pixels. By observing with the slit, not only would the starlight from the central PSF core of the host star be mitigated, but also it would enable subarray readout mode with faster cadences. This strategic observing choice has other consequences which we discuss in greater detail later, including additional exposures/pointings required for the”speckle-side" dataset and loss of representative spatial sampling of the “noise annulus" for estimation of the covariance of the errorbars.

The target we observe, 51 Eridani b, has been studied previously from the ground using GPI, SPHERE, and Keck/NIRC2. [15]–[19] The discovery paper [15] reported a detection of methane absorption around the J band (1.1-1.4 \(\mu\)m) which made it unique among directly imaged objects around the same age and temperature, hinting at the presence of disequilibrium chemistry. While obtaining a dynamical mass measurement of 51 Eridani b from a combination of astrometry and Gaia-Hipparcos acceleration is not yet possible (De Rosa, priv. comm.), evolutionary models can estimate the planetary mass from its age and luminosity. These models indicate that 51 Eridani b is between 2-12 M\(_\textrm{Jup}\), which is one of the lowest mass directly imaged objects known. This makes 51 Eridani b “a bridge between wider-orbit, hotter, more massive planets and planets at Jupiter-like scales" [15].

Previous atmospheric analyses using forward modeling found that a “partly-cloudy" atmosphere [20] was needed to explain the spectral energy distribution. However, retrieval analyses [18], [19] preferred clear models. The authors attribute this to the flexibility of the retrieval framework to”mimic" cloud-like effects in the thermal structure of atmosphere. However, all of these previous analyses had assumed equilibrium chemistry or free chemistry for the retrieval. Our previous work uniquely included modeling of effects from disequilibrium chemistry and clouds together [21], in particular the disequilibrium between carbon monoxide and methane. These effects when combined can resolve the large discrepancy between the models and data for the photometric point from Keck at 4.5 micron in the Ms band. The possibility to directly detect carbon monoxide absorption features motivates the 3-5 \(\mu\)m spectroscopy with JWST / NIRSpec we report in this paper. The partly-cloudy disequilibrium models are consistent with the reported non-detection of 51 Eridani b with cycle 1 JWST NIRCam observations in filters F430M and F460M [22].

The paper is organized as follows. In Section 2 we describe the observations that will be analyzed throughout the paper, in particular the NIRSpec spectra and accompanying ground-based data from GPI and Keck. In Section 3 we describe the methods used to reduce the new JWST / NIRSpec data including the JWST pipeline and the custom modeling package BREADS [23]. Section 4 contains the bulk of the analysis of the NIRSpec data products, including computations of the cross-correlation functions, “leave one molecule out" tests to search for the presence of individual molecules such as carbon monoxide and methane, as well as a comprehensive spectral analysis of the ground and space data together in a joint likelihood framework. Section 5 concludes the paper with a discussion that summarizes our main results and their future implications.

2 Description of Observations↩︎

The new observations that we present in this paper were executed during JWST Cycle 2 as a part of program 3522. The instrument is configured to observe both a “planet-side" and”speckle-side" pointing for both of the 200 mas wide fixed slits S200A1 and S200A2. Since the host star is too bright for target acquisition, we use the target Gaia DR3 3205095194041180160 for wide angle target acquisition and offset from that position. The filter wheel is set to F290LP, which covers \(\lambda \in [2.87, 5.27]\) \(\mu\)m. The disperser of choice is the high resolution (R \(\sim\) 2700) grating G395H. This enables a complete spectral reconstruction across the gap between the two NIRSpec detectors (NRS1 and NRS2), which subtends \(\lambda \in [3.69, 3.79]\) \(\mu\)m in the “A1" slit and \(\lambda \in [3.81, 3.92]\) \(\mu\)m in the”A2" slit. We name these sequences of observations in shorthand as “a1p,"”a1s," “a2p," and”a2s" to concisely refer to specific combinations of slits (a1 and a2) and pointings (planet-side and speckle-side).

The detectors are configured to readout in the fastest possible configuration NRSRAPID, using the subarray modes SUBS200A1 and SUBS200A2 to avoid saturating the detector due the extremely bright nature of the target. The observation sequences are dithered with the maximum possible number of dithers, five primary positions and four additional subpixel dithers at each primary position. This dithering improves the spatial sampling of the observations, with a total of 20 dithers per sequence. Each exposure is programmed to use four groups per integration and 40 integrations per exposure for a total exposure time of \(\sim\)​3.43 hours across all configurations.

In addition to the newly analyzed slit spectra, we also analyze previously published spectra and photometry of 51 Eridani b in a new joint analysis. The previously published results include the original J and H band GPI spectra from the discovery paper [15], additional GPI K1 and K2 spectra as well as Keck/NIRC2 Lp and Ms photometry that was published in [20], and JWST/NIRCam F410M, F430M, and F460M photometry that was published in [22].

3 Data Reduction Methodology↩︎

3.1 JWST Pipeline↩︎

The data from program 3522 are reduced using a combination of the JWST pipeline [24] and BREADS (the Broad Repository for Exoplanet Analysis, Discovery, and Spectroscopy) [8], [23], with minor modifications and customizations specific to this program. First, the raw data or uncal files are downloaded using the helper tool jwst‘_mast‘_query. Next, we use version 1.17.1 of the JWST pipeline to run stages 1 and 2. Stage 1 of the pipeline includes the processing steps group‘_scale, dq‘_init, saturation, superbias, refpix, linearity, dark‘_current, charge‘_migration, jump, clean‘_flicker‘_noise, ramp‘_fitting, and gain‘_scale. All of these steps are run with the default parameters with the exception of saturation, where we pass n‘_pix‘_grow‘_sat = 0 instead of the default of 1. This has the effect of limiting the saturation checking from expanding to neighboring pixels. At the end of stage 1 the pipeline saves the 2 dimensional rate files in units of data numbers per second or DN/s.

There is a consensus in the transit spectroscopy community about the importance of an additional data processing step which we will refer to as the stage 1.5 correction since it occurs between stages 1 and 2. The idea is to correct correlated read noise or “1/f" noise by masking the spectroscopic trace in the rate files and performing a column-independent linear regression on the dark or unilluminated pixels to flatten or zero out their value [25]–[38]. We have experimented with including or not including this additional customization and found that it did not provide an improvement in the quality of our results. We have therefore opted to forgo this step in our final analysis, electing to run the default nsclean step [39] instead. This may be explained in part by the continuously evolving nature of the JWST pipeline software. For reference, our results are using the 1.17.1 version of the pipeline which was just released in January 2025, after the publication of most of the works advocating for stage 1.5.

Next is stage 2 of the JWST pipeline, which includes the data processing steps: assing‘_wcs, nsclean, background, extract‘_2d, srctype, wavecorr, flat‘_field, pathloss, photom, pixel‘_replace, resample‘_spec, and extract‘_1d. All of these steps are run with the default parameters with the exception of extract‘_1d which is skipped. The stage ends by saving the two dimensional cal files in units of MJy, which serve as the starting point for the continued custom post-processing we perform using BREADS.

3.2 Baking BREADS↩︎

The original demonstration of using BREADS to perform high-contrast spectroscopy with JWST data is shown in great detail in [8], where HD 19467 B was imaged and characterized using the NIRSpec IFU. Our analysis follows very similar principles with a few minor modifications that are specific to our observations, which use the NIRSpec fixed slits instead of the IFU. The major concept behind the BREADS framework is to model the data directly in the detector space, where we refer to the four-dimensional object \((\alpha, \delta, \lambda, F)\) as the point-cloud, where \(\alpha\) is the right ascension, \(\delta\) is the declination, \(\lambda\) is the wavelength, and \(F\) is the stage 2 calibrated flux in units of MJy. Each of these objects is implicitly a function of the detector pixel coordinates \(x\) and \(y\). The pixel coordinates are not to be confused with the IFU-aligned coordinate system IFU \(X\) and IFU \(Y\), which is simply the sky coordinates \((\alpha, \delta)\) rotated by the parameter east‘_2V2‘_deg, such that the orientation of the slit is vertical along the IFU \(Y\) coordinate axis. The parameter east‘_2V2‘_deg is derived by taking the difference between the two JWST header objects V3I‘_YANG \(-\) ROLL‘_REF.

By modeling the data in the detector space with a combination of a spline-based continuum and empirical model for high-frequency components of the stellar spectra, we can accurately subtract the stellar contribution to the total flux at the cost of losing the continuum flux from the planet. The continuum-subtracted spectra are then cross-correlated with a template of the planetary spectrum, which contains unique signatures due to molecules such as CH\(_4\), CO\(_2\), and CO. At an effective stellar temperature around 7500 K [40], these molecules do not exist in the atmosphere of 51 Eridani A and reveal the existence of the planetary contribution to the flux.

One of the major free parameters in this analysis is the choice of node spacing for the spline continuum model. Due to the complexity of optimizing the reduction pipeline over this hyperparameter, we opt to only investigate uniform node spacings. This means the node spacing can be characterized by a single integer \(N\) which is the total number of nodes needed to cover the entire wavelength range \(\lambda \in [2.86, 5.13]\). We investigated a range of values for \(N\) from 88 to 152 in increments of 2 and found that our results are somewhat insensitive to this choice, indicating the detections are robust to the choice of hyperparameters. However, the best results were obtained for N = 136, and therefore we chose this value for our final analysis. This value of N results in 58 nodes covering the wavelength ranges in NRS1 and 81 nodes covering the wavelengths in NRS2, corresponding to a wavelength spacing of \(\Delta\lambda = 0.0167\) \(\mu\)m\(/\)node. For comparison, the typical length scale of a molecular feature is 0.00278 \(\pm\) 0.00199 \(\mu\)m which is around six times narrower than the node spacing. A summary plot of the detection significance as a function of this reduction hyperparameter is shown in the left panel of Figure (1).

To describe the specific analysis routines in greater detail, the first step in the BREADS chain is pre-processing the data to improve the bad pixel identification. Each row of the data is high pass filtered using a median filter with a window size of 50 pixels. Then, we compute the median absolute deviation and flag all pixels greater than 50 times that value as bad pixels. The next step is to compute the sky coordinate arrays over the detector pixel indices using the calfile.slits[0].meta.wcs objects populated by the JWST pipeline. This step is one of the customizations needed for fixed slit pipeline. After this pre-processing, the data are selectively masked with NaNs to block out large spectral features due to hydrogen absorption in the stellar spectrum. We find the hydrogen masking improves the final quality of the detection as it was difficult to precisely estimate the exact depth of the hydrogen absorption features from the empirical stellar spectra. We speculate that this is because of the different illumination profiles of slit between the speckle and planet pointings. The exact masking we apply corresponds to central wavelengths \(\lambda \in \{3.0395, 3.297, 3.741, 3.8125, 4.0525, 4.654\}\) \(\mu\)m with a half-width of \(\delta\lambda\) = 0.005 \(\mu\)m. These wavelengths correspond to atomic transitions in the hydrogen atom between the states (\(n_1\), \(n_2\)) corresponding to (5, 10), (5, 9), (5, 8), (6, 16), (4, 5), and (5, 7) respectively, which is mostly the Pfund series of absorption lines.

The next steps are to compute the continuum-normalized stellar spectra, fit and subtract the starlight continuum spline model, and interpolate the residuals onto a regular wavelength grid. Each slit/detector pairing (e.g. a1n1, a1n2, a2n1, and a2n2) is analyzed independently. One unique modification specific to the fixed slit analysis pipeline is that we compute the stellar spectra only from the speckle-side pointing dataset. If a dataset is from the planet-side sequence, we replace the typical dataobj.star‘_func with the equivalent stellar spectra from a speckle-side sequence before computing the starlight subtraction. Another modification specific to the fixed slit pipeline is that we prepare two different kinds of “combined" data objects containing all of the data from the entire observation sequence. One class of combined data objects contains the star-subtracted data from the entire sequence of observations in the planet-side and speckle-side pointings individually. There are a total of eight combined data objects in this set, two for each slit, detector, and pointing. These are used to extract the spectrum of the planet and the spectrum of the noise in the speckle pointing since the starlight continuum has been removed. The other combined data objects contain the continuum-preserving data from both the planet and speckle pointings combined together. These combined data objects are used to apply a coordinate correction / centroiding operation to more precisely calibrate the point cloud coordinates and rely on the starlight continuum remaining to precisely fit for that centroid.

The centroiding / coordinate correction uses the continuum-preserving, planet\(+\)speckle combined sequences and runs the fitpsf module to fit a synthetic WebbPSF to the point cloud as a function of \(\lambda\) and find the best-fit coordinate correction \((\Delta\alpha, \Delta\delta)\). The values range from -0.02 to -0.08 arcsec in both right ascension and declination with a small amount of curvature along the spectral axis. The empirical best-fit coordinate correction as a function of wavelength is regularized by fitting a second-order polynomial to smooth out systematic errors caused by spatial undersampling. The coordinate correction as a function of wavelength is shown in the right plot of Figure (1). It is important to recognize that this coordinate correction is sub-pixel and so represents a relatively small calibration correction to the existing JWST WCS system, but necessary in order to precisely align the expected position of the planet in the center of the sub-pixel dithering pattern, which has offsets of around \(\pm\)​0.025 arcsec in both directions. The planet’s positional uncertainty based on orbit fitting is around .0002 arcsec [41], which is an order of magnitude more precise than the WCS coordinate calibration, since the typical JWST pointing accuracy is a “few mas" (i.e. 0.002-0.003") after target acquisition [42].

After the coordinate correction is applied, the final step in the BREADS chain is to perform a spectral extraction. The fitpsf module is used to find the best fit flux (as a function of position and wavelength) for the starlight-subtracted residual data, using a synthetic WebbPSF. The position locations to perform the fit are an input parameter and it is important to have previously performed the regular wavelength interpolation so that fitpsf can treat each columnar slice of the combined data objects as an independent wavelength channel. This regular wavelength grid interpolation is known to result in a 2 pixel wide covariance in the final extracted errors, which was demonstrated in the appendix of [8]. A major disadvantage of using the slit compared with the IFU is that there are substantially less spatial locations with similar flux that are sampled to reliably estimate the empirical error covariance. The final result of the BREADS chain are the spectral extraction products which are the continuum-subtracted residual fluxes as a function of position and wavelength and their corresponding uncertainties.

4 Spectral Analysis and Model fitting↩︎

4.1 Cross Correlation Analysis↩︎

In order to analyze the extracted spectra, we rely on a template-based cross-correlation technique. We employ substellar atmosphere models from Sonora Elf Owl [43], [44]. We use a combination of the species software package [45] as a backend and the breads.atm‘_utils.miniRGI module for convenience to load and gaussian kernel broaden the model spectra at a resolution of R = 2700. Because the extracted spectra have been effectively high-pass filtered by the spline-based continuum subtraction, it is important that the models are equivalently high-pass filtered with the same spline node configuration. It is possible to accomplish this with the BREADS utility function breads.utils.filter‘_spec‘_with‘_spline. Using a model template with the following parameters: \(T_\mathrm{eff} = 900\) K, \(\log \big( \frac{g}{\mathrm{cm}/\mathrm{s}^2} \big) =\) 3.5, \([\mathrm{M}/\mathrm{H}] = 0.5\), \(\textrm{C}/\textrm{O} = 0.458\), and \(\log \big(\frac{K_\mathrm{zz}}{\mathrm{cm}^2\mathrm{s}}\big) = 2\), we compute the cross correlation function with the following equation: \[\mathrm{CCF}\;(\mathrm{IFU}\;Y, \mathrm{RV}) = \frac{\sum_\lambda \frac{d_\lambda m_\lambda}{e_\lambda^2}}{\sqrt{\sum_\lambda \frac{m_\lambda^2}{e_\lambda^2}}}\] where \(d\) represents the data, \(e\) is the error on the data, and \(m\) represents the model. The data and its error as a function of IFU \(Y\) are directly taken from the spectral extraction products, while the model \(m\) is implicitly a function of the model parameters and the planetary radial velocity. The radial velocity has an effect which is equivalent to stretching the effective wavelengths of the model under the following remapping: \(\lambda \rightarrow (1 + \frac{\mathrm{RV}}{c})\lambda\), which we accomplish by resampling / interpolating the model with the module np.interp. The CCF equation essentially represents fitting the template model to the data in a linear least squares fashion, and taking the ratio of the flux-rescaling constant to the uncertainty on the fit over the IFU \(Y\) - RV phase space (see also Equation (5) in [46]). This means that it is not strictly a cross-correlation function according to the standard mathematical definition, since the standard CCF relies on an additive shift for the dependent variable but instead the RV parameter here multiplicatively stretches the model function in the wavelength domain. Regardless of this nuance, the CCF as we have defined it is a very useful diagnostic to elucidate the statistical correlations between the model and data and is often referred to as a cross-correlation function in the literature. These CCFs are computed for each dataset individually, but can then be

combined using a inverse-variance weighted-mean of the best fit flux and uncertainty. In order to ensure the accuracy of the noise model for the CCF, we rescale by dividing by the standard deviation of the “wings" of the CCF, i.e. where the \(| \mathrm{RV} | > 500\) km/s.

The CCFs for both the joint slit a1 + a2 datasets for both the planet-side and speckle-side pointings are shown in Figure (2). This figure demonstrates the clear detection of 51 Eri b in the planet-side dataset at a level of 4.8\(\sigma\), localized in the IFU \(Y\)-RV phase space at its expected position. If the signal was caused by systematic residuals from the subtraction of the diffracted starlight, we would expect the peak of the CCF to be aligned with the peak of the starlight intensity. We verify that the two are not co-located in the upper right panel of Figure (2), which provides additional evidence that the signal is of planetary origin.

4.2 Molecular Detection Experiments↩︎

In addition to the spatial misalignment between the peak of the stellar flux and the peak of the planetary spectrum CCF, another line of evidence for the signal’s planetary nature is to show that it originates from molecules that are expected in the planetary spectrum. At the effective temperature \(\sim\)​750 K in the NIRSpec wavelengths, these are primarily methane (CH\(_4\)) and carbon monoxide (CO), although there are smaller contributions from water (H\(_2\)O) and carbon dioxide (CO\(_2\)). Other molecules which are possibly relevant at these wavelengths include hydrogen sulfide (H\(_2\)S) and phosphine (PH\(_3\)) although they are expected to be subdominant absorbers due to their lower relative abundance.

To demonstrate the effects these molecules have on the significance of the CCF detection, we perturb the spectral model of the planetary atmosphere using PICASO [47]. We start by loading the Sonora Elf Owl model atmosphere with parameters \(T_\mathrm{eff} = 900\) K, \(\log \big( \frac{g}{\mathrm{cm}/\mathrm{s}^2} \big) =\) 3.5, \([\mathrm{M}/\mathrm{H}] = 0.5\), \(C / O = 0.458\), and \(\log \big(\frac{K_\mathrm{zz}}{\mathrm{cm}^2\mathrm{s}}\big) = 2\). This is accomplished in PICASO using jdi.xr.load_dataset and jdi.input_xarray. It is important that we use the 2025 updated opacity tables all_opacities_0.6_6_R60000.db in order to accurately reproduce the published Elf Owl spectra with relative errors on the order of \(\sim 10^{-5}\). These opacity tables have not been published yet but should be at some point and in the mean time can be requested by contacting the developers (Batalha, priv. comm.). For our simple test, we simply set the relative abundances for an individual molecule to zero at all pressure layers to perform our “leave one molecule out" comparison. These spectra and additional visualizations such as the molecular contribution function are shown in Figure (3), while the CCFs computed using these spectra are shown in Figure (4) alongside two additional tests which restrict the wavelengths to the regions of CH4 and CO dominance.

In the upper panel of Figure (4) one can see how leaving out any individual molecule from the spectra produces a corresponding decrease in the resulting peak of the cross correlation function, which indicates that the features from that molecule are relevant and contribute to the detection of the planetary signal. So the detection of the planetary signal is not just from a single molecule, but from a combination of multiple relevant molecules, all of which could not exist in a stellar atmosphere. In this test, we measure the impact of an individual molecule by computing the difference in significance between the CCF detection using the full model compared to a model with that molecule missing: \(\Delta \sigma = \sigma_\textrm{leave one molecule out} - \sigma_\textrm{full model}\). We find: \(\Delta\sigma_\textrm{CH_4} = -3.35\), \(\Delta\sigma_\textrm{CO} = -1.81\), \(\Delta\sigma_\textrm{H_2O} = -0.32\), and \(\Delta\sigma_\textrm{CO_2} = -0.58\). This test shows how the features from methane and carbon monoxide are the most important in detecting the planetary signal, and for the first time, a direct detection of these two molecules in chemical disequilibrium in 51 Eridani b. The test also shows what we will call a “consistent non-detection" for water and carbon dioxide, in that there is simply not enough evidence at the current level of sensitivity to claim a detection of these molecules on their own, although having them in the model does provide an improved detection.

To make the detection significance of the individual molecules more precise, we perform additional restricted wavelength tests in the two lower panels of Figure (4). The middle panel shows the CCF computed in the region \(3.0 < \lambda < 4.2\) which is the range where CH4 is the dominant absorber, while the bottom panel shows the CCF in the region \(4.4 < \lambda < 5.0\) where CO is the dominant absorber. In these tests one can see how removing the dominant molecule from the spectra greatly reduces the significance of the detection, and removing the irrelevant molecule does not alter the detection significantly. Furthermore, since these molecules are the primary absorber at these wavelengths, we can interpret the detection significance in this restricted wavelength test as the detection signficance of the molecules individually, implying \(\sigma_\mathrm{CH4} = 2.1\) and \(\sigma_\mathrm{CO} = 4.2\).

Perhaps even more interesting than just how much each molecule contributes to the total detection significance is how that is related to the total fractional bandpass where that molecule’s features dominate the absorption spectra. For example, by counting the total number of wavelength channels where a molecule’s contribution function (Figure (3) upper panel) is at a higher altitude than all of the other molecules (or equivalently, lower pressure layer) we can estimate the total fractional bandpass per molecule where that molecule’s features are dominating the spectral shape. This is simply the ratio of the number of channels where its contribution function is extremal compared to the total number of channels. We find the following values for the fractional bandpass domination: \(f_\textrm{CH_4} = 47.8\%\), \(f_\textrm{CO} = 20.9\%\), \(f_\textrm{H_2O} = 22.5\%\), and \(f_\textrm{CO_2} = 8.8\%\). This shows how the methane and carbon monoxide features are expected to play a much larger role in our expectations of what the planetary signal should be when compared with carbon dioxide.

The only real anomaly in this story is the water features, which have a similar fractional bandpass domination as the carbon monoxide, but are not detected at a similar level in the metric of \(\Delta\sigma\). We attribute this fact to the idea that the water features are dominant at the edges of the spectral bandpass, a region of the spectra which is difficult to model accurately with the spline-continuum, due to increased curvature of the spectral trace in the detector. This implies that the systematic errors in those regions of the spectra are larger than in the central wavelengths where the spline-continuum subtraction is more accurate. If we eliminate the edges from our consideration, we then find that \(f_\textrm{H_2O} = 4.3\%\), which puts it at a lower expected relevance than the CO\(_2\), which is consistent with our findings in \(\Delta\sigma\).

The review article [48] discusses detection speed in the context of high resolution observations of exoplanet atmospheres for certain idealized cases. In their “idealized molecular" case the detection speed is linearly proportional to \(\lambda_\textrm{span}\), the bandwidth of the molecule under consideration. This idealized case consists of many narrow and well separated lines and they acknowledge this inexactitude by noting that certain”wavelength regions will be richer in spectral features than others, depending on the type of planet." While this notion corroborates what we have found with a roughly linear relationship between the fractional bandpass domination for individual molecules and their corresponding impact on the CCF detection sigma we computed with \(\Delta\sigma\), it is important to state this caveat that it would not necessarily be true for opacity sources with qualitatively different features. For example clouds have broadband absorption and lack many narrow features which would help identify them. Additionally the CO\(_2\) and H\(_2\)O features appear to be qualitatively flatter when compared with CO, possibly due to an increase in the density of possible molecular state transitions in the wavelengths due to their triatomic nature, which could contribute a similar effect. This can be seen examining the molecular contribution functions in Figure (3).

4.3 Comprehensive Spectral Analysis↩︎

While the NIRSpec data is sufficient to provide a detection of 51 Eridani b on its own and investigate the molecules present in its atmosphere, it is also useful to combine this dataset with additional ground based spectroscopy and photometry to characterize the planet more comprehensively. Using all of the data described in Section 2, we compute model likelihoods over a grid of model parameters: (\(T_\mathrm{eff}\), \(\log g\), \([\mathrm{M}/\mathrm{H}]\), \(\textrm{C}/\textrm{O}\), \(\log K_\mathrm{zz}\), \(R_\mathrm{P}\), \(\hat{e}\), and \(f_\mathrm{hole}\)). The first five are intrinsic to the model grid itself, the effective temperature, surface gravity, metallicity, carbon to oxygen ratio, and vertical eddy diffusion coefficient, while \(R_\mathrm{P}\) is the planet radius, \(f_\mathrm{hole}\) is the cloud hole fraction, and \(\hat{e}\) is a nuisance parameter which represent the errorbar inflation — a rescaling constant which increases the size of the uncertainties of NIRSpec extracted fluxes. The cloud hole fraction represents the relative weighting between two distinct radiative transfer calculations, one with and the other without clouds [49], \[F_\lambda = (1 - f_\mathrm{hole})F_{\lambda\mathrm{,cloudy}} + (f_\mathrm{hole})F_{\lambda\mathrm{,cloudfree}}\] where the clouds are computed using Virga [50] using a sedimentation efficiency of 10. The errorbar inflation modifies the calculation of the log likelihood by introducing an additional constant \[\log \mathcal{L} \propto N_\mathrm{NIRSpec} \frac{1}{\hat{e}} - \frac{1}{2}\chi^2 \label{eq:likelihood}\tag{1}\] where \(N_\mathrm{NIRSpec}\) is the number of points in the NIRSpec data for which the errorbars are inflated. This term originates from the constant term which normalizes the traditional gaussian likelihood being modified \(\frac{1}{\sqrt{(2\pi\sigma^2)}} \rightarrow \frac{1}{\sqrt{2\pi(\hat{e}\sigma)^2}}\). Here \(\chi^2\) is defined using the formalism laid out previously in [21], where the covariances computed in [20] are used to compute \(\chi^2\) for the GPI IFS spectra in the J, H, K1, and K2 bands, no covariance is used for the photometry from Keck/NIRC2 or JWST/NIRCam, and a new term is added for the NIRSpec data in particular.

\[\chi^2 = \sum_\textrm{GPI} v^T C^{-1} v + \sum_\textrm{phot} \Big(\frac{v}{\sigma}\Big)^2 + \chi^2_\textrm{NIRSpec} \label{eq:chi2}\tag{2}\] Here \(v\) represents a vector of residuals between the model and data \[v = \Big(\frac{R_\textrm{planet}}{d_*}\Big)^2 F_{\lambda,\textrm{model}} - F_{\lambda,\textrm{observation}}.\] The NIRSpec term is slightly modified, using a covariance free estimate in the first iteration, and then a simple covariance model in the second iteration, like so: \[\chi^2_\textrm{NIRSpec} = \sum_{\textrm{NIRSpec}}\begin{cases} \Big(\frac{\hat{v}}{\hat{e}\sigma}\Big)^2, \qquad \textrm{iteration 1} \\ \hat{v}^T C^{-1} \hat{v}, \quad \textrm{iteration 2} \end{cases}\] This is because the slit-based observation strategy lacks adequate spatial sampling to reliably estimate the covariance model from noise spectra of empty space observed with equivalent stellar flux, how it was done previously in IFU based observations. Instead, we estimate the covariance from the autocorrelation of the data-model residuals after the first round of computing the likelihood, so the covariance is included only in the second iteration. The residual vector \(\hat{v}\) is also modified to include spline-subtraction of the model, so it can be appropriately compared to the spline subtracted data. \[v = \Big(\frac{R_\textrm{planet}}{d_*}\Big)^2 \mathcal{SF}(F_{\lambda,\textrm{model}}) - F_{\lambda,\textrm{observation}},\] where \(\mathcal{SF}\) is our representation of the spline filter breads.filter_spec_with_spline function which removes the continuum using the spline model defined previously.

In the first iteration of our spectral model analysis, the likelihood defined in Equation (1 ) is calculated over the entire T-type Elf Owl model parameter space for effective temperature, surface gravity, metallicity, carbon to oxygen ratio, and vertical eddy diffusion coefficient. The exact values of these parameters appear in Table (1), in addition to external parameters for the planet radius, cloud hole fraction, and errorbar inflation which we define.

The best fit model is chosen from this selection as the model with the largest value of the posterior probability, which is the product of the likelihood from Equation (1 ) and the prior probability, which is estimated from the Sonora Diamondback evolutionary model prior [51], [52] and appears in Figure (5). We adopt the mass posterior of 2-4 \(M_\mathrm{Jup}\) from [19] and the age posterior of [40] of \(23.2\pm2.0\) Myr.

The best fit model is shown in Figure (6). A couple of things to note about this spectra: first, that it fits the ground based GPI IFS spectra as well as the mid-IR photometry rather well. Second, while it is difficult to visually verify the correlations between the NIRSpec data and spline-filtered model, which is the purpose of the cross-correlation function, the best fit errorbar inflation parameter \(\hat{e}\) is only \(1.75\), which means the photon noise errorbars estimated from the JWST pipeline are less than factor two off the total error. This total error includes not only photon noise but also systematic errors from imperfections in the spline filter subtraction. Thirdly, the best fit radius \(R_\textrm{P} = 1.37\) \(R_\textrm{Jup}\) and cloud hole fraction \(f_\textrm{hole} = 0.5\) are both larger than we found previously in [21] which reported values of \(R_\textrm{P} = 1.0\) \(R_\textrm{Jup}\) and \(f_\textrm{hole} = 0.2\). This can be attributed to the inclusion of the evolutionary prior in particular for the radius, which now enforces a radius \(R_\textrm{p} \geq 1.35\) R\(_\textrm{Jup}\). This also has the consequence of forcing the effective temperature to lower values to maintain the total luminosity. The larger cloud hole fraction could plausibly result from the interaction between opacity sources, as the newest best fit has a higher metallicity of around 0.5 dex, while metallicity was not at all considered in our previous analysis.

This best fit model is used to estimate the covariance on the errorbars for the NIRSpec data in the following manner. For reference, in previous examples applying this technique on IFU data [8], we estimate the error covariance by computing the correlation profile from “noise spectra," that is, spectra extracted from regions of the observation dataset with similar separations but different position angles than the extracted spectra from the region of the planet. However, with NIRSpec in the slit observing mode, the slit itself functionally obstructs these regions which would have similar noise properties, so we estimate the corelation profile from the model-subtracted data as a proxy for the noise spectra. \[\Psi_{\Delta l} = \langle r_l, r_{l + \Delta l} \rangle_l\] This is just a restatement of Equation (E1) in [8] with the average over spatial positions removed (since we only have one realization of the noise spectra instead of many). This results in a correlation profile which is decidedly more noisy. The correlation profile is shown in Figure (7).

An important observation from this figure is that the 1-2 pixel wide correlation which results from the regular wavelength grid interpolation step of the BREADS chain is still measurable, even with only a single instance of noise spectra. However, the typically seen 20-200 pixel scale oscillations which are interpreted as correlated errors from the spline subtraction are not visible. These errors are almost certainly present, but without more available noise spectra to estimate their magnitude, remain unmodeled and contribute to the errorbar inflation \(\hat{e}\). However, the density of our spline nodes are relatively high, which contributes to mitigating their maximum extent. In order to prevent this noisy correlation profile from affecting the covariance matrix, we enforce \(\Psi(\Delta l \geq 3 \textrm{px}) = 0\) to ensure the stability of the matrix inversion. This curve can be seen with the dashed black line in the figure.

This covariance should be an improvement over ignoring the correlated errors completely, as it accounts for the part of the known error correlation due to the regular wavelength grid interpolation. This covariance model is used in the second iteration of model fitting, which also employs sub-model-grid interpolation to better sample the likelihood contours in the region around the best fit model.

We perform the sub-grid interpolation on the discrete parameter vectors defined in the last column of Table (1) using the module breads.atm_utils.broadRGI which internally relies on scipy.RegularGridInterpolator to perform the multi-linear interpolation over the input parameter space. While this methodology is rather commonplace it is also known to be problematic and could potentially be improved by replacing it with a more sophisticated interpolation scheme such as STARFISH [53]. The posterior distributions calculated during the second iteration of the likelihood calculation are shown in Figure (8).

One observation to note is that the best fit model from the second iteration has slightly different parameter values when compared with the first iteration. In this case, we find that \(T_\mathrm{eff}\) = 800 K, \(\log g\) = 3.75, \([\mathrm{M}/\mathrm{H}]\) = 0.7, \(\textrm{C}/\textrm{O}\) = 0.458, \(\log K_\mathrm{zz}\) = 3, \(R_\mathrm{P}\) = 1.36, \(\hat{e}\) = 1.74, and \(f_\mathrm{hole}\) = 0.3. Overall, the temperature, gravity, metallicity, C/O ratio, and eddy diffusion coefficient are all larger while the hole fraction is smaller. While the general structure of the atmosphere models are nearly indistinguishable by eye, this model predicts fainter values for the mid IR photometry, which may be more accurate given the limits set by [22]. From the first iteration, our synthetic photometry of the model predicts a flux value that would have been detected at 2.87\(\sigma\) in F430M, while the second iteration predicts instead a detection at 2.24\(\sigma\). While both of these are below the 3\(\sigma\) threshold that would be prudent to claim a detection, the first is closer to that threshold whereas the second is more marginal. Neither model predicts a significant detection in F460M, with an expectation at 1.05\(\sigma\) and 0.79\(\sigma\) for the first and second iterations respectively.

Regardless, both iterations have substantial preference for models with heightened metallicity relative to solar abundances. The first iteration prefers [M/H] = 0.5, while the second instead prefers [M/H] = 0.7, but with a 1\(\sigma\) range that includes 0.5. This slightly higher metallicity is responsible for the darker color in F430M, as carbon dioxide abundance is very sensitive to metallicity, which is the dominant absorber in this band. [22] did not report flux measurements for these bands, instead opting to report the more conservative upper limit. It is possible that absolute flux measurements at the planet’s position in JWST/NIRCam F430M photometry could provide some additional discriminatory power over the reported upper limits, even if those flux measurements have large errorbars.

While the full triangle shaped posterior is useful to evaluate the possible covariances between model parameters, it is also interesting to ask how the posterior distributions for those model parameters depending on which datasets you choose to include in the calculation of the Likelihood. Equation (2 ) which computes the \(\chi^2\) has three primary terms \(\chi^2_\textrm{GPI}\), \(\chi^2_\textrm{phot}\), and \(\chi^2_\textrm{NIRSpec}\) which we sum to compute the total \(\chi^2\) which then is proportional to the log likelihood. We could have instead chosen different subsets of these datasets to include in this computation, which effectively results in performing the same inference problem on different combinations of the data. Since there are three primary categories of data, there are seven possible combinations of these datasets we consider: \[\begin{align} \log\mathcal{L}_\textrm{all} \propto - 1/2(\chi^2_\textrm{GPI} + \chi^2_\textrm{phot}+\chi^2_\textrm{NIRSpec}) + N_\mathrm{NIRSpec}/\hat{e} \\ \log\mathcal{L}_\textrm{GPI} \propto - 1/2(\chi^2_\textrm{GPI} \color{white}+ \chi^2_\textrm{phot}+\chi^2_\textrm{NIRSpec}\color{black}) \color{white} + N_\mathrm{NIRSpec}/\hat{e}\color{black} \\ \log\mathcal{L}_\textrm{phot} \propto - 1/2(\color{white}\chi^2_\textrm{GPI} +\color{black}\chi^2_\textrm{phot}\color{white}+\chi^2_\textrm{NIRSpec}\color{black})\color{white} + N_\mathrm{NIRSpec}/\hat{e}\color{black} \\ \log\mathcal{L}_\textrm{NIRSpec} \propto - 1/2(\color{white}\chi^2_\textrm{GPI} + \chi^2_\textrm{phot}+\color{black}\chi^2_\textrm{NIRSpec}) + N_\mathrm{NIRSpec}/\hat{e} \\ \log\mathcal{L}_\textrm{GPI+phot} \propto - 1/2(\chi^2_\textrm{GPI} + \chi^2_\textrm{phot}\color{white}+\chi^2_\textrm{NIRSpec}\color{black})\color{white}+ N_\mathrm{NIRSpec}/\hat{e}\color{black} \\ \log\mathcal{L}_\textrm{GPI+NIRSpec} \propto - 1/2(\chi^2_\textrm{GPI} \color{white}+ \chi^2_\textrm{phot} \color{black}+\chi^2_\textrm{NIRSpec}) +N_\mathrm{NIRSpec}/\hat{e} \\ \log\mathcal{L}_\textrm{phot+NIRSpec} \propto - 1/2(\color{white}\chi^2_\textrm{GPI} + \color{black}\chi^2_\textrm{phot}+\chi^2_\textrm{NIRSpec}) +N_\mathrm{NIRSpec}/\hat{e} \\ \end{align}\]

It is important to remark that the additional term \(N_\textrm{NIRSpec}/\hat{e}\) only appears in combinations where the NIRSpec data is relevant, otherwise the \(\hat{e}\) parameter is meaningless since it only applies to the NIRSpec data. For each of these combinations of datasets, we compute the single parameter marginal posterior distributions and show them for comparison in Figure (9).

One of the observations of note within this figure is the tension between the photometry and the NIRSpec posteriors, in particular for temperature, metallicity and cloud hole fraction. On its own, the NIRSpec data seems to prefer a higher effective temperature, lower metallicity, and larger cloud hole fraction, all of which together imply that the NIRSpec data on its own is better fit by a model which is brighter than would be expected when comparing with the photometry. Ultimately, in our likelihood framework, we assumed some unmodeled errorbar inflation \(\hat{e}\) on the NIRSpec data, to deal with our inability to model the error covariance properly. However, this may have the consequence of also giving a degree of freedom to ignore what the NIRSpec data is telling the model would be best. One can speculate that if instead we had parameterized our likelihood with a floating offset parameter on the photometry as a plausible way to instead correct for inter-instrument systematics, we might see a substantially different answer for the best fit model.

The cloud hole fraction degree of freedom gives the model substantial power to adjust the relative intensity of the planet brightness in the near IR (1-3 microns) relative to the mid IR (3-5 microns) because the cloud opacity dominates at the shorter wavelengths, which is easily seen in Figure (2) of [21]. Since the prior we assumed for this analysis (Figure (5) has a relatively narrow distribution for the planet radius, the NIRSpec only posterior is constrained into other parameters which could inflate the total brightness of the planet at these wavelengths, notably temperature, metallicity, and hole fraction.

5 Discussion↩︎

We have demonstrated the capabilities of an uncommon observing mode of JWST / NIRSpec. By using moderate resolution (R \(\sim\) 2700) spectroscopy combined with software based starlight suppression techniques, we are able to detect a signal from 51 Eridani b at 4.8\(\sigma\) in 3.43 hours of exposure time (\(\sim\)​10 hours including overheads) without using a coronagraph. The planet is separated by \(\sim\)​0.3 arcsec from its K = 4.5 magnitude host star with a global (planet-star) contrast ratio of C=1.1\(\times10^{-5}\) and a local (planet-speckle) contrast ratio of Q=1.4\(\times 10^{-3}\), the most challenging target we have detected with this technique. The signal is localized in the slit spatial coordinate (IFU Y) - planetary radial velocity phase space at the expected position of the planet. The peak of the of the cross-correlation function is spatially offset from the peak of the stellar flux, indicating that it is not an artifact of amplified stellar noise. The signal is absent (1.8\(\sigma\)) from the speckle-side pointing dataset. The signal is robust to analysis hyperparameters such as the node spacing, persisting at similar levels of confidence when reducing the data over a wide range of this parameter. “Leave one molecule out" detection tests indicate this signal is primarily driven by the presence of methane (\(\Delta\sigma_\textrm{CH_4} = -3.43\)) and carbon monoxide (\(\Delta\sigma_\textrm{CO} = -1.81\)), which are expected from theoretical models to be the most prominent absorbers at these wavelengths and do not exist in the atmosphere of the stellar host. All of these lines of evidence and testing provide reassurance that the signal is indeed of planetary origin. Additionally, this detection is the first direct evidence of both molecules (CO + CH\(_4\)) coexisting in the planet’s atmosphere, a signature of disequilibrium chemistry due to turbulent vertical mixing. The presence of disequilibrium chemistry has been suspected since the original discovery of methane [15] in its atmosphere, and with our current detection predicated on the existence of CO in the planet’s atmosphere, equilibrium chemistry is no longer a viable explanation. This is in contrast to the previous analysis in [21], which did not have enough statistical evidence to rule out equilibrium chemistry models using only the single photometric point at M band.

However, our analysis is not perfect and our methods could be further refined and improved. Our comprehensive atmospheric model analysis indicates that the best fit NIRSpec errorbar inflation parameter \(\hat{e} \sim 1.75\), which indicates that the photon noise errorbars on the data taken from the JWST pipeline underestimate the total error by roughly a factor of two. The primary source of the extra noise is almost certainly systematic error from the spline-based continuum subtraction which is known to result in spectrally correlated errors from previous work, although additional contributions to the error could come from imperfect calibrations of wavelength and spectral curvature in the NIRSpec instrument model used by the JWST pipeline. Future work on improving the spline modeling and continuum subtraction in general could reduce the systematic error term. In the context of doing similar observations with the NIRSpec IFU instead of the fixed slit, the covariance can be estimated from the data directly. By sampling the “noise spectra" at spatial positions with similar flux levels as the position of the planet (the”noise annulus") and computing their correlation profile, we can estimate the level of this covariance and include that model in our calculation of the likelihood. However, the primary disadvantage of using the fixed slit is that it does not sample those particular representative spatial positions, making this component of the covariance difficult to estimate. It is possible that with further work on this front, it could be possible to build a better covariance model from initial principles, with blind testing on different datasets with known estimations for validation. Application of such a covariance model could provide improvements in the computation of the cross correlation function and the full model likelihood, potentially resolving some discrepancies we observed or increasing the detection significance. We leave this question open for future inquiry. It is also possibly worth considering alternative modes of spline-subtraction which build into the model knowledge of the three dimensional structure of the JWST PSF, instead of our simple one dimensional treatment of the detector rows.

There is another downside to our observing strategy which warrants mention. Using the fixed slit instead of the IFU requires an additional pointing, the “speckle-side" dataset, to compute the stellar spectrum and provide a null hypothesis test. This inflates the observing time required at a fixed SNR by a factor of two as the IFU observation strategy can effectively perform both pointings simultaneously. Our observations were planned with the fixed slit originally because of concerns with the detector saturating due to the extremely bright nature of the target star, however, with additional experience using these modes, we have begun to develop software strategies to mitigate detector saturation and charge bleeding from those saturated pixels. This program was designed during cycle 1, before the initial demonstration of JWST/NIRSpec in this mode on HD 19467 B had been finalized. While we have shown that using the fixed slit arrangement is a viable observing strategy, we would now recommend using the IFU over the fixed slit for the majority of high contrast spectroscopic targets. The charge bleeding effect can be mitigated as long as the planet is in the IFU slice adjacent to the primary slice which the star is contaminating. This was demonstrated in the Cycle 1 analysis on HD 19467 B [8] using a spline-based fit to the dark regions between the IFU traces. Improvements to this charge bleed modeling and subtraction using a Lorentzian-convolution subtraction are explored in (Ruffio & Xuan 2025 in prep.) and continued to be investigated with further work on improving the analysis strategies.

Another interesting open question remains – what would be the best strategy for optimal follow up observations of 51 Eridani b? One molecule that we searched for but only found a “consistent non-detection" for is carbon dioxide (CO\(_2\)). The trade off between doing a deep F430M search with NIRCam to search for the depth of the CO\(_2\) feature and further refine the constraint on the planet’s metallicity could be compared with a longer spectroscopic integration using NIRSpec. How do both the absolute flux continuum level and total fractional bandpass coverage for a molecule compete in the effective detection SNR landscape as a function of exposure time? While NIRcam data was able to infer the existence of this molecule from a non-detection, our model fitting exercise indicates this measurement is possibly close to the edge of detectability, with expected model fluxes in the \(2\)-\(3\sigma\) range in this band at the current NIRCam sensitivity. However, direct measurement of the CO absorption band with NIRCam was not close to the threshold, \(\sim1\sigma\). At the very least, the NIRCam and NIRSpec data are complementary in that they are sensitive to different qualitative aspects of the planet’s spectra. While our analysis presents some tension between the measurements between these two datasets, it is possible this tension arises from the unmodeled error covariance in our data.

Overall, we are just beginning to explore the capabilities of this novel application of direct spectroscopy to exoplanets. Being able to directly detect the presence of carbon monoxide from disequilibrium chemistry is one highlight, and the possibility to explore this technique at other wavelengths searching for other molecules and their implications is an unexplored frontier. While our software analysis techniques continue to be improved, it is possible that future reanalysis of these same datasets could provide improved results.

6 Acknowledgments↩︎

In no particular order, A.M. would like to thank William Balmer, Mathilde Malin, Sagnick Mukherjee, Natasha Batalha, Caroline Morley, and Laurent Pueyo for their helpful conversations. This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. Support for program JWST-GO-03522 (A.M., J.-B.R.) was provided by NASA through grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute. This work benefited from the 2024 Exoplanet Summer program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation.

7 Data/Code Availability↩︎

The JWST/NIRSpec data analyzed in this paper is publicly available in MAST, as well as the raw JWST/NIRCam data. The auto-reduced NIRSpec datasets can be referenced though the DOI http://dx.doi.org/10.17909/zj2p-hf65, while our code contains a script to obtain the level 1 uncalibrated data products using the helper tool jwst_mast_query.py. The reduced GPI, JWST/NIRSpec, and JWST/NIRCam data products, as well as the scripts used to analyze the data and produce the figures in this paper, are available at the following repository: https://github.com/muuud/51-Eri-b-NIRSpec or in Zenodo at the following DOI http://dx.doi.org/10.5281/zenodo.17101634.

8 Software↩︎

BREADS [23], JWST [24], PICASO [47], VIRGA [50], astropy [54]–[56], species [45], matplotlib [57], numpy [58], scipy [59]

9 Appendix: Injection and Recovery Test↩︎

In order to provide another line of evidence for the detection of 51 Eridani b, we perform an artificial injection and recovery test using only the speckle-side datasets. We inject a spline-filtered Elf-Owl model with parameters \(T_\mathrm{eff} = 750\) K, \(\log \big( \frac{g}{\mathrm{cm}/\mathrm{s}^2} \big) =\) 3.5, \([\mathrm{M}/\mathrm{H}] = 0.5\), \(\textrm{C}/\textrm{O} = 0.458\), and \(\log \big(\frac{K_\mathrm{zz}}{\mathrm{cm}^2\mathrm{s}}\big) = 2\), into the speckle-side continuum-subtracted point clouds using a synthetic WebbPSF interpolation after converting the model into the detector units of MJy. Then we perform the spectral extraction and cross-correlation analysis on the perturbed speckle-side point clouds following the same procedure as the main analysis. The resulting phase space CCF from this injection and recovery test is shown in Figure (10) demonstrating that the fake signal is recovered at a significance of \(5.0\sigma\)

This detection significance is very similar to but slightly higher than the real dataset, which could be plausibly attributed to data-model mismatch, as the injected signal is a perfect match to the CCF template used in this analysis which is not necessarily true of the real data. Additionally, the speckle-side only residual noise does not contain minor errors from spectral mismatch that the planet-side dataset has. This spectral mismatch could arise from the difference in the slit illumination profile between the planet- and speckle-side datasets.

10 Appendix: Additional Figures↩︎

References↩︎