The TESS-Keck Survey. VII. A Superdense Sub-Neptune Orbiting TOI-1824

We confirm a massive sub-Neptune-sized planet on a P = 22.8 days orbit around the star TOI-1824 (T eff = 5200 K, V = 9.7 mag). TESS first identified TOI-1824 b (formerly TOI-1824.01) as an object of interest in 2020 April after two transits in Sector 22 were matched with a single transit in Sector 21. TOI-1824 was subsequently targeted for ground-based Doppler monitoring with Keck-HIRES and APF-Levy. Using a joint model of the TESS photometry, radial velocities, and Ca ii H and K emission measurements as an activity indicator, we find that TOI-1824 b is an unusually dense sub-Neptune. The planet has a radius R p = 2.63 ± 0.15 R ⊕ and mass M p = 18.5 ± 3.2 M ⊕, implying a bulk density of 5.6 ± 1.4 g cm−3. TOI-1824 b's mass and radius situate it near a small group of “superdense sub-Neptunes” (R p ≲ 3 R ⊕ and M p ≳ 20 M ⊕). While the formation mechanism of superdense sub-Neptunes is a mystery, one possible explanation is the constructive collision of primordial icy cores; such giant impacts would drive atmospheric escape and could help explain these planets' apparent lack of massive envelopes. We discuss TOI-1824 b in the context of these overdense planets, whose unique location in the exoplanet mass–radius plane make them a potentially valuable tracer of planet formation.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Information on the bulk compositions of sub-Neptunes may inform their origins (e.g., "water worlds" are expected to form beyond the ice line; Kuchner 2003;Bitsch et al. 2019;Izidoro et al. 2021), but degeneracies in theoretical models of sub-Neptune atmospheres and interiors make it impossible to infer their composition from mass and radius measurements alone (e.g., Valencia et al. 2007;Adams et al. 2008;Zeng et al. 2019;Aguichine et al. 2021;Rogers et al. 2023).Transmission spectroscopy of sub-Neptune atmospheres may shed light on their interior composition (e.g., Madhusudhan et al. 2023), but interpreting these observations requires precise knowledge of the planet's surface gravity (Batalha et al. 2019).Consequently, the first step toward better understanding the origins of sub-Neptunes as a whole is the measurement of their masses and radii.
Since 2018, NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) has identified more than 7000 transiting exoplanet candidates orbiting bright, nearby stars across the full sky.These detections have enabled large-scale Doppler follow-up efforts, resulting in precise mass measurements for 86 planets smaller than 4 R ⊕ . 29 The TESS-Keck Survey (TKS; Chontos et al. 2022) has contributed to these follow-up efforts through Doppler monitoring of TESS candidate host stars with the Keck High Resolution Echelle Spectrometer (HIRES) and Automated Planet Finder (APF) Levy.The TKS target list was chosen to investigate four science cases: (1) planet bulk compositions, (2) system architectures and dynamics, (3) planet atmospheres, and (4) evolved planetary systems.TOI-1824, an early K dwarf with a sub-Neptune-sized planet candidate in a P = 22.8 days orbit (TOI-1824.01),was selected for radial velocity (RV) observations under science cases 1 and 2.
In this paper, we present the mass of TOI-1824.01(hereafter TOI-1824 b) as measured from our RV observations of the host star with Keck-HIRES and APF-Levy.Our joint analysis of the TESS photometry, the Keck-HIRES and APF-Levy RVs, and the Keck-HIRES Ca II H and K activity indicator measurements (S HK values; Isaacson & Fischer 2010) indicates that TOI-1824 b has R b = 2.63 ± 0.15 R ⊕ and M b = 18.5 ± 3.2 M ⊕ , which implies a bulk density of ρ b = 5.6 ± 1.4 g cm −3 .This mass makes TOI-1824 b unusually dense for its size and rules out models of bulk composition that include a significant (2% by mass) H/He-dominated envelope (Lopez & Fortney 2014).
The paper is organized as follows.In Section 2, we summarize the TESS observations.In Section 3, we describe the high-resolution imaging (HRI) observations.We describe our determination of TOI-1824's properties in Section 4. We detail our Doppler follow-up in Section 5.In Section 6, we explain our light-curve inspection, cleaning, and initial transit fitting.We discuss stellar activity and its possible impact on the RV measurements in Section 7. In Section 8, we perform an initial analysis of the Keck-HIRES and APF-Levy RVs and S HK values, followed by a joint analysis of the photometry, RVs, and S HK values in Section 9.In Section 10, we discuss TOI-1824 b's possible bulk composition and formation scenarios.We also contextualize the planet among the small number of "superdense sub-Neptunes" (e.g., Akana Murphy et al. 2021) and inspect the nature of the RV data sets used to measure their masses.We conclude in Section 11.

TESS Photometry
TOI-1824 was observed with 2 minutes cadence in TESS Sectors 21,22,41,48,and 49 between 2020 January 21 and 2022 March 11.30TOI-1824 was also observed by TESS in Sector 15 with 30 minutes cadence, but we do not include these data in our photometric analysis.The Sector 15 data contain a single transit of TOI-1824 b.However, as demonstrated for Kepler observations (Borucki et al. 2010), long-cadence photometry does not generally have the time resolution to accurately constrain a transit's impact parameter (b), which can lead to a biased R p /R * measurement (Petigura 2020).To avoid this complication and to simplify the analysis, we elected to only use the 2 minute cadence TESS data.
The photometry was processed by the TESS Science Processing Operations Center pipeline (SPOC; Jenkins et al. 2016).For TOI-1824, there are no individual sources from Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2016Collaboration et al. , 2022) ) within 20″ that cause >1% dilution, nor does the combined flux of all DR3 sources within that radius cause >1% dilution.Furthermore, the SPOC data products we use are already corrected for dilution from Gaia sources per Gaia Data Release 2 (Gaia Collaboration et al. 2018).In Section 3, we present HRI observations that rule out significant dilution from unresolved companions. 31All of the TESS data used in this paper can be found in the Mikulski Archive for Space Telescopes: 10.17909/sp43-n767.
We note that there is another star, TIC 142387022, in the same TESS pixel as TOI-1824, with a sky-projected separation of about 7″.TIC 142387022 is an M dwarf with T eff = 3800 K (Paegert et al. 2021) and is fainter than TOI-1824 by 4.1 mag in the TESS bandpass.Following Equation (7) in Ciardi et al. (2015), the dilution from TIC 142387022 has less than a 1% impact on TOI-1824 b's measured transit depth.This is a negligible effect relative to other sources of error that add to the uncertainty on the planet radius measurement (e.g., the error on the stellar radius measurement).Using parallax and propermotion measurements from Gaia Early Data Release 3 (Gaia Collaboration et al. 2021) and the likelihood ratio test from Oh et al. (2017), Behmard et al. (2022) find that TIC 142387022 is almost certainly gravitationally bound to TOI-1824.Initially, comments from the TESS Follow-up Observing Program (TFOP) suggested that there was a second, even fainter star in the TESS pixel (ΔT = 6.3 mag; TIC 142387024), but this purported source was later determined to be a "phantom" artifact (Paegert et al. 2021).

Palomar-PHARO Observations and Data Reduction
Palomar Observatory HRI observations of TOI-1824 were made with the PHARO instrument on the 5.1 m Hale telescope (Hayward et al. 2001).Palomar-PHARO has a pixel scale of 0 025 pixel −1 for a total field of view of about 25″.Observations of TOI-1824 were taken on 2021 June 22 in both the Brγ (λ 0 = 2.169 μm and Δλ = 0.0323 μm) and the H-continuum (λ 0 = 1.668 μm and Δλ = 0.018 μm) narrowband filters.Observations were acquired using the natural guide 29 For planets discovered by TESS with R p < 4 R ⊕ and M p /s > 5 Mp , according to the NASA Exoplanet Archive as accessed on 2024 February 1.
star adaptive optics (AO) system P3K (Dekany et al. 2013) in the standard five-point quincunx dither pattern with steps of 5″.Each dither position was observed three times, with 0 5 positional offsets between each observation, for a total of 15 frames.In both filters, each frame had an integration time of 1.4 s, amounting to a total on-source time of 21 s.The Brγ image has a resolution of 0 09 and a 5σ contrast at 0 5 of Δ7.9 mag.The H-continuum image has a resolution of 0 08 and a 5σ contrast at 0 5 of Δ7.0 mag.
The Palomar-PHARO data were reduced as follows.First, the science frames were flat-fielded and sky-subtracted.The flat fields were generated from a median average of dark subtracted flats taken on-sky.The flats were normalized such that the median value of the flats is unity.The sky frames were generated from the median average of the dithered science frames; each science image was then sky-subtracted and flatfielded.The reduced science frames were combined into a single coadded image using an intrapixel interpolation that conserves flux, shifts the individual dithered frames by the appropriate fractional pixels, and median-coadds the frames.The final resolutions of the combined dithers were determined from the full width at half-maximum (FWHM) of the pointspread functions (PSFs) in the corresponding filter.
The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20°at separations of integer multiples of the central source's FWHM (Furlan et al. 2017).The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance.The resulting brightness of the injected sources relative to each target set the contrast limits at that injection location.The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation.The uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance.The final sensitivity curve for the Palomar-PHARO Brγ image is shown in Figure 1.No stellar companions were detected within 1″.

TIC 142387022: A Stellar Companion to TOI-1824
In Brγ, TIC 142387022 has an on-sky separation from TOI-1824 of 7 2 at a position angle of 161°E of N. Using Gaia EDR3 parallax and proper-motion measurements and the likelihood ratio test from Oh et al. (2017), Behmard et al. (2022) find that TOI-1824 and TIC 142387022 are comoving. 32At TOI-1824's distance of 59 pc, an on-sky separation of 7 2 translates to a sky-projected separation of about 420 au.The TESS Input Catalog version 8.2 (TICv8.2;Paegert et al. 2021) lists TIC 142387022 as an M dwarf with T eff = 3800 K and = g log 4.9 10 .

Determination of Stellar Properties
We used HIRES (Vogt et al. 1994) on the 10 m Keck I telescope at the W. M. Keck Observatory on Maunakea to obtain an iodine-free spectrum of TOI-1824 at high resolution and signal-to-noise ratio (S/N; a "template" spectrum).The template spectrum was collected with the B3 decker (14″ × 0 574, R = 60,000), the length of which allows for effective sky subtraction.The template was obtained at an air mass of 1.57 on UT 2020 July 11 with an exposure time of 590 s, resulting in an S/N of 215 pixel −1 at 5500 Å.The template was obtained while the Moon was below the horizon.This observation was used to create a deconvolved stellar spectral template (DSST) for TOI-1824.Triple-shot exposures of rapidly rotating B stars were taken with the iodine cell in the light path immediately before and after the high-resolution template was collected to precisely constrain the instrumental PSF.The data collection and reduction followed the methods of the California Planet Search (CPS) as described in Howard et al. (2010).
We performed an initial stellar characterization using SpecMatch-Emp (Yee et al. 2017)  To estimate the posteriors of TOI-1824's fundamental parameters, we used isoclassify (Huber et al. 2017;Berger et al. 2020) in grid mode.isoclassify infers marginal posteriors for stellar properties by integrating over a grid of MIST isochrones (Choi et al. 2016).To inform the isoclassify analysis, we input priors stemming from our SpecMatch-Emp results, TOI-1824's parallax from Gaia DR3, and Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) JHK s magnitudes.Following Tayar et al. (2022), to account for model-dependent systematic uncertainties, we inflated the errors on TOI-1824's mass and radius by adding an additional 5% and 4% uncertainty, respectively, in quadrature with the measurement error reported by isoclassify.The final stellar parameters are summarized in Table 7.Note that in the postage stamp image (inset), the field of view shown is too small to include the stellar companion, TIC 142387022 (which has a separation of about 7″).North is up, and east is left.

Keck-HIRES
We obtained 58 high-resolution spectra of TOI-1824 with Keck-HIRES between UT 2020 June 2 and UT 2022 July 31.The observations have a median exposure time of 307 s and a median S/N of 155 pixel −1 at 5500 Å.The majority of the observations were taken using the B5 decker (3 5 × 0 861, R = 45,000), while five spectra were collected with the C2 decker (14″ × 0 861, R = 45,000).The longer length of the C2 decker is more useful when sky subtraction is important (generally for stars with V > 10 mag), but for TOI-1824, there is little difference between observations with B5 versus C2.
The Keck-HIRES spectra were obtained with a warm (50°C) cell of molecular iodine in the light path (Butler et al. 1996), and RVs were computed following the methods of Howard et al. (2010).In brief, the superposition of the iodine absorption lines on the stellar spectrum provides both a fiducial wavelength solution and a precise, observation-specific characterization of the instrument's PSF.Each iodine-in spectrum was then modeled as the sum of the DSST and the instrumental PSF convolved with an "infinite-resolution" laboratory spectrum of the molecular iodine.For each spectrum, we also measure the Ca II H and K emission strength (S HK values; Isaacson & Fischer 2010;H. Isaacson et al. 2024, in preparation).

APF-Levy
We also obtained spectra of TOI-1824 using the Levy spectrograph mounted on the 2.4 m APF telescope (Vogt et al. 2014) at Lick Observatory.We obtained 275 high-resolution spectra of TOI-1824 with APF-Levy between UT 2020 April 28 and UT 2023 February 1.The observations had a median exposure time of 1800 s and a median S/N of approximately 57 pixel −1 at 5500 Å.The observations were taken with the W decker (1″ × 3″, R = 95,000).
The reduction pipeline used to compute RVs from APF-Levy spectra follows the methods of Howard et al. (2010).As with our Keck-HIRES observations, spectra were obtained with a warm cell of molecular iodine in the light path.We also obtained three iodine-free template spectra with APF-Levy on UT 2021 November 2. Each iodine-free spectrum had an S/N of 31 pixel −1 at 5500 Å.Although we obtained iodine-free observations with APF-Levy, we found that computing the APF-Levy RVs using the Keck-HIRES template resulted in systematically smaller scatter and measurement uncertainties than when using the APF-Levy template.Keck-HIRES templates have been shown to serve as effective replacements for APF-Levy templates in the CPS Doppler reduction pipeline (e.g., Dai et al. 2020;MacDougall et al. 2021;Dalba et al. 2022;Lubin et al. 2022).They also provide an efficient alternative to the long exposure that would otherwise be required to achieve a similar S/N on an iodine-free APF-Levy observation.
To avoid using low-quality data, we set a minimum threshold of 1000 counts at 5500 Å for our APF-Levy spectra (typical TOI-1824 APF-Levy exposures have about 3300 counts at 5500 Å).This threshold removed several outlier RVs that all had >9 m s −1 measurement uncertainty.Table 1 summarizes the RV observations, and Table 2 contains the RV measurements.

Light-curve Inspection, Cleaning, and Initial Transit Fitting
We followed the methods of Akana Murphy et al. (2023) to inspect, clean, and fit an initial transit model to the TESS light curve.Using lightkurve (Lightkurve Collaboration et al. 2018), we downloaded all of the 2 minute cadence Simple Aperture Photometry (SAP; Twicken et al. 2010;Morris et al. 2020) data for TOI-1824, excluding data with NaN values or data quality flags.We then normalized the data on a sector-bysector basis.We chose to use the SAP light curve for our analysis instead of the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) light curve because TOI-1824 exhibits a strong (∼1% peak-to-peak variability) spot modulation signal in the TESS photometry.The PDC algorithm is known to suppress stellar activity signals with P  10 days.We found that applying the following analysis to TOI-1824's PDCSAP light curve resulted in measured transit parameters that were consistent with the SAP analysis, but the stellar variability signal was obscured by the PDC algorithm.The Gaussian process (GP; e.g., Rasmussen & Williams 2006) that we use to flatten the light curve (see Section 6.2) is flexible enough to  account for both the stellar activity signal and any spacecraft systematics that are present in the SAP light curve.

Transit Search
We conducted a blind transit search of the TESS data before applying any models to the photometry.We searched for transits in the TESS SAP light curve using the box leastsquares method (BLS; Kovács et al. 2002), as implemented in astropy (Astropy Collaboration et al. 2018).We recovered a 22.81 day signal with S/N = 13, which corresponds to the period of TOI-1824 b as reported in the TESS object of interest (TOI) catalog (Guerrero et al. 2021).After masking the five transits associated with this signal in the 2 minute TESS data, we performed another BLS search but found no indication of additional transit-like events.

Light-curve Cleaning and Initial Transit Fitting
After identifying the planet transits, we cleaned the SAP light curve with an outlier rejection scheme.For each sector, we smoothed the normalized TESS SAP data in bins of 0.3 days with a cubic Savitzky-Golay filter (Savitzky & Golay 1964) and iteratively removed out-of-transit (OoT), 5σ outliers until convergence.We used the SPOC-reported orbital period, time of transit, and transit duration for TOI-1824 b to mask the planet transits.Figure 2 illustrates the results of the Savitzky-Golay filtering for TOI-1824's Sector 21 SAP data.The Savitzky-Golay filtering removed a total of 30 outliers across all sectors.The iterative filtering routine converged in three iterations.
With the light curve cleaned, we simultaneously fit a transit model and a GP to the cleaned data.The GP was used to remove low-frequency stellar variability and instrumental systematics.The transit model uses the quadratic limbdarkening law (Kipping 2013) from starry (Luger et al. 2019), and the GP was constructed in celerite2 (Foreman-Mackey 2018).Following Kipping (2013), the limb-darkening coefficients are parameterized as , where u 1 and u 2 are the usual quadratic limb-darkening coefficients.The transit model is parameterized using P ln , T c , R R ln p * , b, and T ln dur .
The parameters and priors of this initial, photometry-only model are generally the same as for our joint model of the photometry and RVs (see Section 9.1 and Table 5).The main difference between the two is that the joint model does not assume a circular orbit-rather, it explicitly uses w e cos and w e sin to parameterize the transit instead of T ln dur .For this initial, photometry-only model, we placed a broad Gaussian prior on T ln dur , the center of which was the logarithm of the transit duration as reported in the TOI catalog when accessed on 2022 October 4 and whose width was ln 10 days.This initial transit model also assumed no information about the stellar mass, since assuming a circular orbit and fitting in terms of T dur implies a stellar density.
The kernel of the GP used to flatten the light curve is in the form of an overdamped stochastic harmonic oscillator (SHO).The power spectral density (PSD) of the SHO kernel can be written as where ω f is the angular frequency, ω 0 is the undamped fundamental angular frequency, S 0 is the power at ω 0 , and Q is the quality factor of the oscillation.Following the reparameterization for the SHO PSD from the celerite2 documentation, we define where ρ is interpreted as the undamped fundamental period of the oscillator, τ is the characteristic timescale of the damping, and η scales the amplitude of the GP (i.e., η 2 populates the diagonal of the GP covariance matrix).Rewriting Equation (1) in terms of ρ, τ, and η, we have The parameters and priors for this GP kernel are the same as used for the GP that flattens the light curve in the joint model (see Section 9.1 and Table 5).We note that a lower bound of 1 day was placed on ρ and τ to prevent the GP from overfitting the transits (see Figure 6).We also note that this GP component is not meant to serve as a physically interpretable model of the stellar activity signal (which we investigate in later sections of the paper) but rather as a flexible regressor that can effectively flatten the light curve while preserving the intrinsic shape of each planet transit.In this regard, the GP component serves a similar purpose to other methods such as flattening light curves with a basis spline (e.g., Vanderburg & Johnson 2014).
We iteratively fit this model to the data and then identified and removed 7σ outliers in the residuals.From two iterations of fitting and outlier removal, we removed a total of 32 outliers, after which no other outliers were identified.The maximum

Search for Transit Timing Variations
In the interest of completeness, we searched for transit timing variations (TTVs; e.g., Hadden & Lithwick 2017) across the five transits found in the 2 minute cadence TESS data.We used the best-fitting transit time and orbital period for TOI-1824 b from our initial photometry-only transit model (above) as a reference for the expected transit times.We performed a MAP fit of the photometry that was analogous to the initial transit model, but now P ln and T c were replaced with free parameters for the midpoint of each individual transit.We placed a Gaussian prior on each of the observed transit times centered at the expected time with a width of 1 day.
To produce posterior estimates for the difference between the observed and expected transit times (O − C), we used No-U-Turn sampling (NUTS; Hoffman & Gelman 2014), an adaptive form of Hamiltonian Monte Carlo (HMC; Duane et al. 1987;Neal 2011), as implemented with exoplanet and pymc3 (Salvatier et al. 2016a).We employed the same method of posterior estimation for our joint model of the photometry, RVs, and activity indices (Section 9), so we defer a slightly more detailed description of the sampling to Section 9.4.
The results of our TTV analysis are shown in Figure 3. Any differences in the observed and expected transit times for TOI-1824 b are consistent with zero.Given the lack of obvious evidence for TTVs in the TESS data, we excluded any accounting for them in our subsequent analysis.

Photometric Variability
Stellar activity can confound efforts to measure the masses of small planets via Doppler monitoring, especially when the stellar rotation period-or one of its harmonics-is close to a planet's orbital period (e.g., Vanderburg et al. 2016).It is clear from the TESS SAP light curve that TOI-1824 is spotted-the rotation of these spots in and out of view induces as much as a 1% modulation in the star's observed flux.A generalized Lomb-Scargle (GLS;Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009) periodogram of the OoT TESS SAP light curve suggests a stellar rotation period of P = 26.8days, though the width of this peak in the GLS power spectrum spans P ≈ 20-40 days (see the top panel of Figure 10).
Scattered light from the Earth and Moon is known to contaminate TESS observations (e.g., Dalba et al. 2020).While the GLS periodogram of the 2 minute cadence OoT SAP photometry suggests a stellar rotation period for TOI-1824 of P = 26.8days, this value is suspiciously close to the length of a TESS sector (27 days, which comprises two TESS orbits).A natural question then follows: is the periodicity seen in the SAP light curve related to spacecraft systematics or the stellar rotation period?To address this issue, we computed a systematics-insensitive periodogram (SIP) of TOI-1824's 2 minute cadence SAP data using the TESS-SIP tool (Hedges et al. 2020).TESS-SIP computes a Lomb-Scargle periodogram while simultaneously detrending the TESS SAP photometry using principal component analysis to mitigate contamination from systematic signals.The approach in Hedges et al. (2020) is similar to the methods used by Angus et al. (2016) to detrend K2 photometry (Howell et al. 2014).TESS-SIP also estimates the contribution to the SAP flux from pixels outside of the SAP aperture.The SIP and background flux periodogram for the TOI-1824 SAP light curve are shown in Figure 4. Just like the GLS periodogram of the SAP data, the SIP shows a peak near P = 26.8days.There is no such peak in the periodogram of the background flux.These results suggest that the modulation seen in the TESS SAP light curve is indeed astrophysical.
In addition to using a periodogram-based approach to estimate the rotation period of TOI-1824, we also estimated the rotation signal in the OoT TESS SAP photometry by calculating the data's autocorrelation function (ACF) with SpinSpotter (Holcomb et al. 2022).To mitigate the influence of the TESS window function, we applied SpinSpotter in three different cases: (1) to the Sector 21 and 22 OoT data only, (2) to the Sector 48 and 49 OoT data only, and (3) to all of the OoT TESS 2 minute cadence data (Sectors 21,22,41,48,and 49).The highest peak in the ACF and its half-width at half-maximum was P = 21.5 ± 3.1 days,  P = 8.1 ± 1.4 days, and P = 20.0 ± 4.0 days for cases 1, 2, and 3, respectively.We note that in case 2, the peak at P = 8.1 days is accompanied by clear peaks at P = 12.2 days and P = 24.3days.Given the GLS periodogram's preference for a stellar rotation period in the range of P ≈ 20-40 days and the results of cases 1 and 3, the period identified by SpinSpotter when applied to the OoT Sector 48 and 49 data may be the second harmonic of the true stellar rotation period.The ACF may prefer a harmonic of the true rotation period in this case given the larger data gaps in Sectors 48 and 49 compared to Sectors 21 and 22.How might this photometric variability influence the RV time series?Equation (2) from Vanderburg et al. (2016) provides a means of estimating the amplitude of RV variations induced by starspots as a function of photometric variability, assuming that the photometric and spectroscopic observations are taken in roughly the same bandpass.This is true for the TESS photometry and the Keck-HIRES and APF-Levy RVs, which are all measured in the optical.The relationship from Vanderburg et al. (2016) states: where RV pp is the peak-to-peak RV variation caused by starspots, F pp is the peak-to-peak flux variation, and v i sin * is the sky-projected stellar rotational velocity.Applying Spec-Match-Syn (Petigura et al. 2017) to our high-resolution, high-S/N, iodine-free template of TOI-1824, we find v i sin * < 2 km s −1 .We measure F pp following the methods of McQuillan et al. (2014), as described below.
For Kepler targets, McQuillan et al. (2014) defined the periodic photometric amplitude as the range between the 5th and 95th percentiles of the median-divided, unity-subtracted, 10 hr boxcar-smoothed light curves.We applied the same procedure to the TESS SAP light curve of TOI-1824 to estimate the amplitude of the photometric spot modulation, F pp .Since the amplitude of the modulation seems to dampen in Sectors 41, 48, and 49 compared to the first two sectors of observations, we measured this amplitude twice, once for the Sector 21 and 22 data and then again for the Sector 41, 48, and 49 data.We find F pp = 11.5 parts per thousand (ppt) for Sectors 21 and 22.For Sectors 41, 48, and 49, we find F pp = 4.6 ppt.Combining these results with our upper limit on v i sin * , Equation (6) suggests that RV pp  23 m s −1 for Sectors 21 and 22, and RV pp  9.2 m s −1 for Sectors 41, 48, and 49.

Ca II H and K Emission
In addition to its spot-induced photometric modulation, near-UV emission indicates that TOI-1824 is chromospherically active.With each of our Keck-HIRES spectra, we measure Ca II H and K emission via ¢ R log 10 HK (Middelkoop 1982;Noyes et al. 1984) and S HK indices (Isaacson & Fischer 2010).From these measurements, we find that TOI-1824 has ¢ R log 10 HK = −4.58± 0.05 and S HK = 0.36 ± 0.01.For reference, over its magnetic cycle, the Sun oscillates between ¢ R log 10 HK = −5.05 at the solar minimum and ¢ R log 10 HK = -4.84 at the solar maximum (Meunier et al. 2010).Applying the activity-rotation relation from Noyes et al. (1984) to TOI-1824, B − V = 0.83 and ¢ R log 10 HK = −4.58± 0.05 imply a rotation period of P rot ≈ 21 ± 3 days.We note that the uncertainty on this estimate of P rot reflects only the uncertainty on our measurement of ¢ R log 10 HK .
Examining the GLS periodograms of the Keck-HIRES ¢ R log 10 HK and S HK measurements, we find that in each case, the highest peak in the range P = 2-100 days occurs at P = 12.1 days.For both the ¢ R log 10 HK and the S HK periodograms, the peak at P = 12.1 days reaches a false-alarm probability (FAP) level of 5% (Baluev 2008).The GLS periodogram of the Keck-HIRES S HK values can be seen in the second panel from the bottom of Figure 10.
To investigate the potential impact of this activity on the RV measurements, we examined the linear correlation between the Keck-HIRES S HK values and the Keck-HIRES RVs with the Keplerian signal of TOI-1824 b removed (where the model of TOI-1824 b's RV orbit used to produce the residuals was taken from our RV signal detection results, as described in Section 8.1).We calculated the Spearman rank-order correlation coefficient (r Spearman ) and the Pearson correlation coefficient (r Pearson ) for the S HK values and Keck-HIRES RV residuals (e.g., Press et al. 1992).The data appear to have a slight linear correlation (r Spearman = 0.24 and r Pearson = 0.30), indicating that stellar activity may be contaminating the RVs (see Figure 5).
In summary, given the OoT TESS SAP photometry and TOI-1824's Ca II H and K emission, it seems reasonable that the stellar rotation is in the neighborhood of P rot ≈ 24 days.Our GLS periodogram analysis of the OoT TESS SAP photometry suggests P rot ≈ 20-40 days, with the highest peak in the periodogram occurring at P = 26.8days.Our ACF-based analysis of the photometry suggests P rot = 21.5 ± 3.1 days, P rot = 8.1 ± 1.4 days (which may represent the second harmonic of the true period), and P rot = 20.0 ± 4.0 days for different partitions of the data (see Section 7.1).The activityrotation relation from Noyes et al. (1984) implies P rot = 21 ± 3 days.GLS periodograms of the Keck-HIRES ¢ R log 10 HK and S HK measurements both show peaks rising above the 5% FAP level at P = 12.1 days, which could represent the first harmonic of the stellar rotation period.Therefore, including a component of the RV model to address stellar activity, such as a GP (e.g., Robertson et al. 2013;Haywood et al. 2014;Grunblatt et al. 2015;Kosiarek & Crossfield 2020), seems warranted given TOI-1824's high levels of Ca II H and K emission, the proximity of the suspected stellar rotation signal to the orbital period of the transiting planet (P = 22.8 days), and the slight linear correlation of the RV residuals and the S HK values.

Signal Detection
Before modeling the spectroscopic orbit of TOI-1824 b, we first attempted to independently identify the planetary signal in the RV time series.We also searched for additional signals such as those of nontransiting planets and/or the manifestation of stellar activity.We used RVSearch (Rosenthal et al. 2021) to look for periodic signals in our RV data.RVSearch iteratively fits Keplerian orbits to RV data using RadVel (Fulton et al. 2018) and subsequently searches for significant periodicity in the residuals.RVSearch uses an empirical FAP threshold of 0.1%, which is computed using the Bayesian information criterion (BIC; Schwarz 1978): where n is the number of observations, k is the number of free parameters in the model, and  is the maximum of the likelihood function with respect to the model parameters.
A blind search of the RVs with RVSearch identifies a Keplerian signal at the FAP < 0.1% level that matches the period of TOI-1824 b (P = 22.8 days) and has K = 4.5 m s −1 with e = 0.2.After removing this signal, RVSearch also identifies a signal at P = 12.3 days.The nature of this second signal is not immediately clear.However, based on our discussion in Section 7, P = 12.3 days could feasibly be the first harmonic of the stellar rotation period, meaning that interpreting this additional signal as a nontransiting planet should be treated with skepticism.If the P = 12.3 day signal were indeed a nontransiting planet, assuming a circular orbit, b > 1 would imply that i p < 87°.8.For reference, we find that TOI-1824 b has i p = 89.4± 0.2.

Keplerian-only Modeling
After searching for signals in the RV time series with RVSearch, we modeled the RVs with RadVel.Our RadVel model describes the orbit of TOI-1824 b using five parameters: orbital period (P), time of inferior conjunction (T c ), RV semiamplitude (K ), orbital eccentricity (e), and argument of periastron passage (ω).The latter two variables are encoded using the parameterization w e cos and w e sin .We also include an offset (γ) and a jitter term (σ) for each RV instrument, where the latter is added in quadrature with the pointwise RV measurement error in the likelihood function.Finally, we include a linear RV trend (g  ).RadVel computes the MAP values for these parameters before estimating their posterior distributions with Markov Chain Monte Carlo (MCMC) sampling (as implemented in emcee; Foreman- Mackey et al. 2013).
Table 3 contains our Keplerian-only RadVel model of the RVs.We fixed P and T c to the values from our initial photometric model (Section 6.2) since they are primarily constrained by the transit observations.We placed a uniform prior of  [−1, 1] on w e cos and w e sin in order to enforce e < 1.The RV instrument offsets were fixed to their MAP values for the MCMC sampling.We placed a uniform prior of  [0, 20] m s −1 on the RV instrument jitter.No prior was imposed on the RV semiamplitude.
We used the AICc, the small sample-corrected version of the Akaike information criterion (AIC; Akaike 1974), to compare different RadVel models of the RVs-namely, we compared models to assess the evidence for including orbital eccentricity and/or a linear RV trend.The AICc is defined as where the AIC is defined as and where n, k, and  are the same as in the definition of the BIC (Equation ( 7)).In general, smaller values for these metrics are more favorable.In the following text, we refer to the "AIC" for ease of reference, but the same also applies to the AICc.Let D º -AIC AIC AIC i i min , where AIC i is the AIC of the ith model under consideration and AIC min is the lowest AIC value of all models considered.While there are no hard-and-fast rules that dictate model selection with the BIC and/or AIC, Note.BTJD ≡ BJD − 2457000.Burnham & Anderson (2004) provide the following guidelines for interpreting ΔAIC values.
The AICc favors a RadVel model that includes a planet at the orbital period of TOI-1824 b, e and ω, and a linear RV trend.A model that assumes a circular orbit and includes a linear RV trend is nearly indistinguishable (ΔAICc < 2) from the favored model.Eccentric and circular orbit models without linear trends are both slightly disfavored (ΔAICc = 3.0 and ΔAICc = 5.3, respectively).While the AICc slightly favors a model that includes e, ω, and g  over simpler models (i.e., a circular model with no trend), the derived posterior on orbital eccentricity is consistent with zero, and the linear trend is only detected at the 2σ level (see Table 3).For comparison, a Keplerian-only model that assumes a circular orbit and no linear trend returns M p = -+ 16.40 2.66 2.60 M ⊕ , which is consistent with the planet mass measurement from the AICc-favored model shown in Table 3.Because e, ω, and g  are not confidently detected in the RVs, we exclude them from our RV models moving forward.Furthermore, for low-to-moderate-S/N RV orbits, including e is known to systematically inflate K (Shen & Turner 2008).
As mentioned in Section 8.1, in addition to the RV signal of TOI-1824 b, RVSearch identifies a signal in the periodogram of the residuals at P = 12.3 days.The TESS SAP photometry and Keck-HIRES S HK values suggest that this signal is in fact the first harmonic of the stellar rotation period (P rot ≈ 24 days; see Section 7).For completeness, we also fit a RadVel model with two Keplerian signals, one for the orbit of TOI-1824 b and one for the signal at P = 12.3 days, where P and T c for the latter signal were taken from the RVSearch results.We find that our mass measurement for TOI-1824 b is consistent between both the one-and two-Keplerian solutions.

Keplerian plus GP Modeling
Given our discussion of stellar activity, we also explored RadVel models of the data that included a GP to account for the manifestation of this activity in the RV time series.We added a GP to a one-planet Keplerian model of the RVs.We used a "quasiperiodic" GP kernel (e.g., Grunblatt et al. 2015;Kosiarek et al. 2021).For instrument i, the kernel quantifies covariance between data observed as times t and ¢ t as The hyperparameters are η 1−4 : η 1,i represents the amplitude of the covariance for instrument i, η 2 is interpreted as the evolutionary timescale of active stellar regions, η 3 is interpreted as the stellar rotation period, and η 4 is the length scale of the covariance's periodicity.The hyperparameters are shared between instruments except for the amplitudes η 1,i .
The model fitting consisted of first "training" the GP by fitting it to the Keck-HIRES S HK values and then using the resulting posteriors on η 2 , η 3 , and η 4 as numerical priors for a subsequent fit to the RVs.During the training, we also included an instrumental offset and jitter parameter for the Keck-HIRES S HK values.The priors for the GP training were as follows.Both the Keck-HIRES S HK instrument offset and the S HK instrument jitter had a prior of  [0, 1].η 1 had a prior of  [0, 1].We placed a broad Jeffreys prior (Jeffreys 1946) on η 2 of  [30, 10,000] days, where the lower bound was chosen to be slightly larger than the suspected stellar rotation period (Kosiarek & Crossfield 2020).
We also placed a Gaussian prior on η 3 of  (24.6, 5) days, where the center was chosen to be twice the period identified in the RV residuals by RVSearch.Finally, we placed an informed Gaussian prior on η 4 of  (0.5, 0.05) per Haywood et al. (2018).
With the GP training complete, we then fit the APF-Levy and Keck-HIRES RVs with the combined Keplerian plus GP model.The Keplerian and instrument parameter portions of the model were the same as described in Section 8.2.The only differences were that (1) we enforced a circular orbit for TOI-1824 b and (2) we excluded the linear RV trend.We simplified the RV model in this way because both the orbital eccentricity of TOI-1824 b and the linear RV trend were not detected at high significance in our Keplerian-only model.Furthermore, since GPs are highly flexible, we elected to simplify the non-GP portion of the model as much as possible to avoid overfitting.For the GP, we placed a uniform prior of  [0.1, 20] m s −1 on h 1,HIRES and η 1,APF .We used the posteriors from the GP training as numerical priors for η 2 , η 3 , and η 4 .The model parameters, priors, and posterior estimates are summarized in Table 4.The AICc strongly favors a model (ΔAICc > 15) that includes both the GP and the Keplerian signal for TOI-1824 b.This model finds η 3 = 24.9 ± 0.4 days, which agrees with our general discussion of the stellar rotation period in Section 7.

Joint Photometry, RV, and Stellar Activity Modeling
In addition to modeling the photometry, RVs, and stellar activity indicators independently, we also applied a joint model to the data, following the methods of Akana Murphy et al. (2023).The code for reproducing these results is publicly available online (Akana Murphy 2023).We summarize the joint model in Tables 5  and 6.In general, all model parameters had relatively broad priors, save for the stellar mass and radius, whose informed Gaussian priors stemmed from our stellar characterization (see Section 4).The likelihood function of the joint model is the product of the likelihood of the transit, RV, and S HK components, all of which assume Gaussian residuals.

Transits
We parameterize the transit portion of the joint model in terms of P ln , T c , R R ln p * , b, and w e cos and w e sin .As in our initial transit modeling, we use the quadratic limbdarkening law from Kipping (2013).As we did with our Log GP undamped period r ln phot ln days  (ln 10, ln 50)[ln 1, ln 200] F Log GP damping timescale t ln phot ln days  (ln 10, ln 50)[ln 1, ln 200] F Notes. (X, Y) refers to a Gaussian distribution with mean X and standard deviation Y.  (X, Y)[A, B] refers to a bounded Gaussian with mean X, standard deviation Y, and hard bounds at A and B.  [X, Y] refers to a uniform distribution inclusive on the interval X and Y. A: σ phot is treated as a uniform pointwise flux measurement error.s phot refers to the sample standard deviation of the SAP light-curve flux.s RV,i refers to the same for the RVs of instrument i. B: the parameterization º where u 1 and u 2 are the usual quadratic limb-darkening coefficients, follows the prescription by Kipping (2013).C: the bounded Gaussian priors on stellar mass and radius have centers and widths corresponding to our derivation of the stellar parameters in Section 4. D: for some parameter, x, x TOI refers to the value of that parameter as reported in the TOI catalog when accessed on 2022 October 4. E: (ξ 1 , ξ 2 )[0, 1] refers to a uniform distribution over the unit disk (i.e., x x + 1 1 2 2 2  ).VE(e|θ) refers to the mixture distribution from Van Eylen et al. (2019), which is used as a prior on e and whose hyperparameters, θ, are fixed to the posterior medians from that work.F: the hyperparameters of the GP used to flatten the light curve, which has a kernel whose PSD is in the form of an SHO (see Equation (1)).ρ phot and τ phot , the undamped period and damping timescale of the SHO, respectively, are forced to be >1 day to prevent the GP from overfitting low-S/N transits.
initial transit modeling, we fit the transit model simultaneously with a GP using a kernel in the form of an overdamped SHO (Equation ( 5)) in order to flatten the light curve.To prevent the GP from absorbing part of the transit signal, we enforced that the GP's undamped period (ρ) and damping timescale (τ) both be >1 day.We also visually inspected each transit to ensure that the GP's prediction was sufficiently smooth across the transit duration.Figure 6 illustrates the simultaneous transit and GP fitting for TOI-1824 b's transit in Sector 41.

RVs
We used P ln , T c , w e cos and w e sin , and K ln to describe the RV curve of TOI-1824 b's orbit, where all but K ln were shared with the transit model.For each RV instrument, we also included an offset (γ) and an RV jitter term (σ), where the latter is added in quadrature with the pointwise RV measurement errors.We used the AIC to compare joint models with and without a linear RV trend (g  ).We find that the AIC slightly disfavors a joint model that includes a linear RV trend compared to one without it (ΔAIC = 3.2 between the two), so we do not include g  in our adopted solution.

GP Modeling of Stellar Activity
Motivated by TOI-1824's indications of stellar activity, we added a multidimensional GP to our joint model following the methods of Akana Murphy et al. (2023).The GP is fit simultaneously to the RVs and the Keck-HIRES S HK values.Each instrument (APF-Levy RVs, Keck-HIRES RVs, and Keck-HIRES S HK ) is assigned its own GP kernel that shares all hyperparameters with the other kernels, save for the GP amplitude (which we denote with η).In addition to the GP, the Keck-HIRES S HK values are also fit with an offset and jitter term.
The GP kernel is a mixture of three terms, each of which has a PSD in the form of an SHO.The first term is an overdamped oscillator, meant to describe exponentially decaying behavior, such as spot evolution.This term is the same as the kernel of the GP used to flatten the light curve (Equation ( 5)).The only difference is that we fix the quality factor to º Q 1 2 (which effectively sets τ), since this gives the SHO the same PSD as stellar granulation (Harvey 1985;Kallinger et al. 2014).Plugging into and rearranging Equation (5), for instrument i,   we have The second and third terms of the kernel are underdamped SHOs, with fundamental frequencies corresponding to the stellar rotation period and its first harmonic, respectively.For instrument i, the PSDs of these terms can be written as where ω 0 is the undamped fundamental angular frequency, S 0 is the power at ω 0 , and Q is the quality factor of oscillation.The hyperparameters of S P i where η rot,i is the amplitude of + S S P i P i , 2 , rot rot relative to S dec,i , Q 0 is the quality factor minus 1/2 for the oscillator at P rot /2, δQ is the difference between the quality factors of the oscillators at P rot and P rot /2, P rot is the primary period of variability (meant to represent the stellar rotation period), and f is the fractional amplitude of the SHO at P rot /2 relative to the SHO at P rot .
For instrument i, the PSD of the GP kernel can be summarized as the sum of a term describing exponentially decaying behavior (S dec,i ) and a term describing periodic behavior (S rot,i ): The GP parameters and priors are summarized in Table 6.In general, all GP hyperparameters are given relatively broad priors.

Posterior Estimation
We followed the methods of Akana Murphy et al. (2023) to estimate the posteriors of the joint model parameters.The posterior estimation employs NUTS (Hoffman & Gelman 2014), an adaptive form of HMC (Duane et al. 1987;Neal 2011), as implemented with exoplanet and pymc3 (Salvatier et al. 2016a).A NUTS sampler ran eight parallel chains, with each chain taking 8000 tuning steps before drawing 6000 samples.Samples drawn during the tuning period were discarded, similar to how various MCMC methods discard burn-in samples.The chains were concatenated to produce a total of N = 4.8 × 10 4 samples from the marginal posteriors of each model parameter.
We assessed convergence of the HMC sampling through multiple diagnostic statistics.Following Vehtari et al. (2021), we confirmed that the rank-normalized R ˆstatistic for each parameter was <1.01.Also following Vehtari et al. (2021), we verified that the rank-normalized effective sample size in the bulk and tails of each marginal posterior was >400 (typically, we find N eff  1000 in both the bulk and the tails for each model parameter).
The results of our joint model are summarized in Figures 7, 8, 9, and 10.Posterior estimates for stellar and planetary parameters are summarized in Table 7.

Cross-validation to Assess Overfitting
GPs are a common tool for mitigating correlated noise in spectroscopic time series data.However, because of their flexibility, GPs can tend to overfit (Blunt et al. 2023).This issue is of interest in the case of TOI-1824 since the suspected stellar rotation period is close to the orbital period of the planet -though it is unclear in which direction overfitting by the GP would serve to influence the measured value of K, if at all.To assess whether or not the RV component of our adopted joint model is susceptible to overfitting, we performed a series of cross-validation (CV) experiments following the methods described in Blunt et al. (2023).
CV consists of conditioning a model on a randomly selected subset of the data (the "training set") and predicting on the remaining held-out data (the "evaluation set"); a typical training/evaluation split is 80%/20%.If a model is overfitting, the evaluation data residuals should have systematically higher scatter than the training data residuals and/or be offset from zero.For examples, see Figure 5 in Blunt et al. (2023) and Figure 8 in Beard et al. (2024).
We conducted our CV analysis for two different models of the RVs.The first model was a replica of the RV portion of our adopted joint model, except we excluded the GP component.The second model was the RV plus S HK portion of our adopted joint model.For each model, we randomly selected 80% of the APF-Levy data and 80% of the Keck-HIRES data (for the GP-enabled model, the randomly selected Keck-HIRES data included both the RVs and the corresponding S HK values).We fit the model to the training set, made predictions at the time stamps of the evaluation data, and computed the residuals for both.We repeated this process for 100 iterations.
The results of our CV analysis are summarized in Figure 11.The distribution of the residuals of the evaluation data for the GPenabled model is slightly wider than the distribution of the residuals of the training data, but the former is still seemingly normally distributed around zero.Compared to the Keplerian-only model, the GP-enabled model may slightly overfit the training data, but it generally predicts the evaluation data accurately, if not slightly imprecisely.Our CV experiment provides confidence that our adopted joint model of the data is not biased and our reported measurement uncertainties are not grossly underestimated.2021)-hereafter LF14, Z16, and A21 respectively.We followed the methods of Akana Murphy et al. (2023) to estimate the posteriors of the smint model parameters and assess sampling convergence.We placed informed Gaussian priors resulting from our joint model on planet properties where applicable (e.g., on planet mass and instellation flux).
The LF14 grid assumes that a planet is composed of an Earth-like core surrounded by a H/He-dominated envelope (which we choose to have solar metallicity).To interpolate over the LF14 grid, smint takes inputs of planet mass, instellation flux, and system age and determines a H/He envelope mass fraction ( f env ) that best matches the observed planet radius.Assuming an Earth-like interior composition, we find that TOI-1824 b's density is consistent with hosting a - + 1.5 0.6 0.8 % H/Hedominated envelope by mass.
The grid from Z16 models planets as a mixture of liquid H 2 O, high-pressure H 2 O ice, and silicates.Such a composition resembles that of the solar system's icy moons.To interpolate over the Z16 grid, smint fits the planet mass and core H 2 O mass fraction ( f H O 2 ) to match the observed planet radius.
Interpolating over the Z16 grid, we find that TOI-1824 b's mass and radius are consistent with f H O 2 = 45% ± 25%.The Z16 model assumes that all of a planet's H 2 O mass budget is contained in the solid or liquid phases, but this is not necessarily true for warm transiting planets orbiting close to their host stars.Alternatively, the irradiated ocean-world models from A21 describe the planets as an Earth-like core surrounded by a supercritical H 2 O fluid layer and a steam-dominated envelope.To interpolate on the A21 grid, we fit the planet iron core mass fraction ( f Fe ; the fraction of the planet's refractory core that is iron, with the rest of the core being made up of silicates), the planet water mass fraction ( f H O 2 ; which includes both the supercritical fluid and steam envelope components), the irradiation temperature (for which we use T eq assuming zero Bond albedo and full day-night heat redistribution), and the planet mass to match the observed planet radius.We find f Fe = 46% ± 30% and f H O 2 = 35% ± 15%, which implies a silicate mantle mass fraction of about 19%.The inferred value for f H O 2 is smaller when using the A21 grid than compared to the Z16 grid because the water is allowed to exist in an extended steamdominated envelope.

TOI-1824 b Formation and Evolution Scenarios
As described above, TOI-1824 b's mass and radius are consistent with the planet either being an Earth-like core surrounded by a thin (<2% by mass) H/He-dominated envelope or a mixture of water and silicates, where the mixing ratio depends on the phase of the water.While degenerate in mass and radius, the compositions imply distinct formation and evolution mechanisms.
Assuming the Earth-like composition beneath a thin atmosphere, it is difficult to simultaneously explain the planet's large core mass and apparent lack of a substantial envelope, since runaway gas accretion should be triggered once the core reaches about 10 M ⊕ (e.g., Emsenhuber et al. 2021).Furthermore, following the methods of Akana Murphy et al. (2021), it seems unlikely that photoevaporation alone (e.g.,  (Baluev 2008).The posterior median for the periodic hyperparameter of the GP (P rot ) used to model stellar activity in the RVs is marked with the dark red vertical line, and its first harmonic is marked with the light red line.The widths of the vertical lines correspond to the 68% confidence interval on the posterior of P rot and its first harmonic.Our joint model finds P rot = 24.8 ± 0.9 days.We compute the GLS periodogram of the RV window function with a minimum period slightly longer than 1 day to avoid the strong nightly alias that otherwise skews the y-axis scale.Owen & Wu 2013) could have removed the primordial envelope accreted by an 18 M ⊕ core.Using Equation (15) from Lecavelier Des Etangs (2007), we find that at its current orbit, TOI-1824 b could only lose up to 0.1 M ⊕ over 10 Gyr due to extreme ultraviolet radiation from its K dwarf host.TOI-1824 b's mass also seems to rule out formation in situ at its current orbit, since forming an 18 M ⊕ planet at 0.14 au would require a disk mass enhancement factor to the minimum-mass solar nebula of about 40 (Schlichting 2014).Alternatively, a water-rich composition for TOI-1824 b would suggest that the planet formed beyond the ice line and subsequently migrated inward to its present orbit, alleviating the issue of insufficient planet-forming material in the inner disk.This explanation also agrees with the fact that the RVs do not show evidence for a giant planet in the system, which would potentially inhibit water-world formation (e.g., see Figure 1 of Bitsch et al. 2021).Agnostic to the migration mechanism, TOI-1824 b could have formed beyond the ice line via the constructive collisions of primordial icy cores, since water-rich cores are more likely to collide constructively than dry ones, and the collisions are violent enough to drive atmospheric mass loss (Inamdar & Schlichting 2015;Zeng et al. 2019).
Is TOI-1824 b a massive water world?While the existence of water worlds orbiting M dwarf hosts has been debated (Luque & Pallé 2022;Rogers et al. 2023), the picture is even less clear for more massive stars.This is in part an observational biasfor the same planet mass, a young water world orbiting a K dwarf will have less time to migrate to P < 100 days (i.e., be detected in transit) than one orbiting an M dwarf, since the K dwarf's disk will disperse more quickly (Burn et al. 2021).Therefore, only the most massive water worlds will be found on close-in orbits around FGK stars.While the massive-waterworld scenario seems like an appealing explanation for TOI-1824 b's mass and radius, measurements of the planet's atmospheric metallicity are needed before we can break the degeneracy in bulk composition and shed light on the The right column shows the distribution of the rms of the residuals in units of m s −1 for all 100 CV iterations.A vertical black dashed (solid green) line shows the raw rms of the Keck-HIRES (APF-Levy) RVs before applying any models.Bottom row: the same as the top row but now for the model that includes the GP.For this row, the left panel's x-axis units also include the uncertainty in the GP prediction.The spread of the distribution of the residuals for the Keplerian-only model is consistent for both telescopes between the training and evaluation data, suggesting that the model is not overfitting.For the GP-enabled model, the distribution of the residuals for the training data is slightly narrower compared to the distribution for the evaluation data.The rms of the evaluation data residuals is systematically larger than the rms of the training data residuals for the GP-enabled model-though, as seen in the distribution of the residuals, this scatter is typically within 2σ of zero.Takeaway: our CV tests suggest that a GP-enabled model of the RVs might be slightly overfitting the data compared to a Keplerian-only model, but the GP is not leading to model inaccuracy as the residuals of the evaluation data are still centered around zero and typically fall within 1σ-2σ.formation mechanism at work (e.g., Rogers & Seager 2010;Kite et al. 2020).

Are the Mass Measurements of Superdense Sub-Neptunes
Accurate?
Superdense sub-Neptunes (R p  3 R ⊕ and M p  20 M ⊕ ) are a potentially valuable tracer of planet formation because their unique location in the mass-radius plane can a priori rule out certain formation and evolution scenarios.However, are their unusually large mass measurements to be believed?Publication bias is known to favor high-mass planets; for a fixed mass measurement uncertainty, a higher planet mass will yield higher measurement precision and is therefore more likely to be published (e.g., Wolfgang et al. 2016;Burt et al. 2018).Furthermore, with limited telescope resources, if an intermediate analysis reveals a planet detection at high significance, Doppler monitoring of the system may be discontinued so as to allocate exposure time elsewhere.This may result in a concerningly small amount of RV data underlying a highsignificance detection.Below, we consider the reliability of the mass measurements of sub-Neptunes from the NASA Exoplanet Archive with R p < 3 R ⊕ and M p > 12 M ⊕ , where M p ≈ 12 M ⊕ approximately marks the high-mass boundary of the mode of the sub-Neptune mass-radius distribution.
From the sample of Doppler-confirmed planets shown in Figure 12, there are 26 planets, including TOI-1824 b, with M p > 12 M ⊕ and R p < 3 R ⊕ .For reference, there are nine such planets confirmed via TTVs in this mass and radius range, but we do not consider them in the discussion that follows.Of the 26 Doppler-confirmed planets, 18 have been confirmed using >30 RV measurements.For the eight planets confirmed using N RV < 30, we caution that additional follow-up should be carried out to confirm their mass estimates. 33An example of one such planet is K2-182 b (Akana Murphy et al. 2021), a sub-Neptune orbiting an active K dwarf purported to have M p = 20 ± 5 M ⊕ but that was confirmed with only 12 RV measurements.We find that for planets with M p 18 M ⊕ (i.e., about the mass of TOI-1824 b), six out of the 11 planets were confirmed with less than 30 RVs. 34 For planets with M p between 12 M ⊕ and 18 M ⊕ , there are two out of 14 such planets. 35uperdense sub-Neptunes are more likely to be confirmed with small RV data sets (N RV  30) compared to their lowermass counterparts.This may be, at least in part, a symptom of publication bias, whereby inflated Doppler semiamplitude measurements and underestimated errors conspire to produce higher-significance detections.In fact, simulations suggest that Doppler semiamplitudes are preferentially overestimated for small RV sample sizes (Shen & Turner 2008).Additional Doppler follow-up should be undertaken to verify the mass measurements of these planets.

Conclusion
In this work, we reported the discovery and confirmation of a superdense sub-Neptune orbiting the K dwarf star TOI-1824.Using TESS photometry, Keck-HIRES and APF-Levy RVs, and Keck-HIRES Ca II H and K stellar activity indicator measurements, we find that TOI-1824 b has P = 22.80857 ± 0.00006 days, R p = 2.63 ± 0.15 R ⊕ , M p = 18.5 ± 3.2 M ⊕ , and ρ = 5.6 ± 1.4 g cm −3 .TOI-1824 b joins a small population of superdense sub-Neptunes with R p  3 R ⊕ and M p  20 M ⊕ .It remains unclear whether these planets are large, rocky cores surrounded by thin volatile envelopes or massive analogs of the water worlds purported to exist around M dwarfs (Luque & Pallé 2022).Transmission spectroscopy may hold the key to disentangling the composition of these planets, though such observations will be expensive given their high surface gravities.The majority of Doppler-confirmed planets in the literature with R p < 3 R ⊕ and M p 18 M ⊕ were published using N RV < 30.We encourage additional follow-up of these targets to verify the accuracy of their mass measurements.

Figure 1 .
Figure 1.The Palomar-PHARO contrast curve for TOI-1824 in the Brγ filter.Note that in the postage stamp image (inset), the field of view shown is too small to include the stellar companion, TIC 142387022 (which has a separation of about 7″).North is up, and east is left.

Figure 2 .
Figure 2.An example of our Savitzky-Golay filtering procedure for TOI-1824's Sector 21 SAP photometry.The black points are the SAP data, the green line is the data after being smoothed by the Savitzky-Golay filter, and the blue points are the in-transit data, which are not subjected to the outlier rejection.Outliers are marked in red.

Figure 3 .
Figure 3.The results of our exploratory TTV analysis for TOI-1824 b.Differences in the observed and expected transit times are shown as the blue points with 1σ error bars for the five transits in the 2 minute cadence TESS SAP data.We find no evidence of TTVs for TOI-1824 b.Figure 4. Top: the SIP is shown in black.The vertical red dashed line represents the highest peak found in the GLS periodogram of the SAP light curve.The SIP shows a similar peak to the one seen in the GLS periodogram, suggesting that the feature is not necessarily a product of TESS systematics.Bottom: the periodogram of the flux from the pixels outside of the SAP aperture.There is little power in the background flux at the supposed stellar rotation period.

Figure 4 .
Figure 3.The results of our exploratory TTV analysis for TOI-1824 b.Differences in the observed and expected transit times are shown as the blue points with 1σ error bars for the five transits in the 2 minute cadence TESS SAP data.We find no evidence of TTVs for TOI-1824 b.Figure 4. Top: the SIP is shown in black.The vertical red dashed line represents the highest peak found in the GLS periodogram of the SAP light curve.The SIP shows a similar peak to the one seen in the GLS periodogram, suggesting that the feature is not necessarily a product of TESS systematics.Bottom: the periodogram of the flux from the pixels outside of the SAP aperture.There is little power in the background flux at the supposed stellar rotation period.

Figure 5 .
Figure 5. Keck-HIRES RVs-with the Keplerian signal of TOI-1824 b removed-plotted as a function of their associated S HK value.The RV residuals and the S HK values have a slight linear correlation, suggesting that the RV time series is contaminated by stellar activity.
and S HK Shared Hyperparameters for GP Rotation Term Log GP rotation period P ln rot ln days  ( P ln rot, phot , ln 15)[ln 1, ln 60] C Log quality factor of secondary mode lnQ 0 L  (0, 2) L Log quality factor offset between primary and secondary modes dQ ln L  (0, 2) L Fractional amplitude of secondary mode relative to primary mode f L  [0.01, 1] L RV and S HK Shared Hyperparameters for GP Exponential Decay Term Log undamped period of exponential decay term r ln dec ln days  (ln 10, ln 50)[ln 1, ln 100] L Quality factor of exponential decay term Q dec L º1 2 D Note.Notation in this table mirrors that in Table 5. A: s SHK, HIRES is the sample standard deviation of the Keck-HIRES S HK activity indices.B: Inv-Γ refers to the inverse gamma distribution, the parameters of which have been chosen to define the tails of the distribution such that p(x < 0.001) < 0.01 and p(x > 1) < 0.01.This prior helps keep the amplitude of this GP component positive though with a lighter tail near zero as opposed to a gamma distribution.C: P rot, phot is assigned according to the period with the peak power in a GLS periodogram of the TESS photometry.From the TESS SAP photometry, we find P rot, phot = 26.8days.D: fixing º Q 1 2 dec gives this (overdamped) SHO the same PSD as stellar granulation (Harvey 1985; Kallinger et al. 2014).

Figure 6 .
Figure 6.An example of our simultaneous transit and GP fitting for the transit of TOI-1824 b in Sector 21 (i.e., the same sector that is shown in Figure 2).The SAP data are shown in black, and data in 20 minute bins are shown in red.The GP prediction across the transit (plus a small global offset fit to the data) is shown as the green dashed line, and the best-fitting transit model is shown as the blue solid line.We visually inspected each transit to ensure that the GP prediction did not absorb any of the transit signal.

Figure 7 .
Figure 7. Top: the TESS SAP light curve is shown in black, and the GP used to flatten the light curve is overplotted in green.TOI-1824 exhibits peak-to-peak spot modulation of close to 1% during Sectors 21 and 22, though the amplitude of this rotation signal is damped in the subsequent sectors.Blue triangles mark TOI-1824 b transits-the fifth and seventh triangles mark transits that fall in data gaps.Middle: the flattened SAP light curve is shown in black, and the TOI-1824 b transit model is overplotted in blue.Residuals are shown below.Bottom: the flattened SAP light curve folded to TOI-1824 b's orbit is shown in black, and data in 0.5 hr bins are shown in red.The pointwise uncertainty of the (unbinned) data is annotated in the bottom left of the panel.The MAP transit model along with 25 realizations of the model generated from random posterior draws is shown in blue.Residuals are shown below in units of ppt.

Figure 8 .
Figure 8. Top: APF-Levy and Keck-HIRES RVs are shown as the green diamonds and black circles, respectively, where the markers' 1σ error bars represent the quadrature sum of the RV measurement error and an instrument-specific RV jitter term that is fit to the data.The GP plus Keplerian MAP models for both instruments are shown in blue, with error envelopes in corresponding colors representing the 1σ GP prediction uncertainty.Residuals are shown below in units of m s −1 .Bottom: data minus the GP component of the RV model are shown as the green diamonds and black circles.Data from both APF-Levy and Keck-HIRES binned in 0.125 units of phase are shown in red.The Keplerian component of the MAP model along with 25 realizations of the model generated from random posterior draws is shown in blue.Residuals are shown below in units of m s −1 .

Figure 9 .
Figure9.The S HK portion of the joint model, in which GPs, with some shared hyperparameters, are simultaneously fit to the RVs and the Keck-HIRES S HK values.The GP kernel is described in Section 9.3.The posterior median and 68% confidence interval for the GP period is P rot = 24.8 ± 0.9 days.See Figures7 and 8for the photometry and RV portions of the joint model.The GP posterior prediction is shown as the blue line surrounded by a 1σ error envelope.The Keck-HIRES S HK values are shown as the black points with 1σ errors, where the measurement error on each S HK value has been added in quadrature with a jitter term that has been fit to the data.Residuals are shown in the bottom panel.
The sub-Neptune regime of the mass-radius diagram is host to numerous theoretical models of planet composition.
Figure 12 shows TOI-1824 b in the mass-radius diagram alongside the population of confirmed planets and multiple composition curves (Lopez & Fortney 2014; Zeng et al. 2016).TOI-1824 b's density is consistent with that of either a massive "water world" (e.g., Luque & Pallé 2022) or an Earth-like core surrounded by a thin (<2% by mass) H/He-dominated envelope.To infer TOI-1824 b's composition in more detail, we used the Structure Model INTerpolator tool (smint; Piaulet et al. 2021) to interpolate over the theoretical grids of planet composition from Lopez & Fortney (2014), Zeng et al. (2016), and Aguichine et al. (

Figure 10 .
Figure 10.GLS periodograms for TOI-1824 are shown in black.The vertical blue line marks the orbital period of TOI-1824 b (P = 22.8 days).Dashed horizontal lines indicate the 0.1% FAP threshold(Baluev 2008).The posterior median for the periodic hyperparameter of the GP (P rot ) used to model stellar activity in the RVs is marked with the dark red vertical line, and its first harmonic is marked with the light red line.The widths of the vertical lines correspond to the 68% confidence interval on the posterior of P rot and its first harmonic.Our joint model finds P rot = 24.8 ± 0.9 days.We compute the GLS periodogram of the RV window function with a minimum period slightly longer than 1 day to avoid the strong nightly alias that otherwise skews the y-axis scale.

Figure 11 .
Figure 11.Top row: the results of our CV analysis for a Keplerian-only model of the RVs.The left column shows the distribution of the residuals in units of the total model uncertainty (measurement error and telescope jitter added in quadrature) for all 100 CV iterations for the the Keck-HIRES training (solid black) and evaluation data (open hatched gray) and the APF-Levy training (semitransparent green) and evaluation data (open bright green).The right column shows the distribution of the rms of the residuals in units of m s −1 for all 100 CV iterations.A vertical black dashed (solid green) line shows the raw rms of the Keck-HIRES (APF-Levy) RVs before applying any models.Bottom row: the same as the top row but now for the model that includes the GP.For this row, the left panel's x-axis units also include the uncertainty in the GP prediction.The spread of the distribution of the residuals for the Keplerian-only model is consistent for both telescopes between the training and evaluation data, suggesting that the model is not overfitting.For the GP-enabled model, the distribution of the residuals for the training data is slightly narrower compared to the distribution for the evaluation data.The rms of the evaluation data residuals is systematically larger than the rms of the training data residuals for the GP-enabled model-though, as seen in the distribution of the residuals, this scatter is typically within 2σ of zero.Takeaway: our CV tests suggest that a GP-enabled model of the RVs might be slightly overfitting the data compared to a Keplerian-only model, but the GP is not leading to model inaccuracy as the residuals of the evaluation data are still centered around zero and typically fall within 1σ-2σ.

Figure 12 .
Figure 12.The mass-radius diagram in the sub-Neptune regime.Gray dots show planets from the NASA Exoplanet Archive (NEA; accessed on 2023 June 8; NASA Exoplanet Archive 2023) that were confirmed using RVs and have better than 50% and 15% fractional measurement precision in mass and radius, respectively.TOI-1824 b is shown as the green star.Blue triangles represent planets with M p > 12 M ⊕ and R p < 3 R ⊕ whose RV data sets have >30 measurements.Red octagons represent planets in the same mass and radius range whose RV data sets have <30 measurements.We plot various density profiles for reference.The dashed-dotted gray lines represent Earth-like cores surrounded by a 0.1%, 1%, or 2% solar metallicity H/He envelope by mass for a 10 Gyr old planet receiving 10× Earth's instellation (Lopez & Fortney 2014).The dashed blue lines represent either a planet composed of 100% condensed H 2 O or a planet made up of a 1:1 mixture of condensed H 2 O and Earth-like silicates (Zeng et al. 2016).Both are calculated under the assumption of constant specific entropy for a fixed temperature of 700 K at 100 bars.Additional composition curves come from Zeng et al. (2016).TOI-1824 b is emblematic of the degenerate nature of inferring sub-Neptune bulk composition -its mass and radius are consistent with either a 1:1 mixture of water and rock or an Earth-like core surrounded by a 1%-2% H/He envelope by mass.

Table 5
Joint Model of the Photometry and RVs

Table 6
Joint GP Model of the RVs and Keck-HIRES S HK Values