Detecting Exomoons from Radial Velocity Measurements of Self-luminous Planets: Application to Observations of HR 7672 B and Future Prospects

Moons similar to those around Jupiter are expected to form in circumplanetary disks (CPDs) as a by-product of planet formation (Batygin & Morbidelli 2020). The typical CPD total dust mass relative to the planet is around 10⁻⁴ (Canup & Ward 2006) commensurate with the mass ratios of solar system satellites around the gas giants listed in Table 1. This is consistent with the measured value of the CPD around PDS 70 c from Atacama Large Millimeter/submillimeter Array observations (Benisty et al. 2021), which is 0.031 M_Earth and corresponds to a mass ratio of about 5 × 10⁻⁵, assuming a 2 M_Jup planet (Benisty et al. 2021). It is also possible to form larger moons from the merger of Galilean-like multiple systems. This is the proposed scenario to explain the high eccentricity and large mass of Saturn's moon Titan (Asphaug & Emsenhuber 2018). Alternative formation pathways include the capture of satellites (e.g., Neptune's moon Triton; Agnor & Hamilton 2006), collisions with protoplanets (e.g., the Moon; Canup & Asphaug 2001), or even gravitational instability like in the formation of brown-dwarf binaries (Lazzoni et al. 2020). The detection of satellites that formed in a CPD remains challenging with current instrumentation, but binary planets and brown dwarfs are already accessible with various techniques depending on their separation (Lazzoni et al. 2022). Characterizing the different populations with different mass ratios, such as binary brown-dwarf companions (i.e., triple systems) and smaller CPD moons could help inform the formation pathways of directly imaged planets. Indeed, binary companions could only occur in a top-down scenario, while CPD formation could occur in all cases. An important lesson from early exoplanet discoveries is that the solar system is not a good predictor of exoplanet demographics. Planet formation theories also often struggle to account for the diversity of new discoveries. For example, the first discoveries of hot Jupiters and the ubiquity of super-Earths and mini-Neptunes were initially a surprise to the community (Batalha 2014). As a corollary, it would be unwise to assume that exomoon searches should be any different (Kipping et al. 2015). It is therefore important that we keep pushing the discovery space with new observational methods.

For the past decade, transiting surveys have unequivocally dominated the landscape of exomoon searches (Kipping et al. 2012) through the analysis of transit timing variations and additional transit signal from the moon. They have placed the first constraints on exomoon occurrence rates and shown that high-mass-ratio satellites are not common around short-period exoplanets (Kipping et al. 2015; Teachey et al. 2018). Other detection techniques have been used to look for exomoons around directly imaged exoplanets such as astrometry or direct imaging. While these terms also refer to planet-detection techniques, in this context, direct imaging means to spatially resolve the satellite from the planet. Astrometric detections refers to the measurement of the astrometric wobble of a planet caused by orbiting moons with precise interferometric instruments such as the Very Large Telescope Interferometer (VLTI)/GRAVITY (Gravity Collaboration et al. 2021). To date, only a handful of exomoon candidates have been proposed: for example, two around transiting planets (Kipping et al. 2022; Teachey & Kipping 2018), one orbiting a directly imaged brown dwarf (Lazzoni et al. 2020), and another around an isolated planetary-mass object (Limbach et al. 2021). None have been confirmed. Most notably, the exomoon candidate Kepler 1708 b-i is a transiting 2.6 Earth radii object orbiting a Jupiter-sized planet, which, if confirmed, would be several orders of magnitude larger than the Galilean moons in terms of mass ratio (Kipping et al. 2022).

Transiting planets generally have short orbital periods and smaller Hill spheres. These conditions could be less favorable to moon formation and retention, while observed transits of long-period exoplanets (>1 au) are intrinsically rare. For example, it has been suggested that planets could lose their satellites as they migrate inward (Spalding et al. 2016). The detection of exomoons around imaged planets with astrometry or direct imaging is more sensitive to longer-period moons (Lazzoni et al. 2022). These two techniques might not be well suited for satellites with mass ratios and separations (<30 R_Jup) similar to the ones orbiting the solar system gas giants.

Another technique has been proposed to look for moons around directly imaged planets using radial velocity (RV) measurements of the planet itself (Vanderburg & Rodriguez 2021; Vanderburg et al. 2018). By measuring the wobble of planets caused by orbiting satellites, planetary RV surveys are a promising alternative for finding Galilean moon analogs. Directly imaged companions are likely to have a different formation and migration history compared to transiting exoplanets. They are also generally more massive and farther away, resulting in much more extended Hill spheres. Models suggest that larger planets form even larger moons following the scaling m ∝ M^3/2, with m and M the masses of the moon and the planet, respectively (based on Equation (43) in Batygin & Morbidelli 2020). As another formation pathway, binary systems forming as the tail end of stellar formation through gravitational instability of the protostellar cloud or the protoplanetary disk could lead to a population of easily detectable high-mass ratio satellites and binary companions (Lazzoni et al. 2022). In summary, directly imaged companions could be more likely to host larger moons, which could be detected from RV monitoring of the planets themselves.

There are two aspects to optimizing the choice of a target for exomoon searches: the RV precision that can be achieved and the probability of the object to host a moon. In this work, we focus on the former and study the detectable mass ratios for exomoons as a function of the instrument, the telescope, and the planet or brown-dwarf properties. Although the RV detection of exomoons remains challenging due to the intrinsic faintness of planets and the light contamination from the glare of the host star, the expertize obtained from 30 years of stellar RV exoplanet detections is an invaluable asset. Indeed, stable high-resolution spectrographs and data-analysis techniques are already demonstrating stability and performance in excess of the level of precision needed for exomoons. Vanderburg & Rodriguez (2021) derived the first exomoon mass upper limits with this technique around the HR 8799 planets based on the planetary RV time series reported by Ruffio et al. (2021) using OSIRIS, an R = 4000 spectrograph, at the W. M. Keck Observatory. This study ruled out moons with mass >1 M_Jup and period less than 1 day that would be orbiting the 7 M_Jup planet HR 8799 c.

Recent technological advances in infrared high-resolution spectroscopy for high-contrast companions are enabling the first planetary RV searches for exomoons (Delorme et al. 2021; Jovanovic et al. 2017; Otten et al. 2021; Snellen et al. 2015). The Keck Planet Imager and Characterizer (KPIC) recently demonstrated R = 35,000 K-band spectroscopy of directly imaged exoplanets, including the HR 8799 system, measuring their RVs and obtaining spin measurements for the first time (Wang et al. 2021b). KPIC is the first implementation of a new class of spectrographs that combines the power of the Keck II adaptive optics (AO) systems, the stability and starlight suppression of single-mode fibers, and the high spectral resolution of the NIRSpec spectrograph for the detailed study of directly imaged planets (Delorme et al. 2021). We observed a bright brown-dwarf companion (HR 7672 B, K = 13.04, ∼07; Boccaletti et al. 2003; Liu et al. 2002) as part of the commissioning and science verification of KPIC (Delorme et al. 2021; Wang et al. 2022a). While it is at the boundary of the stellar regime, HR 7672 B is an interesting benchmark companion, because it has a dynamically measured mass of 72.7 ± 0.8 M_Jup (Brandt et al. 2019; Crepp et al. 2012) and its composition should be similar to that of the star due to its assumed formation history from gravitational instability. Using this target, Wang et al. (2022a) showed that accurate atmospheric compositions could be retrieved using KPIC's high resolving power and angular resolution by demonstrating a 1.5σしぐま consistency between the composition of HR 7672 B and its host star (see also Xuan et al. 2022). HR 7672 B was also observed for a full night with KPIC as a test case for variability studies. This time series can be used to put the deepest limits to date on the mass of an orbiting satellite around the substellar companion, which we are demonstrating in this work. KPIC is already undergoing several upgrades, including a laser frequency comb that will enable precise RV science (Jovanovic et al. 2020; Yi et al. 2016). The expected doubling of the instrumental throughput will significantly improve its sensitivity (Echeverri et al. 2022; Jovanovic et al. 2020). The next generation of high-contrast, high-resolution spectrographs such as the High-resolution Infrared Spectrograph for Exoplanet Characterization (HISPEC) at the W. M. Keck Observatory and the Multi-Object Diffraction-limited High-resolution Infrared Spectrograph (MODHIS) on the future Thirty Meter Telescope (TMT) will undoubtedly open new frontiers in this field by allowing 0.98–2.46 μみゅーm simultaneous coverage at an average spectral resolution R > 100,000 (Mawet et al. 2019, 2022).

In Section 2, we present exomoon RV detection limits for the brown-dwarf companion HR 7672 B using KPIC. In Section 3, we then simulate observations of the same brown dwarf and the planet HR 8799 c with next-generation facilities and compare their sensitivity to the moons in the solar system. In Section 4, we explore the parameter space of satellites that could be detected with TMT/MODHIS as a function of planet properties. Finally, we conclude on the prospects for RV detections of exomoons in Section 5.

The brown-dwarf companion HR 7672 B was observed three times in 2020 and then for a full night on 2021 July 4 with KPIC (R ∼ 35,000) in the K band (1.9–2.4 μみゅーm); Delorme et al. 2021; Mawet et al. 2017). These observations are detailed in Table 2. The first three epochs included 1 to 2 hr of on-target exposures per night and were already published in Wang et al. (2022a) and Delorme et al. (2021). Unfortunately, the conditions on 2021 July 4 were well below average, with the companion undetectable in some individual 5 minute exposures. During this one night specifically, we used an ABAB pattern to nod the companion between two KPIC fibers, Fiber 1 and Fiber 2, to limit or identify any fiber-specific biases. There was no nodding during the other epochs. The data were reduced with the KPIC data-reduction pipeline (DRP) ¹⁹ following the same approach described in Wang et al. (2021b, 2022a). The first steps include background subtraction, bad pixel correction, and the calibration of the fiber trace location and width on the detector for each NIRSpec spectroscopic order. Optimal extraction is then used to extract the spectra and the wavelength solution is derived from the telluric and stellar lines of an M giant, namely HIP 81497, taken on the same night. For this purpose, the telluric model is generated with the Planetary Spectrum Generator (Villanueva et al. 2018) and the wavelength calibrator star is modeled by a Phoenix model: $\mathrm{log}(g/[1\,\mathrm{cm}\,{{\rm{s}}}^{-2}])=1;$ T_eff = 3600 K (Husser et al.2013).

Table 2. K‐band Observations of HR 7672 B with KPIC

Object	Date	Exposure Time	Seeing	Throughput
HR 7672 B	2020-06-08	11 × 10 minutes	04	1%
HR 7672 B	2020-06-09	10 × 10 minutes	06	1.5%
HR 7672 B	2020-09-28	7 × 10 minutes	04	2.7%
HR 7672 B	2021-07-04	61 × 5 minutes	1''	2%

Note. The quoted throughput is the end-to-end from the top of the atmosphere, which is a better proxy of performance than the seeing for KPIC.

Download table as:  ASCIITypeset image

We use a forward-modeling approach similar to Wang et al. (2021b) and Ruffio et al. (2021) to measure the RV of HR 7672 B, which includes a joint modeling of the starlight and the companion signal. Wang et al. (2021c) showed that the continuum could be included in the forward model with a fourth-order polynomial, therefore not requiring the data to be high-pass filtered nor continuum normalized. In this work, we model the continuum using a spline-based linear model, which can be analytically marginalized using the general purpose Python module breads ²⁰ (Broad Repository for Exoplanet Analysis, Discovery, and Spectroscopy) based on the formalism in Ruffio et al. (2019). The spline forward modeling has the advantage of being more robust to bad pixels than a Fourier-based high-pass filter and avoids the nonlinearity of a sliding-window median filter. The spline parameters are also easier to optimize than the coefficients of a high-order polynomial, for example.

We define the forward model as

$$\begin{eqnarray}&&{\boldsymbol{d}}={{\boldsymbol{M}}}_{\mathrm{RV}}{\boldsymbol{\phi }}+{\boldsymbol{n}},\end{eqnarray} \tag{ 1 }$$

where d is the data vector of size N_d , M _RV is the linear model, ϕ are the linear parameters, and n is a random vector of the noise with a diagonal covariance matrix Σしぐま. Off-diagonal elements in the covariance matrix are neglected here, but subsequent data-processing steps would correct for this inaccuracy. The different column vectors of the linear model are illustrated in Figure 1. The data vector and the standard deviation of the noise used to define Σしぐま₀ are direct outputs of the KPIC data-reduction pipeline. The variance of the noise is multiplied by a free parameter scaling factor s² that is fitted to account for any underestimation of the noise.

Zoom In Zoom Out Reset image size

KPIC includes four single-mode fibers separated by 08 on a line. We can therefore acquire simultaneous spectra of the companion and the host star, more specifically the speckle field, by rotating the field of view using the Keck II AO system front-end K-mirror rotator. The observations of the speckle field are used to derive simultaneous empirical models of the transmission and the starlight spectra used in the forward model. The starlight is used to model the speckle noise leaking into the fiber at the position of the companion. The wavelength calibration is different in each fiber so the spectra are linearly interpolated to match the sampling of the science fiber. The planet model is defined as the spin-broadened best-fit model from Wang et al. (2022a) using petitRADTRANS (Mollière et al. 2019) multiplied by the empirical telluric and instrument transmission profile. The continuum of both the planet and the speckle are modulated by a third-order spline model. Ten spline nodes are used in each spectral order (Δでるたλらむだ ∼ 0.05 μみゅーm) for the planet model to manage any inaccuracies in the continuum due to imperfections in the atmosphere model fit. This number of nodes is analogous to a 200 pixel wide high-pass filter. The number of nodes was chosen as a trade-off between the number of additional parameters and the optimal high-pass filter scale of 100 pixels found in Xuan et al. (2022). The speckle continuum is modeled with three spline nodes to model any speckle crossing the fiber location as the wavelength changes. This results in 13 linear parameters per spectral order representing the values of the continua at the location of the nodes (see Figure 1). This defines the linear model M _RV with dimensions N_d × 13, which is also a function of the RV of the planet, the only nonlinear parameter fitted for here.

KPIC data feature strong spectral fringing due to the Fabry–Pérot cavities formed by the transmissive optics inside the NIRSpec spectrograph (Hsu et al. 2021) and within the KPIC fiber-injection unit (Finnerty et al. 2022). This effect is made worse by the high spatial coherence of the wave front in KPIC. We therefore apply a Fourier filter to the data and the forward model by zeroing frequencies corresponding to the fringes. A physical model of the fringing, such as in Cale et al. (2019), could be explored in the future.

The likelihood function is defined from a multivariate Gaussian distribution as

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }(\mathrm{RV},{\boldsymbol{\phi }},{s}^{2})=\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{{N}_{d}}| {{\boldsymbol{\Sigma }}}_{0}| {s}^{2{N}_{d}}}}\\ \quad \exp \left[-\displaystyle \frac{1}{2{s}^{2}}{\left({\boldsymbol{d}}-{{\boldsymbol{M}}}_{\mathrm{RV}}{\boldsymbol{\phi }}\right)}^{\top }{{\boldsymbol{\Sigma }}}_{0}^{-1}({\boldsymbol{d}}-{{\boldsymbol{M}}}_{\mathrm{RV}}{\boldsymbol{\phi }})\right].\end{array}\end{eqnarray} \tag{ 2 }$$

The likelihood is maximized using a linear least-square solver on a grid of RV values from −400 to 400 km s⁻¹ in steps of 0.2 km s⁻¹. The 1σしぐま RV uncertainties are derived from the RV posterior calculated analytically according to Equation (10) in Ruffio et al. (2021) on this RV sampling. This method analytically marginalized the RV posterior for the modulation of the continuum and the noise scaling factor. The linear spline parameters used to fit the continuum are forced to be positive. This is theoretically inconsistent with the framework, which assumes unconstrained parameters, but it does not appear to significantly impact the RV time series.

Only the three reddest orders, out of nine in the K band, are used in this analysis. The bluest three orders (numbered 39–37; 1.94–2.09 μみゅーm) were discarded because they feature strong saturated CO₂ telluric lines that are generally harder to model, but also make for an unstable fit due to overlapping frequencies with the fringing and the simple Fourier filter. The middle three orders (2.10–2.27 μみゅーm) lack sufficient stellar and telluric spectral lines to calibrate the wavelength precisely enough. Thus, only the remaining three orders are used in this analysis: 2.29–2.34 μみゅーm (order 33), 2.36–2.41 μみゅーm (order 32), and 2.44–2.49 μみゅーm (order 31). Order 33 includes the CO bandhead and therefore results in the strongest signal-to-noise ratio (S/N) and the most precise RV measurement. Each NIRSpec spectral order is fitted separately resulting in three RV estimates for each exposure.

The barycentric corrected RV measurements for the four epochs and three orders are shown in Figure 2. Following the method described in Section 2.2, the median RV uncertainties in 5 minute exposures are 2.5 km s⁻¹, 4.0 km s⁻¹, and 4.9 km s⁻¹ for orders 6, 7, and 8, respectively. We overplot the predicted RV of the brown dwarf from orbital fits to the relative astrometry from Crepp et al. (2012) and RV measurements of the host star (Crepp et al. 2012; Rosenthal et al. 2021). The orbit fits were done with orbitize! (Blunt et al. 2020) following its RV tutorial ²¹ and using the emcee (Foreman-Mackey et al. 2013) sampler to obtain a posterior of allowed orbits. This orbital RV of the companion in each epoch is predicted from this orbit fit and is subsequently subtracted from the estimated RV of the planet when running the exomoon search. Similarly to fitting the centroid of a Gaussian (King 1983), the RV precision goes as the typical line width in the spectrum divided by the total S/N of the detection. In the case of HR 7672 B, the large spin with $v\sin i=45.0\pm 0.5\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Wang et al. 2022a) is a limiting factor in deriving more precise RVs. The impact on the exomoon sensitivity of other fundamental parameters such as the brightness, age, mass, and separation from the star are discussed in Section 4 in the context of TMT/MODHIS.

Zoom In Zoom Out Reset image size

The open-source Python package RVSearch ²² (Rosenthal et al. 2021) is used to look for possible exomoons around HR 7672 B and derive the sensitivity of our KPIC RV time series. RVSearch is a planet-search algorithm that was developed by the California Legacy Survey for high-precision RV surveys (Fulton et al. 2021; Howard & Fulton 2016; Rosenthal et al. 2021). Planets are detected from periodograms, which are expressed as the difference in Bayesian information criterion (BIC) between a model including the planet and a model without it (Rosenthal et al. 2021). The ΔでるたBIC can be used to select the model that best represents the data, or, in other words, determine if a planet is necessary to explain the observations. Planet candidates are detected by iteratively adding additional planet signal to the model (Rosenthal et al. 2021). For each iterative search, the algorithm fits a detection threshold to the periodogram using the power-law noise model described in Howard & Fulton (2016). To characterize the search completeness of a data set, RVSearch performs injection-recovery tests, drawing many synthetic planet signals, injecting them in the data, and checking whether their signals surpass the last detection threshold. The simulated signals were injected as described in Rosenthal et al. (2021) with period and $M\sin i$ from log-uniform distributions, and eccentricity from an empirically calibrated beta distribution (Kipping 2013).

RVSearch is directly applicable to the search for exomoons by replacing the properties of the star by the ones of the planet. We assume that each spectral order in NIRSpec has a different zero RV point due to possible inconsistencies between them. This can be done with RVSearch, which linearly solves for offsets between subsets of RVs, and uses a wide, Gaussian, uninformative prior on white noise for each subset. This feature is usually used to fit data from different instruments. Two analyses are performed, first only using the long night of observations (2021 April 7) and then all the available data. The latter provides a longer time baseline. The resulting periodograms and exomoon completeness are shown in Figure 3. By combining the four epochs, the observations are sensitive to satellites with a mass ratio of 1% at semimajor axes similar to that of Io (6 R_Jup) around Jupiter or 4% at the distance of Callisto (15 R_Jup). While these are encouraging results, the smallest detectable satellites would be as large as Jupiter due to the already large mass of HR 7672 B. As shown in Section 4, targeting smaller brown dwarfs and planets does not generally allow the detection of moons with smaller absolute masses, because the S/N drops faster than the mass of the object due to the decreasing brightness. If satellites around HR 7672 B were to orbit within ∼10 R_Jup of the brown dwarf, they would likely fall within the Roche radius (see Figure 4). Such satellites would be tidally disrupted and likely result in the formation of rings around the planet. It is possible that this issue would prevent the formation of a resonant chain of satellites if the inner edge of the decretion disk falls within the Roche limit. This is, for example, cited as a possibility to explain the difference between the Galilean and the Saturnian satellite systems in Batygin & Morbidelli (2020). At the other end of possible satellite semimajor axes, stable orbits can generally exist up to one half of the Hill sphere for prograde orbits (Shen & Tremaine 2008). The Hill sphere of HR 7672 B being r_H ≈ 5.6 au = 1.2 × 10⁴ R_Jup, time series like these ones will not be sensitive to the vast majority of possible orbits without observations spanning years or decades.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we simulate observations from current and future instrumentation at the Keck Observatory and the TMT to estimate the properties of putative satellites that should be detectable using planetary RVs. We use an instrument and observation simulator called PSIsim, ²³ which was first developed for the Planet Systems Imager (PSI; Fitzgerald et al. 2022) instrument concept for TMT and then expanded to include other instruments and telescopes. PSIsim is first used to estimate the RV precision. Then, RV times series are simulated assuming six full nights of observations over 25 days, and the exomoon sensitivity is finally computed using RVSearch. These simulations are meant to represent an ideal scenario in terms of instrument performance and telescope time allocation.

We simulate observations of two substellar companions, the brown-dwarf companion HR 7672 B and the planet HR 8799 c, with four generations of instruments. An exhaustive analysis of all directly imaged companions is beyond the scope of this work so HR 8799 c was chosen as a representative example of the field with a planetary mass. HR 8799 is also the only other high-contrast system with published RV time series and exomoon upper limits (Vanderburg & Rodriguez 2021). The four instruments considered in this work are Keck/KPIC I, Keck/KPIC II, Keck/HISPEC, and TMT/MODHIS. KPIC I corresponds to observations carried out pre-2022A (Delorme et al. 2021). KPIC II refers to the series of upgrades started during the first semester of 2022 with the primary goal of doubling the instrument throughput (Jovanovic et al. 2020; Echeverri et al. 2022). The HISPEC is expected to provide Y–K (0.98–2.46 μみゅーm) spectroscopy at a spectral resolution of R > 100,000 (Mawet et al. 2019). The MODHIS is a similar instrument to HISPEC planned for the future TMT. A broader range of exoplanet masses is explored in Section 4 for this latter TMT instrument.

PSIsim includes full budgets of the throughput and thermal background for each instrument, telescope, and the Earth atmosphere. The Strehl ratio is calculated based on an empirically calibrated model of the AO's performance under median seeing conditions for Maunakea. For KPIC I and KPIC II, we assumed Keck AO's current performance with the infrared Pyramid wave front sensor described in Bond et al. (2020). For HISPEC, we assumed extreme-AO performance as predicted for the upcoming HAKA high-density deformable mirror upgrade (W. M. Keck Observatory, private communication). The star is modeled with a PHOENIX model (Husser et al. 2013) and the substellar companion with a BT-Settl atmospheric model grid ²⁴ (Allard et al. 2012a). Table 3 includes the input parameters and the predicted RV precision for these simulations. The simulations include a level of systematics at 1% of the continuum, which is modeled by an additional white Gaussian noise. Otherwise, the estimated RV precision assumes a perfect data reduction.

Table 3. Radial Velocity Precision Simulations of Current and Future Instrumentation for Two Substellar Companions, HR 7672 B and HR 8799 c

Parameters
Star – Phoenix Model	HR 7672	HR 8799
Apparent K mag	4.4 a	5.2 a
Effective temperature (Teff)	6000 K b	7400 K c
Surface gravity ($\mathrm{log}(g)$)	4.5 b	4.5 c
Spin ($v\sin (i);$ km s−1)	5.6 d	49 e
Companion — BT-Settl model	HR 7672 B	HR 8799 c
Mass	73 MJup b	7 MJup f
Apparent K mag	13.0 b	16.1 g
Effective temperature (Teff)	1800 K b	1200 K h
Surface gravity ($\mathrm{log}(g)$)	5.5 b	4.0 h
Spin ($v\sin (i)$)	45 km s−1 b	10 km s−1 i
Separation	072 j	095 j
Telescope and instrument
Airmass	1.2
Water vapor column	1.5 mm
Integration time (tint)	$5\,{\rm{\min }}$
Predicted RV sensitivity (ms−1)
assuming 06–10 seeing
	HR 7672 B	HR 8799 c
Keck/KPIC I (measured)	∼2000 k	∼7000 i
Keck/KPIC I (simulated)	800–1400	3000–5000
Keck/KPIC II	500–800	2000–3000
Keck/HISPEC	∼200	100–200
TMT/MODHIS	30–40	10–20

Notes. Top: representative parameters for the telescope, instrument, star, and companions used in the PSIsim simulations. Bottom: predicted RV sensitivity for values of seeing ranging from 06 to 10.

^a Cutri et al. (2003). ^b Wang et al. (2022a). ^c Wang et al. (2020). ^d Luck (2017). ^e Royer et al. (2007). ^f Wang et al. (2018). ^g Currie et al. (2011). ^h Wang et al. (2018). ⁱ Wang et al. (2021b). ^j http://whereistheplanet.com/ (Wang et al. 2021a). ^k This work.

Download table as:  ASCIITypeset image

The predictions from PSIsim are about a factor 2 more sensitive than existing measurements with KPIC I (see Table 3). This difference can first be explained by uncorrected wave front errors reducing the throughput, both noncommon path aberrations and uncorrected atmospheric turbulence. Then, our current data-analysis framework remains limited in its ability to model KPIC systematics. As explained in Section 2.2, only the redder orders of NIRSpec are being reduced due to strong telluric lines in the bluer orders, and an imperfect Fourier filtering is used to remove the fringing. The gap between the simulations and the measurements should decrease as observing strategies and data-reduction frameworks are improved.

The final expected exomoon sensitivity of the four instruments is shown for the two companions in Figure 4. For a fixed time sampling of the RV series, the minimum detectable mass ratio is approximately proportional to the RV semi-amplitude of the signal, which is also proportional to the RV precision of the instrument, so the improvement for each generation of instrument can be read from the simulated RV precision shown at the bottom of Table 3. These simulations are compared to other detection techniques in the Appendix, specifically astrometric monitoring of the companion or spatially resolving the moon through imaging. We separately discuss the possibility of detecting transiting exomoons using the Rossiter–McLaughlin (RM) effect in Section 5.2.

The mass ratios of the largest gas giant satellites in the solar system are also shown in Figure 4 for comparison. The higher planet masses, M, of directly imaged planets and brown dwarfs compared to the solar system could yield significantly bigger moons, so we also include scaled-up mass ratios, q, according to $q\propto \sqrt{M}$ (Batygin & Morbidelli 2020). While the CPD does scale with the Hill sphere, we do not expect the semimajor axis of satellites to depend on this parameter. Indeed, young moons are thought to migrate toward the planet during their formation due to the interaction with the gas. The migration is stopped at the inner radius of the CPD, which is set by the magnetic field of the planet (Batygin & Morbidelli 2020). In this work, we therefore keep the semimajor axis of the solar system satellites constant. A caveat is that large moons could be susceptible to tidal forces if they form or migrate too close to the planet within the Roche limit. The Roche limit is calculated using the mass–radius relationship from Chen & Kipping (2017) and their associated Python package. ²⁵ However, this relationship does not account for the fact that young objects are likely inflated.

Using KPIC, we derive the most sensitive upper limits on the mass ratio of satellites orbiting a high-contrast substellar companion. We rule out satellites larger than 1%–4% of the mass of the brown dwarf HR 7672 B at separations similar to the Galilean moons. Based on end-to-end simulations, we predict that instruments such as TMT/MODHIS could be two orders of magnitude more sensitive. This would be sufficient to detect moons forming in the CPD of a planet with mass ratios of ∼10⁻⁴, albeit with a substantial investment in observing time. If the satellite-to-planet mass ratio grows as $q\propto \sqrt{M}$, with M the mass of the planet, the Keck/HISPEC should be sensitive to these objects around brown dwarfs. Any detection with HISPEC, or lack thereof, will therefore already be capable of constraining CPD formation models. In order to validate our instrument simulations, we compared them with existing observations. The gap in sensitivity can be explained by imperfections in the data reductions. A continued investment in more accurate data-processing algorithms or observing strategies is therefore required in order to realize these predictions. Planet variability will also be a challenge to overcome using the different timescales and the wavelength dependence of the variability compared to an exomoon signal, for example (Vanderburg et al. 2018). Measuring the variability of substellar companions would in fact be an important result of exomoon surveys to better understand the physics of their atmospheres (Biller 2017).

Binary formation processes favor high-mass ratios so they would be more easily detectable than the smaller satellites forming by accretion in the CPD. The majority of multiplicity surveys for isolated brown dwarfs (Fontanive et al. 2018) or companion brown dwarfs (Burgasser et al. 2005; Lazzoni et al. 2020) have searched for visual companions, leaving the separation regime of <1 au underexplored. Figure 5(b) shows that unresolved binary substellar companion would be detectable with RV precision between 0.1 and 1 km s⁻¹, which is already routinely achieved with KPIC. As an example, the measured dynamical mass of the brown-dwarf companion HD 47127 B suggests that it could be a binary (Bowler et al. 2021), but this specific companion is too faint (K ∼ 18.4) to be a practical target for KPIC.

From the Appendix and Figure 7, we conclude that the different detection techniques are sensitive to distinct regions of the parameter space, and therefore complementary, not unlike exoplanet searches. If exomoons follow the model of solar system gas giant satellites, RV searches could be the most promising approach due to their sensitivity to short-period moons. However, unless the theoretical prediction that bigger planets form even bigger moons holds true, small satellites with mass ratios ∼10⁻⁴ might only be detectable around brown dwarfs.

Zoom In Zoom Out Reset image size

As suggested in Heller & Albrecht (2014), an alternative strategy to look for exomoons around directly imaged planets using RV measurements could be to look for transiting moons through the RM (Gaudi & Winn 2007) effect on the planet. Precise photometric calibration and stability of high-constrast instruments is notoriously difficult (Wang et al. 2022b), so detecting a RM event during a transit could be easier than detecting its photometric counterpart.

An RM event consists of the subsequent masking of a portion of the blue and redshifted areas of the surface of a spinning object, therefore leading to large and very distinct deviations of the measured RV. The amplitude of the RV signal can be hundreds of times larger than the RV semi-amplitude due to the orbital motion of the moon. Its amplitude is proportional to the spin of the planet, which could make it an interesting alternative to detect the smallest moons around rapidly rotating planets and brown dwarfs. Indeed, the RV uncertainties scale with the spin of the object, so detecting the orbital signal of small exomoons could be more challenging.

The Galilean moons have rather small orbital periods from days to weeks. Assuming a random inclination distribution, the transit probability of a moon (P) is given by the ratio of the planet radius (R_p ) and the moon semimajor axis (d_m ), P = R_p /d_m (Borucki & Summers 1984). Therefore, the probability of a transit of a moon at the separation of Io around Jupiter is 1:6, and 1:27 for the farthest Galilean moon Callisto. Assuming a full 8 hr night of observations, we estimate the probability of observing an RM event for Galilean-like moons around a Jupiter-like planet to be around 3% for Io, 1% for Europa, 0.3% for Ganymede, and 0.07% for Callisto. However, the orbital periods of the moons would be even shorter around larger substellar companions, which would increase the probabilities up to 17% for Io, 8% for Europa, 2.6% for Ganymede, and 0.6% for Callisto . The transits would last between ∼2 and 5 hr for the Galilean moons around Jupiter, but they would only last 15–30 minutes for similar moons around HR 7672 B.

As an example, a satellite around HR 7672 B with a mass of 1 M_⊕ (q = 5 × 10⁻⁵) would generate a RM signal of ∼300 m s⁻¹ compared to the ∼0.5 m s⁻¹ generated by the orbital motion (Gaudi & Winn 2007). The amplitude would be ∼5 km s⁻¹ for a Neptune-sized moon. Multiple satellite systems would increase the probability of a detection. The detection probability remains relatively low so RM searches could be carried out in synergy with other science cases such as brown-dwarf variability (Biller 2017). For example, Doppler spectroscopy also favors long observations of rapidly rotating objects, which would make for ideal data sets for exomoon RM searches.

Estimating the occurrence rate of Earth-sized exoplanets in the habitable zone (HZ) of a Sun-like star, called ηいーた_⊕, has been an important goal of exoplanet surveys. While such planets remain challenging to detect, the best estimates of ηいーた_⊕ range between 5% and 50% to date (Gaudi 2021). However, these are not the only Earth-sized objects that could harbor life in the HZ of their stars. Any rocky satellites orbiting HZ gas giant planets could also provide suitable conditions for life. Close-in exomoons can be protected from stellar radiation by the strong magnetic field of Jovian-mass planets (Heller & Zuluaga 2013).

Integrating the distribution of gas giants with an incident flux between 0.3 and 1.5 times the solar irradiance on Earth for an optimistic HZ, or 0.3–1 for a conservative HZ (Kasting & Harman 2013), yields about five to seven giant planets per 100 FGKM-type stars. This is using the giant planet (30–6000 ${M}_{\oplus }\sin i$) occurrence rates derived from the California Legacy Survey as a function of stellar irradiation (see Figure 11 of Fulton et al. 2021). Given that each planet can have multiple satellites, this could represent a significant number of habitable Earth-sized moons that are not accounted for in ηいーた_⊕. The occurence rate of habitable exomoons could be constrained by measuring the population of satellites around more distant directly imaged planets and brown dwarfs.