TESS-Keck Survey. V. Twin Sub-Neptunes Transiting the Nearby G Star HD 63935

Nicholas Scarsdale; Joseph M. Akana Murphy; Natalie M. Batalha; Ian J. M. Crossfield; Courtney D. Dressing; Benjamin Fulton; Andrew W. Howard; Daniel Huber; Howard Isaacson; Stephen R. Kane; Erik A. Petigura; Paul Robertson; Arpita Roy; Lauren M. Weiss; Corey Beard; Aida Behmard; Ashley Chontos; Jessie L. Christiansen; David R. Ciardi; Zachary R. Claytor; Karen A. Collins; Kevin I. Collins; Fei Dai; Paul A. Dalba; Diana Dragomir; Tara Fetherolf; Akihiko Fukui; Steven Giacalone; Erica J. Gonzales; Michelle L. Hill; Lea A. Hirsch; Eric L. N. Jensen; Molly R. Kosiarek; Jerome P. de Leon; Jack Lubin; Michael B. Lund; Rafael Luque; Andrew W. Mayo; Teo Močnik; Mayuko Mori; Norio Narita; Grzegorz Nowak; Enric Pallé; Markus Rabus; Lee J. Rosenthal; Ryan A. Rubenzahl; Joshua E. Schlieder; Avi Shporer; Keivan G. Stassun; Joe Twicken; Gavin Wang; Daniel A. Yahalomi; Jon Jenkins; David W. Latham; George R. Ricker; S. Seager; Roland Vanderspek; Joshua N. Winn

doi:10.3847/1538-3881/ac18cb

1. Introduction

Extremely rapid growth in the number of known transiting planets, thanks in large part to the Kepler and K2 missions (e.g., Batalha et al. 2013; Crossfield et al. 2016; Thompson et al. 2018; Kruse et al. 2019), has enabled population-level studies based on bulk properties like orbital period and size. One of the striking features of the population is that the mode of the radius distribution lies near 2.5 R_⊕—a class of planets lacking in the architecture of our own solar system. These super-Earths and/or mini-Neptunes bridge the size domain of terrestrials and ice giants in our solar system.

The occurrence rate distributions of these "bridge" planets exhibit additional features in the period–radius plane. For example, a feature of particular relevance to this work is the so-called "radius cliff," a steep drop in planet occurrence between 2.5 and 4.0 R_⊕ (Borucki et al. 2011; Howard et al. 2012; Fulton et al. 2017) for planets interior to 1 au. Both planets described in this paper fall into this radius regime. Also, the so-called super-Earth (1–1.8 R_⊕) and sub-Neptune (1.8–4 R_⊕) exoplanets are the modes of the known planet radius distribution within these bridge planets (Fulton et al. 2017), a conclusion supported by modeling (e.g., Owen & Wu 2017). Description of the population in relation to the radius cliff and valley has become ubiquitous since their discovery.

As theoretical understandings have matured with the acquisition of mass information, a variety of origins for the radius cliff have been proposed. Kite et al. (2019) proposed atmospheric sequestration into magma as the cause, as larger atmospheres achieve the critical base pressure necessary to dissolve H₂ from the atmosphere into the core. Work is still ongoing to understand how atmospheric observables vary across these features in order to better grasp their underlying physics. Atmospheric characterization of the planets described in this paper could help provide support for theories of the causes of the radius cliff.

Though confirmed multiplanet systems are only a subset of the full planet sample, substantial work has been done to understand the properties of such systems. Weiss et al. (2018) identified the "peas in a pod" phenomenon in the Kepler sample, which describes the fact that planets of a given size are more likely to have neighbors of a similar size than of a random size. The nearly identical planets in the HD 63935 system conform to this trend. Weiss et al. also identified a trend toward denser inner planets, which may be related to photoevaporation.

Mass measurements are resource-intensive to obtain, particularly for the dim host stars common in the Kepler sample, but they are crucial to understanding the population in detail. Early mass measurements by Weiss & Marcy (2014) demonstrated that Kepler planets above ∼1.8 R_⊕ whose stars are sufficiently bright for radial velocity (RV) follow-up appear to retain substantial H/He envelopes. This is in agreement with theoretical predictions for that size regime (Lopez & Fortney 2014). There has since been a great deal of additional work establishing relations between planetary mass and radius, though for any given planetary radius there is a wide spread in masses, implying substantial compositional diversity (e.g., Rogers 2015; Chen & Kipping 2016; Wolfgang et al. 2016; Zeng et al. 2019).

Efforts to understand the sub-Neptune population are being supported by the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2014). Because TESS is observing bright, nearby stars, among the planets it identifies we expect many high-quality atmospheric targets. One of the TESS Level 1 science goals is to ensure that the masses of 50 planets with radii less than 4 R_⊕ are measured. In doing so, we will begin to better understand the processes shaping the underlying distributions of exoplanets. The first TESS catalog paper was recently released, documenting 2241 transiting candidates from the mission (Guerrero et al. 2021). Ongoing follow-up work, including that undertaken by our group, the TESS-Keck Survey (TKS), continues to release well-determined masses for TESS planets.

TKS is a consortium performing precise RV (PRV) follow-up of TESS planet candidates (Dalba et al. 2020; Dai et al. 2020; Rubenzahl et al. 2021; Weiss et al. 2021). One of our group's primary science goals is to measure a diverse set of planet masses at high enough precision to be suitable for atmospheric characterization (Batalha et al. 2019), particularly with the James Webb Space Telescope (JWST), which is what led us to observe HD 63935. One relevant axis of diversity is in host star spectral type. So far, only one sub-Neptune-sized planet around a G star has been the subject of an atmospheric characterization study (HD 3167 c; Mikal-Evans et al. 2021). Because of their relatively small transit signals, such planets represent more challenging targets than similar-sized targets around M or K dwarfs. However, the brightness of the host stars can compensate for this, and a number of compelling targets for atmospheric characterization with JWST around G star hosts have already emerged from TESS (e.g., Gandolfi et al. 2018; Mann et al. 2020; Weiss et al. 2021, Lubin et al. 2021; Turtelboom et al. 2021, in preparation).

HD 63935 b, the subject of this paper, is one of the planets we have identified to be a compelling atmospheric target around a G-type host star. HD 63935 is a bright (V_mag = 8.58) G5 star at a distance of 49 pc. Study of planets around this class of stars is valuable for understanding how differences in host star characteristics shape planet formation. With atmospheric data, we will be able to test theoretical predictions, like those of Lopez & Fortney (2014), that H/He mass fraction is primarily a function of radius. Characterizing these planets will also be valuable for comparative studies with our own solar system. Although we have not discovered any planetary systems that closely resemble our own, the reasons for this are likely observational (Martin & Livio 2015). A better understanding of what makes our solar system unique (or not) is important to the search for life.

In this paper, we present the confirmation of the sub-Neptune planets HD 63935 b and c. Planet b is uniquely well suited to atmospheric characterization, being the second-best target on the radius cliff and the best in its niche of sub-Neptune-sized (2.6 R_⊕ < R_p < 4 R_⊕), moderately irradiated (100 F_⊕ < F_p < 1000 F_⊕) planets around G stars. Planet c is also amenable to atmospheric characterization. We also discuss evidence for a longer-period planetary-mass companion to the two confirmed planets, though we ultimately adopt a two-planet model. In Section 2, we describe the selection algorithm that our consortium's atmospheres working group uses to select high-quality targets. In Section 3 we provide a description of our efforts to characterize the planets' host star. In Section 4, we describe the observations we undertook to confirm this planetary system. In Section 5, we provide details of our analysis and results, and in Section 6 we discuss their implications.

2. TKS: Atmospheric Target Selection

This work is based on data obtained as part of TKS, which performs PRV follow-up of TESS planet candidates using the Keck telescope on Maunakea and the Automated Planet Finder (APF) telescope at Lick Observatory. As part of TKS, we are interested in selecting the TESS planet candidates that would represent the best prospects for atmospheric characterization. TESS has produced (and continues to produce) far too many promising atmospheric targets for one consortium to follow up. Consequently, we attempted to develop an algorithm to prioritize targets to add to our PRV prioritized observing list; see Chontos et al. (2021) for more details about the general selection of TKS targets.

Our algorithm aims to find high-quality atmospheric targets in regions of parameter space mostly bereft of them. We select mostly planets in the sub-Neptune regime, as many highly observable giant planets are already known and terrestrial planets are, with a few exceptions, not accessible with JWST. As a quantification of "underpopulated parameter space," our selection algorithm bins planets in stellar effective temperature, planet radius, and insolation flux. We then select targets that stand out in bins without any characterized planets. A detailed explanation of the algorithm follows.

The inputs to our algorithm are the star and planet properties from the NASA Exoplanet Archive⁴³ and the TESS object of interest (TOI) list.⁴⁴ We subject these lists to certain cuts as well as manual inspections of the data validation reports (Twicken et al. 2018). We exclude TOIs with decl. < −20° and V_mag > 12 for visibility at our facilities and to ensure acceptable signal-to-noise ratios (S/N), respectively. We also cut planets with R_p > 10 R_⊕, as the Jovian population is already reasonably well sampled (for the atmospheres science case only—TKS as a whole does follow up some large planet candidates). We also exclude stars with T_eff > 6500 K.

After our initial culling of the sample, we calculate an estimated mass for each TESS planet candidate based on the fitting formulae of Chen & Kipping (2016) and Weiss & Marcy (2014). With that, we calculate a transmission spectroscopy metric (TSM) value (Kempton et al. 2018). This value is an estimate of the S/N of a planet candidate's atmosphere as observed with the NIRISS instrument on JWST. For our purposes, it serves as a proxy for relative observability.

Our selection algorithm then computes a final parameter, which is equal to the TSM value normalized by the expected exposure time on the High Resolution Echelle Spectrometer (HIRES) that would be required to obtain a mass precision better than 20%, estimated based on Plavchan et al. (2015). This metric is chosen as our ranking parameter in order to select a reasonably large sample of planets. Ranking by the TSM alone results in time spent on M stars with the HIRES, which is suboptimal for a visual light spectrograph. Cool stars are better suited to characterization by instruments at other facilities like MAROON-X (Seifahrt et al. 2018) and CARMENES (Quirrenbach & Consortium 2020), which have more sensitivity in the red. See Figure 1, which plots the most promising atmospheric targets from TESS with their estimated time for measurement of the mass with a precision of 20%, for a visual representation of this.

**Figure 1.** Sub-Neptune TESS planet candidates with estimated TSM values higher than 84 (the value suggested by Kempton et al. (2018) as the cutoff for follow-up in this size regime), colored by their estimated time to a mass precision of 20% with HIRES. The figure also plots known planets with published atmospheric transmission spectrum observations (most from Hubble's Wide Field Camera 3) as black open circles. Population-level studies are more attainable for the sub-Neptunes with G star hosts, since we can obtain high-precision mass measurements on such planets quickly. TOIs included in the TKS atmospheric sample (sub-Neptune-sized planets only) are identified with stars and labeled with their TOI numbers, highlighting the sample's focus on G-type host stars. HD 63935 b is colored green. HD 63935 c is not shown as it was not known as a planet candidate when our target list was finalized.
Download figure:
Standard image High-resolution image

With all relevant metrics calculated, our algorithm divides planets into bins along three axes in parameter space: planet radius, stellar effective temperature, and insolation flux. We use five log-uniform radius bins (which conveniently include edges at 1.7 R_⊕, approximately the location of the radius gap, and at 4 R_⊕, dividing the sub- and super-Neptune populations), five log-uniform insolation flux bins, and three stellar effective temperature bins, resulting in 75 bins total, some of which are unpopulated. TOIs that have higher X metric values than other TOIs and known planets in their bin are prioritized for our PRV observing list. There are enough planet candidates that meet this criterion that we cannot observe them all, and as a result we focus on those with the highest X metrics. Our algorithm permits the removal of candidates deemed observationally unsuitable, e.g., by spectral signatures suggestive of an eclipsing binary, or by substantial stellar activity making continued observations infeasible. This selection process is what leads us to observe HD 63935, the subject of this paper, which remains the highest-ranked target in its bin by the X metric.

3. Host Star Characterization

We obtained high-resolution spectra of HD 63935 (HIP 38374, TIC 453211454; other aliases in Table 1). We used SpecMatch-Syn (Petigura et al. 2017) to obtain the stellar effective temperature, $\mathrm{log}g$ , and metallicity. We then used these values, combined with the luminosity and parallax from Gaia (Prusti et al. 2016; Gaia Collaboration et al. 2021), as priors to obtain the stellar mass, radius, and age using isoclassify (Huber et al. 2017; Berger et al. 2020; Table 3). isoclassify functions by using the input parameters stated above to fit the star to a curve of constant age (isochrone) in the relevant parameter space, allowing the calculation of a radius, mass, and age value. We added a 4% and 5% systematic uncertainty to the stellar radius and mass values, respectively, to account for isochrone grid uncertainty, following Tayar et al. (2020). Note that our derived age does not account for these uncertainties. The values we obtained for mass and radius are consistent with those provided by SpecMatch-Syn, those derived by the Gaia mission (Gaia Collaboration et al. 2021), and those from our spectral energy distribution (SED) fitting. The star is slightly smaller than the Sun (R_* = 0.959 R_⊙, M_* = 0.933 R_⊙), and its isoclassify-derived age suggests that the system is older as well, at 6.8 ± 1.8 Gyr. This is consistent with our nominally low value of $v\sin i={0.24}_{-0.24}^{+1.0}$ km s⁻¹ and low activity indicator logR'_HK = −5.06, suggestive of an old, relatively inactive star. The two sectors of TESS photometry do not provide enough information to obtain a reliable rotation period estimate, though we discuss other methods for obtaining this value in Section 3.2.

Table 1. HD 63935 Identifiers and Gaia Solution

Aliases
HIP ID	38374
TIC ID	453211454
Tycho ID	783-536-1
Gaia EDR3 ID	3145754895088191744

Gaia 6D Solution
R.A.	07^h51^m4204
Decl.	+09°23'1140
Parallax	20.470 ± 0.019 mas
R.A. Proper Motion	−78.696 ± 0.022 mas yr⁻¹
Decl. Proper Motion	−188.512 ± 0.013 mas yr⁻¹
RV	−20.34 ± 0.19 km s⁻¹

Download table as: ASCII Typeset image

Table 2. RVs

Time (BJD)	RV (m s⁻¹)	RV Unc. (m s⁻¹)	Inst.
2,458,733.13862	−6.81	1.29	HIRES
2,458,744.13959	−8.04	1.11	HIRES
2,458,777.06304	0.28	1.16	HIRES
2,458,788.11234	−11.38	1.06	HIRES
2,458,795.00493	−5.68	1.00	HIRES

Note. A sample of the RVs, uncertainties, and instruments for our data on HD 63935. The full table of RV data is available online, which includes the Mount Wilson s-value activity indicators.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

3.1. SED and Activity

As an independent determination of the stellar parameters, we also performed an analysis of the broadband SED of the star together with the Gaia EDR3 parallaxes, in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016) and Stassun et al. (2017, 2018). We obtained the B_T V_T magnitudes from Tycho-2, the Strömgren uvby magnitudes from Paunzen (2015), the JHK_S magnitudes from the Two Micron All Sky Survey, the W1–W4 magnitudes from the Wide-field Infrared Survey Explorer, the GG_BP G_RP magnitudes from Gaia, and the far-UV and near-UV magnitudes from the Galaxy Evolution Explorer. All together, the available photometry spans the full stellar SED over the wavelength range 0.15–22 μm (see Figure 2).

**Figure 2.** SED of HD 63935. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).
Download figure:
Standard image High-resolution image

We performed a fit using Kurucz stellar atmosphere models, with the effective temperature (T_eff), metallicity ([Fe/H]), and surface gravity ( $\mathrm{log}g$ ) adopted from the spectroscopic analysis. The only additional free parameter is the extinction (A_V), which we restricted to the maximum line-of-sight value from the dust maps of Schlegel et al. (1998). The resulting fit is very good (Figure 2) with a reduced χ² of 1.1 and a best-fit A_V = 0.02 ± 0.02. Integrating the (unreddened) model SED gives the bolometric flux at Earth, F_bol = 1.060 ± 0.012 × 10⁻⁸ erg s⁻¹ cm⁻². Taking F_bol and T_eff together with the Gaia EDR3 parallax gives the stellar radius, R_⋆ = 0.967 ± 0.035 R_⊙. In addition, we can use R_⋆ together with the spectroscopic $\mathrm{log}g$ to obtain an empirical mass estimate of M_⋆ = 0.82 ± 0.20 M_⊙, which is consistent with that obtained via the empirical relations of Torres et al. (2010), M_⋆ = 1.02 ± 0.06 M_⊙. These parameters are also consistent with those we derived from isoclassify. Finally, R_⋆ and M_⋆ together yield a mean stellar density of ρ_⋆ = 1.59 ± 0.20 g cm⁻³.

3.2. Predicted Rotation Period from Gyrochronology

Before photometry confirmed the existence of HD 63935 c as a transiting planet, we were interested in obtaining the stellar rotation period in order to rule that out as the source of the RV signal at 21 days. Although TESS Sector 34 photometry has since confirmed that candidate, we include the following gyrochronology analysis in the text as it provides novel information about the host star. We used Markov Chain Monte Carlo (MCMC) within kiauhoku (Claytor et al. 2020) to obtain a posterior probability distribution of stellar parameters for HD 63935. For input we used Gaussian priors based on the spectroscopic effective temperature and metallicity, as well as on the isoclassify-derived age. Assuming a gyrochronological model, we were then able to predict the rotation period. We performed MCMC using two different braking laws on the same grid of stellar models: one ("fastlaunch") uses the magnetic braking law presented by van Saders & Pinsonneault (2013), while the other ("rocrit") uses the stalled braking law of van Saders et al. (2016). With the fastlaunch model, we predicted P_rot = 32.4 ± 6.5 days, while we predicted P_rot = 31.1 ± 4.3 days using the rocrit model. Both of these are consistent with a star somewhat older than the Sun, as implied by our measured R'_HK value and derived via isoclassify, as well as our nondetection of a rotation period in the TESS data.

4. Observations

4.1. TESS Photometry

HD 63935 (TIC 453211454, TOI-509) was observed by the TESS mission in Sector 7 between UT 2019 January 7 and 2019 February 2, and in Sector 34 between UT 2021 January 13 and 2021 February 9. The star was imaged by CCD 4 of Camera 1. The data consist of 33,208 data points with integration times of 2 minutes each. The Science Processing Operations Center (SPOC; Jenkins et al. 2016) processed the data, generated light curves using Simple Aperture Photometry (SAP; Twicken et al. 2010; Morris et al. 2020), and removed known instrumental systematics using the Presearch Data Conditioning SAP (PDCSAP) algorithm (Smith et al. 2012; Stumpe et al. 2012, 2014). The Sector 7 data contained two transits of planet b, and the Sector 34 data contained two additional transits of planet b and two transits of planet c. The transit of planet c that occurred during Sector 7 happened during a gap in the TESS light curve (see Section 5.1). For the analysis described here, we downloaded the PDCSAP flux data from the publicly accessible Mikulski Archive for Space Telescopes (MAST).⁴⁵ The full light curve is plotted in Figure 3 and the phase-folded transits are plotted in Figure 4.

**Figure 3.** The systematic-corrected TESS light curve of HD 63935 in TESS Sectors 7 and 34. Planet b's transits are marked red and planet c's blue. Note that the Sector 7 data gap contains a transit of both planets. The transit-like event near time = 1506 in Sector 7 is an artifact introduced into the light curve during background subtraction.
Download figure:
Standard image High-resolution image

**Figure 4.** Phase-folded transits of HD 63935 b and c. The normalized PDCSAP flux is shown in gray. Our best-fit transit model is shown in red. Residuals from each fit are plotted in the lower panels.
Download figure:
Standard image High-resolution image

4.2. Adaptive Optics Imaging

As part of our standard process for validating transiting exoplanets to assess the possible contamination of bound or unbound companions on the derived planetary radii (Ciardi et al. 2015), we observed TOI-509 with high-resolution near-infrared adaptive optics (AO) imaging at Palomar and Keck Observatories.

The Palomar Observatory observations were made with the PHARO instrument (Hayward et al. 2001) behind the natural guide star AO system P3K (Dekany et al. 2013) on 2019 April 18 UT in a standard five-point quincunx dither pattern with steps of 5'' in the narrowband Br − γ filter (λ_o = 2.1686 μm, Δλ = 0.0326 μm). Each dither position was observed three times, offset from each other by 0 farcs 5 for a total of 15 frames; with an integration time of 1.4 s per frame, the total on-source time was 21 s on-target. PHARO has a pixel scale of 0 farcs 025 per pixel for a total field of view (FOV) of ∼25''. These observations were taken at an airmass of 1.1858.

The Keck Observatory observations were made with the NIRC2 instrument on Keck II behind the natural guide star AO system (Wizinowich et al. 2000) on 2019 March 25 UT in the standard three-point dither pattern that is used with NIRC2 to avoid the left lower quadrant of the detector, which is typically noisier than the other three quadrants. The dither pattern step size was 3'' and was repeated twice, with each dither offset from the previous dither by 0 farcs 5. NIRC2 was used in the narrow-angle mode with a full FOV of ∼10'' and a pixel scale of approximately 0 farcs 0099442 per pixel. The Keck observations were made in the narrowband Br – γ filter (λ_o = 2.1686 μm, Δλ = 0.0326 μm) with an integration time of 0.5 s for a total of 4.5 s on-target. The observations were taken in narrow camera mode with a 1024'' × 1024'' FOV and at an airmass of 1.43.

The AO data were processed and analyzed with a custom set of IDL tools. The science frames were flat-fielded and sky-subtracted. The flat fields were generated from the median average of dark-subtracted flats taken on-sky. The flats were normalized such that their median value was unity. The sky frames were generated from the median average of the 15 dithered science frames; each science image was then sky-subtracted and flat-fielded. The reduced science frames were combined into a single image using an intrapixel interpolation that conserves flux, shifts the individual dithered frames by the appropriate fractional pixels, and median-coadds the frames (Figure 5). The final resolution of the combined dithers was determined from the FWHM of the point-spread function, 0 farcs 094 and 0 farcs 050 for the Palomar and Keck observations, respectively.

**Figure 5.** Final combined full-FOV dithers of the Palomar observations showing no companions within the TESS pixels.
Download figure:
Standard image High-resolution image

4.3. Ground-based Photometry

The TESS pixel scale is ∼21'' per pixel, and photometric apertures typically extend out to roughly 1', which generally results in multiple stars blending in the TESS aperture. An eclipsing binary in one of the nearby blended stars could mimic a transit-like event in the large TESS aperture. We conducted ground-based photometric follow-up observations as part of the TESS Follow-up Observing Program⁴⁶ with much higher spatial resolution to confirm that the transit signal of HD 63935 b is occurring on-target, or in a star so close to HD 63935 that it was not detected by Gaia DR2.

4.3.1. MuSCAT

We observed one partial transit of HD 63935 b on 2019 March 24 from 10:29 to 15:09 in UTC covering the expected egress, with the multicolor simultaneous camera MuSCAT (Narita et al. 2015), which is mounted on the 1.88 m telescope of the Okayama Astronomical Observatory in Okayama, Japan. MuSCAT has three optical channels each equipped with a 1024 × 1024 pixel CCD camera, enabling g-, r-, and z_s-band simultaneous imaging. Each camera has a pixel scale of 0 farcs 358 per pixel, providing an FOV of 6.1'' × 6.1''. The exposure times were 10, 3, and 3 s for the g, r, and z_s bands, respectively.

We performed standard aperture photometry using the custom photometry pipeline described in detail in Fukui et al. (2011). The adopted aperture was 20 pixels or 7'', which excludes any nearby stars as the source of the signal of HD 63935 b. Our precision was not enough to detect the expected 0.09% transit signal on-target in each band, but the data ruled out deep eclipses in all nearby stars within the FOV that are consistent with the transit depth from TESS.

4.3.2. Las Cumbres Observatory Global Telescope

We observed a full transit of HD 63935 b in the Pan-STARRS Y filter (central wavelength 1004 nm) on UTC 2020 November 20 from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network node at McDonald Observatory. We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations. The 4096 × 4096 LCOGT SINISTRO cameras have an image scale of 0 farcs 389 per pixel, resulting in a 26' × 26' FOV. The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018), and photometric data were extracted with AstroImageJ (Collins et al. 2017). The images were focused and had typical stellar point-spread functions with an FWHM of ∼1 farcs 6, and circular apertures with radii ∼7 farcs 8 were used to extract differential photometry. The light curve is presented in Figure 6.

**Figure 6.** Transit of HD 63935 b observed by the LCOGT 1 m telescope at McDonald Observatory on UT 2020 November 20. Top: Observed differential photometry (gray) with data binned in 5 minute intervals shown in color. Bottom: Transit model fit, with detrending on airmass and BJD; the shaded region shows the 68.3% credible interval of the posterior from the MCMC fit (see Section 5.1.1). The error bars in the lower plot include a jitter term determined by the fit.
Download figure:
Standard image High-resolution image

4.4. Ground-based Spectroscopy

4.4.1. Las Cumbres Observatory's NRES

We obtained a spectrum of the target with the automated LCOGT 1 m/NRES optical (380–860) spectrograph (Brown et al. 2013; Siverd et al. 2018), in order to characterize the star and look for signs of a stellar binary system. The observation was done on UT 2019 March 22, at the McDonald Observatory node of the LCOGT Network. We observed the target with two consecutive 20 minute exposures that were processed by the Las Cumbres Observatory (LCO) BANZAI reduction pipeline (McCully et al. 2018) and then stacked together for a spectrum with an effective 40 minute exposure time and an S/N of 73. The reduced spectrum was processed by the SpecMatch-Syn pipeline (Petigura et al. 2017), where spectral and stellar parameters were derived while accounting for the target's distance derived from the Gaia DR2 parallax (Gaia Collaboration et al. 2018). The spectrum did not show evidence for a second set of lines, and the SpecMatch-Syn analysis showed the target is a slowly rotating ( $v\sin i\,={0.24}_{-0.24}^{+1.0}$ km s⁻¹) Sun-like star with an absolute RV (−20.6 ± 0.1 km s⁻¹) consistent with the Gaia DR2 RV (−20.34 ± 0.19 km s⁻¹).

4.4.2. Tillinghast Reflector Echelle Spectrograph

We obtained two spectra, on UT 2019 March 28 and 2019 April 4, using the Tillinghast Reflector Echelle Spectrograph (TRES) on the 1.5 m telescope at the Whipple Observatory on Mount Hopkins in Arizona. TRES is an optical echelle spectrograph with a wavelength range of 385–910 nm and a resolution of R = 44,000 (Fűrész 2008; Mink 2011). The two TRES observations were well separated in the orbital phase (4.21 and 4.59) of the photometric ephemeris and were used to derive relative RVs. Using the strongest observed spectrum as a template, the second spectrum was cross-correlated order-by-order in the wavelength range 426–628 nm. The observed template spectrum was assigned a velocity of zero and the small velocity difference between the two spectra was 13 m s⁻¹, which ruled out a stellar or brown-dwarf companion as the source of the transit-like events. The TRES observations also revealed a spectrum very similar to that of the Sun, with line broadening due to rotation of less than 4 km s⁻¹ and no indication of surface activity, such as emission at Ca ii H and K, thus confirming that this target was well suited for PRV work.

The stellar effective temperature (T_eff), metallicity ([Fe/H]), surface gravity ( $\mathrm{log}g$ ), and rotation ( $v\sin i$ ) were also determined using the Stellar Parameter Classification (SPC) tool on the TRES spectra (Buchhave et al. 2012). SPC cross-correlates observed spectra against a library of synthetic spectra calculated using Kurucz model atmospheres (Kurucz 1993) and performs a multidimensional fit for the stellar parameters that give the highest peak correlation value. These stellar parameter estimates are in 1σ agreement with the results from the SpecMatch-Syn analysis of the PRV observations.

4.4.3. Telescopio Nazionale Galileo/HARPS-N

Between UT 2019 April 2 and 2019 April 29 we collected 11 spectra of TOI-509 with the HARPS-N spectrograph (Cosentino et al. 2012; 383–693 nm, R ≈ 115,000) mounted at the 3.58 m Telescopio Nazionale Galileo (TNG) of Roque de los Muchachos Observatory in La Palma, Spain, under the observing programs CAT19A_162 (PI: Nowak) and CAT19A_96 (PI: Pallé). The exposure time was set to 600–900 s, based on weather conditions and scheduling constraints, leading to an S/N per pixel of 59–119 at 5500 Å. The spectra were extracted using the off-line version of the HARPS-N DRS pipeline (Cosentino et al. 2014), version 3.7.

Doppler measurements and spectral activity indicators (bisector inverse slope, FWHM (CCF_FWHM) contrast (CCF_CTR) of the cross-correlation function, Mount Wilson S-index, and logR_HK index) were obtained using an on-line version of the DRS, the YABI tool,⁴⁷ by cross-correlating the extracted spectra with a G2 mask (Baranne et al. 1996). We measured a logR'_HK value of −4.94. We also used serval ⁴⁸ (Zechmeister et al. 2018) to measure the relative RVs, chromatic index, differential line width, and Hα index. The uncertainties of the relative RVs measured with serval are in the range 0.5–1.4 m s⁻¹, with a mean value of 0.83 m s⁻¹.

4.4.4. Keck HIRES

Between 2019 August and 2021 March, we obtained 51 high-resolution spectra of HD 63935 with the HIRES instrument (Vogt et al. 1994) on the 10 m Keck I telescope at the W. M. Keck Observatory on Maunakea, Hawai'i. We obtained spectra with the C2 decker, which has dimensions of 14'' × 0 farcs 86 and spectral resolution R ≈ 60,000 at 500 nm. The chosen exposure meter setting regulates the S/N at 200 photons per pixel and the resulting median exposure time is 185 s. We obtained RV measurements from the spectra using the method described in Howard et al. (2010). The rms value of our RVs before fitting for any planets was 4.22 m s⁻¹, and the median internal uncertainty was 1.1 m s⁻¹. Our measured logR'_HK value is −5.04, indicating a relatively low-activity star and consistent with the value measured by HARPS-N.

4.4.5. APF-Levy

Between 2019 August and 2021 February, we obtained 100 spectra of HD 63935 with the Levy Spectrograph, a high-resolution slit-fed optical (500–620 nm) echelle spectrograph (Radovan et al. 2010) on the APF telescope at Lick Observatory (Vogt et al. 2014). We observed the star using the W decker, which has dimensions of 3'' × 1'' and R ≈ 114,000 between 374 and 970 nm. Our median exposure time with APF was 1200 s. We acquired one or two observations per night. Nightly observations were binned to improve the RV precision. The rms value of our RVs before fitting for any planets was 8.27 m s⁻¹, and the median uncertainty was 1.8 m s⁻¹. We excluded data points with RV uncertainties >5 m s⁻¹, which resulted in the removal of five points. All such points were clear outliers and had <800 counts on the detector, indicating low data quality.

5. Analysis and Results

5.1. TESS Photometry Analysis

We used juliet (Espinoza et al. 2019) to model the light-curve data available for HD 63935. juliet serves as a wrapper for a variety of existing publicly available tools. Our transit modeling used the functions based on batman (Kreidberg 2015) for transit fitting and on PyMultiNest (Buchner et al. 2014) for the sampling of parameter space. We fit for both planets' periods, crossing times, transit depths, and impact parameters, as well as for two quadratic limb darkening parameters following Kipping (2013), the out-of-transit flux, and a jitter term from TESS. We kept all other parameters fixed, including the mean stellar density, for which we fixed the value to that derived in Section 3. We imposed normal priors centered on the SPOC values for period and crossing time (Jenkins 2002; Jenkins et al. 2010; Li et al. 2019), all with widths of 0.1 days, and constrained the occultation fraction and impact parameter to a uniform range between 0 and 1.

When we began to study this system, the gap in the TESS Sector 7 photometry (Sector 34 data had not been obtained yet) combined with only two transits meant that two periods (∼18 days, the one initially reported by SPOC, and ∼9 days, which we ultimately selected based on RV measurements and which was later confirmed by TESS Sector 34) were possible for planet b. The Sector 7 light curve also provided no indication of the existence of planet c, as its transit fell in the data gap as well. We selected the correct period (9 days) and identified planet c as a candidate based on our RV observations (see Section 5.3). Both of our predictions were validated by the Sector 34 light curve, which confirmed the 9 day period as the correct one for planet b and provided two transits of planet c at almost exactly the same period predicted by our RV data (21.40 days compared to our predicted value of 21.35 days). Both planets have high-S/N transits (18.6 and 23.0 for planets b and c, respectively). The results of these fits are reported in Table 3.

Table 3. Complete Table of Properties Used in This Analysis

Parameter	Symbol	Value		Units
Stellar Parameters
Mass^a	M_*	0.933 ± 0.054		M_⊙
Radius^a	R_*	0.959 ± 0.042		R_⊙
Age^a	...	${6.8}_{-1.9}^{+1.8}$		Gyr
Stellar Effective Temperature^b	T_eff	5534 ± 100		K
Surface Gravity^b	$\mathrm{log}g$	4.38 ± 0.1		cm s⁻¹
Metallicity^a	[Fe/H]	0.07 ± 0.06		dex
Activity Index^b	logR'_HK	–5.06		...
V-band Magnitude^c	V_mag	8.58		...
J-band Magnitude^c	J_mag	7.30		...
K-band Magnitude^c	K_mag	6.88		...
Distance^d	d	48.8 ± 0.8		pc
Luminosity^d	L	0.798 ± 0.002		L_⊙
Limb Darkening Parameters^e	q₁, q₂	${0.27}_{-0.15}^{+0.28}$ , ${0.26}_{-0.18}^{+0.32}$		...

Transit Parameters ^e		Planet b	Planet c
Period	P	${9.058811}_{-0.000016}^{+0.000017}$	${21.4023}_{-0.00194}^{+0.00189}$	days
Transit Crossing Time	T₀	1494.4462 ± 0.0010	2231.8280 ± 0.0014	TJD^f
Occultation Fraction	R_p/R_*	0.0285 ± 0.0004	0.0277 ± 0.0004	...
Orbital Separation	a/R_*	18.64 ± 1.1	33.06 ± 2.0	...
Inclination	i	88.49 ± 0.0018	${88.241}_{-0.0013}^{+0.0011}$	deg
Transit Duration	T_d	3.36 ± 0.10	4.85 ± 0.33	hr
Transit S/N	S/N	18.6	23.0	...

RV Parameters ^g		Planet b	Planet c
Planet Semiamplitude	K_amp	${3.18}_{-0.52}^{+0.55}$	${2.71}_{-0.49}^{+0.52}$	m s⁻¹
Eccentricity	e	0 (fixed)	0 (fixed)	rad
Periastron Passage	ω	0 (fixed)	0 (fixed)	rad
Instrumental and GP Parameters
Linear HIRES Offset	${\gamma }_{{\mathrm{hires}}_{j}}$	${3.55}_{-0.40}^{+0.46}$		m s⁻¹
Linear APF Offset	γ_apf	${5.49}_{-0.4}^{+0.44}$		m s⁻¹
Linear HARPS-N Offset	γ_HARPS-N	${0.22}_{-3.1}^{+3.0}$		m s⁻¹

Derived Parameters		Planet b	Planet c
Planet Radius	R_p	2.99 ± 0.14	2.90 ± 0.13	R_⊕
Impact Parameter	b	0.49 ± 0.02	${0.30}_{-0.04}^{+0.03}$	...
Planet Mass	M_p	10.8 ± 1.8	11.1 ± 2.4	M_⊕
Planet Density	ρ_p	2.2 ± 0.5	2.5 ± 0.6	g cm⁻³
Insolation Flux	F_p	115.6 ± 4.6	36.7 ± 1.4	S_⊕
Equilibrium Temperature^h	T_eq	911 ± 27	684 ± 21	K
Transmission Spectroscopy Metric	TSM	108.8 ± 35.8	72.5 ± 21.3	...
Emission Spectroscopy Metric	ESM	10.5 ± 1.6	4.8 ± 0.7	...

Notes.

^a isoclassify ^bSpecMatch-Syn^cexofop^dGaia^e juliet ^fBJD 2,457,000^g RadVel ^hAssumes zero albedo and full day–night heat redistribution.

Download table as: ASCII Typeset image

Table 4. Mass Comparison for Different Models

N_planets	M_b (M_⊕)	M_c (M_⊕)	${M}_{{\rm{d}}}\,\sin \,i$ (M_⊕)
2	10.8 ± 1.8	11.1 ± 2.4	n/a
3 (P_d = 58.7 days)	${11.1}_{-1.5}^{+1.6}$	${12.8}_{-2.0}^{+2.1}$	${16.8}_{-2.6}^{+2.7}$
3 (P_d = 101.7 days)	${9.9}_{-1.7}^{+1.8}$	${11.2}_{-2.0}^{+2.2}$	20 ± 4

Notes. These are the mass values for planets b and c for different possible models, as well as the $M\sin i$ for a possible third planet where relevant. Note that the masses for the two transiting planets are consistent within 1σ no matter which model is selected.

Download table as: ASCII Typeset image

5.1.1. LCO Photometry Analysis

As a check on the orbital solution for planet b (particularly on the orbital period), we performed an independent fit of the LCO light curve (Figure 5, Section 4.3.2) using exoplanet (Foreman-Mackey et al. 2021). We fixed T₀ at the value found from the TESS light-curve fit. The fit to the LCO light curve gives 68% credible intervals on the posterior values found for $P={9.058875}_{-0.000034}^{+0.000021}$ , $b={0.60}_{-0.41}^{+0.20}$ , T_d = 3.39 ± 0.14 hr, and ${R}_{{\rm{p}}}/{R}_{* }={0.0325}_{-0.0043}^{+0.0040}$ . These are consistent with the values from the TESS fit (Table 3), though the period values are consistent only within the 2σ errors, and not within the 1σ errors.

5.2. High-resolution Imaging

Within the limits of the AO observations, no stellar companions were detected. 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; M. Lund et. al. 2021, in preparation). 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 TOI-509 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 and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance (Figure 7).

**Figure 7.** Companion sensitivity for the near-infrared AO imaging at Palomar (above) and Keck (below). The black points represent the 5σ limits and are separated in steps of 1 FWHM; the purple area represents the azimuthal dispersion (1σ) of the contrast determinations (see text). The inset image is of the primary target showing no additional close-in companions.
Download figure:
Standard image High-resolution image

5.3. RV Analysis

We used the RadVel ⁴⁹ package (Fulton et al. 2018) to model the RV measurements of HD 63935. RadVel uses the MCMC sampler emcee (Foreman-Mackey et al. 2013) to sample the posterior space of the model's parameters. In these fits, we fixed the period and time of inferior conjunction to the values derived from photometry (previous section). We enforced circular orbits in our fits. Allowing eccentricity to vary produced a fit that was not preferred by the information criterion analysis. Varying e_b or e_c only was somewhat disfavored with ΔAIC_b = 2.28 or ΔAIC_c = 3.34, respectively, while varying both eccentricities was strongly disfavored with ΔAIC_c = 6.02. If the eccentricities were allowed to vary, we would find 1σ upper limits of 0.16 and 0.29 for planets b and c, respectively.

Our preliminary single-planet fits of this system were dominated by a signal at ∼21 days; we identified this signal as corresponding to an additional planet candidate, HD 63935 c, which was confirmed by two transits in TESS Sector 34. The final results of our RV fits are displayed in Figure 8.

**Figure 8.** Our RV points obtained for HD 63935. (a) The complete RV time series, including data from HARPS-N (maroon), APF (green), and HIRES (yellow) and our best-fit model (blue). (b) Our model residuals. Note that the residuals appear to exhibit some structure; we discuss this in Section 5.4. (c) Phase-folded RV curve of HD 63935 b. (d) Phase-folded RV curve of HD 63935 c.
Download figure:
Standard image High-resolution image

5.4. Is There a Third Planet?

The photometry provides clear evidence for two planet candidates, whose presence we confirmed with RV follow-up. The residuals of a two-planet fit, however, show substantial structure. Because the HIRES and APF points in the residuals show similar behavior to each other, instrumental effects are unlikely to be the cause. This suggests that there is something still unaccounted for in our model, which could be a third planet. To test this, we generated a Lomb–Scargle periodogram of the RV residuals from our two-planet fit (Figure 9). There are two primary visible peaks at longer periods, at ∼59 and ∼102 days. Because the peaks in the periodograms of the RV residuals and s-values (activity indicators) do not correspond to each other, we consider it unlikely that this signal is caused by stellar activity. The lack of significant periodicity in the s-values suggests no motivation to adopt a Gaussian process model for our data, consistent with our low value of logR'_HK. We have also performed an independent analysis (Section 3.2) to identify the stellar rotation period, arriving at the conclusion that this period is roughly 30–35 days, and therefore inconsistent with both of the longer-period RV periodogram peaks.

**Figure 9.** Lomb–Scargle periodograms of the RV data for HD 63935, the residuals of the RV data with the signals corresponding to the transiting planets removed, and the Mount Wilson s-values (activity indicators) for the system. False-alarm probabilities of 0.1, 0.01, and 0.001 are shown as horizontal gray lines in each plot. The signals corresponding to a potential long-period companion (purple) do not rise to the required level of statistical significance, but they also do not correspond to peaks in the s-values. The period near 17 days in the middle panel is ruled out as a transiting planet by the TESS photometry but may be an alias of the stellar rotation period.
Download figure:
Standard image High-resolution image

As a method of examining the significance of these peaks, we performed a bootstrap analysis. In this analysis we repeatedly resample our entire RV data set with replacement, calculated the power at the locations in period space where the peak is highest in our true data set, and repeated this 10,000 times. Selecting a p-value of 0.05, we are unable to reject the null hypothesis for either signal. Therefore, in this paper we have adopted the two-planet model, as the evidence for a third planet does not rise to the level of statistical significance we require. Important to note, however, is that this choice does not substantially impact the mass precision of our two confirmed planets, and the resulting masses are consistent to each other to 1σ for both models. Further RV monitoring or photometric follow-up could provide clarity as to the source of this additional signal.

6. Discussion

6.1. Examining Plausible Compositions

In order to better understand this planetary system, we investigated a range of possible compositions based on the planets' bulk density and corresponding positions in mass–radius space. This is of particular interest given the substantial degeneracies in composition that exist for planets in this radius range, primarily between ice/volatile-dominated planets and rock-dominated interiors with substantial H/He envelopes (Lopez & Fortney 2014). To do this, we used the public tool smint ⁵⁰ (Piaulet et al. 2021), developed by Caroline Piaulet, which utilizes the models of Lopez & Fortney (2014) and Zeng et al. (2016) to do an MCMC exploration of the posterior space of planetary interior compositions based on mass, radius, age, and insolation flux. We used 25 MCMC walkers and 10,000 steps for each in this analysis. Our results from smint are that HD 63935 b and c have a 3.6% ± 0.8% and 3.4% ± 0.9% mass in H/He, respectively. Both planets could have cores that are intermediate between ice-dominated and rock-dominated. However, neither planet is sufficiently dense to be a pure "water world"; in other words, both are expected to have substantial (few-percent mass) H/He envelopes. The compositional degeneracies that exist in this part of parameter space are one of the reasons why sub-Neptune-sized planets are so compelling for atmospheric characterization. Study of these planets' atmospheres can potentially reveal more about their composition, which in turn can contribute to our understanding of how these planets form and why no analogs exist in our own solar system.

6.2. Nearly Twins: How Does HD 63935 Fit in the "Peas in a Pod" Structure?

Weiss et al. (2018) identified the phenomenon of Kepler planets in multiplanet systems being more likely to be similar to each other than to be drawn from a random distribution. They referred to this as "peas in a pod." The two confirmed planets in the HD 63935 system appear to be in line with this trend, as both planets have masses and radii that are consistent to each other within 1σ. Higher-precision measurements have the potential to distinguish differences between the two planets. In particular, if planet c is discovered to be of higher density (as nominally appears to be the case, though the difference is at present not statistically significant), this will be interesting, because Weiss & Marcy (2014) noted that the larger and lower-density planet tended to be the exterior in their sample. They suggested that this result could be explained by photoevaporation. In this case, however, the nominally denser planet is the cooler and less irradiated one, implying that a different explanation is required. We calculated the Λ parameter described by Fossati et al. (2017), which estimates whether atmospheric erosion is relevant for a given planet. Regions with Λ < 25 at T_eq = 1000 K and Λ < 35 at T_eq = 500 K (Figure 4 in Fossati et al. 2017) are the relevant regions of atmospheric erosion for G star hosts. The respective Λ values for HD 63935 b and c are 30 and 42, suggesting that neither planet is likely to experience significant atmospheric Jeans escape. In the absence of atmospheric loss, Zeng et al. (2019) proposed that denser outer planets could be explained via impacts of ice-rich planetesimals (see also Marcus et al. 2009).

Also potentially of interest are the differing equilibrium temperatures of these planets (∼911 K and ∼684 K for b and c, respectively). Given their otherwise similar properties, they could serve as an experimental testing ground for the role of insolation in planetary composition and nature. As we describe in more detail in the following section, planets cooler than 1000 K are likely to have decreasing spectral feature amplitude in transmission spectra as a result of increased haze formation (Gao et al. 2020). Observation of such a phenomenon in the HD 63935 system could be used as further evidence for this trend and could be of particular significance because the planets orbit the same star (eliminating a possible confounding variable).

6.3. Assessing Atmospheric Observability

As described above in Section 2, we identified HD 63935 b as the best target for atmospheric characterization follow-up in its region of parameter space (between 2.6 and 4 Earth radii, between 10 and 100 times the Earth insolation flux, and stellar effective temperatures between 5200 and 6500 K). The TSM value for HD 63935 b, incorporating our measured mass, is 108.8 ± 35.8. In addition to its uniqueness within our algorithm's defined parameter space bins, planet b is of interest as a tool to probe the radius cliff, which refers to the apparent drop-off in planet occurrence rate above ∼2.5 R_⊕. In Figure 10, we identify planets with TSM values equal to or greater than that of HD 63935 b, and find that only one (HD 191939 b; Lubin et al. 2021) falls on the radius cliff. The physical causes of the radius cliff are hypothesized to be related to atmospheric sequestration (Kite et al. 2019), signs of which may be visible in an atmosphere, further enhancing the target's desirability for characterization. Note also that of the planets in that figure, four (HD 219134 b, HD 219134 c, 55 Cnc e, and π Mensae c) have host star magnitudes that saturate JWST, meaning they are not suitable targets for transmission spectroscopy with JWST.

**Figure 10.** The position of HD 63935 b in mass–radius space. Top: The planet sample from Zeng et al. (2019), colored by equilibrium temperature, as well as a simple subset of composition curves from Zeng et al. (2016) to contextualize the image. Bottom: A density estimation of the same planet sample, with the host star spectral type indicated by a marker symbol (F, G, K, or M). The bimodal distribution of planets around the radius gap at ∼1.8 R_⊕ is clearly visible. HD 63935 b and c are emphasized and colored red, with planet b the upper of the two red points. Planets with TSM values higher than that of HD 63935 b are colored cyan (from top to bottom, GJ 1214 b, HD 97658 b, π Mensae c, GJ 9827 d, 55 Cancri e, HD 219134 b, and HD 219134 c), and planets that have been subjected to atmospheric characterization in the past are circled in black. Those that have published atmospheric spectra but lower TSM values (from top to bottom, HD 3167c, K2-18 b, and HD 3167 b) are colored white. Planets meeting neither criterion are faded gray. The dark red lines correspond to a subset of plausible composition curves for HD 63935 b, from Lopez & Fortney (2014) and Zeng et al. (2016). This figure emphasizes the planet's uniqueness as a quality atmospheric target in parameter space.
Download figure:
Standard image High-resolution image

To validate the observability of HD 63935 b, we simulated planetary transmission spectra using CHIMERA (Line et al. 2012, 2013, 2014), then simulated transit observations of the planets using PandEXO (Batalha et al. 2017). A sample of our simulations (cloud-free cases only) are shown in Figure 11.

**Figure 11.** This figure shows the simulated near-infrared transmission spectrum of HD 63935 b based on a single transit, with the NIRISS SOSS and NIRSpec G395H instruments (note that this means that the full data set shown in this plot would require two transits to obtain, one for each instrument and corresponding wavelength range). Different colors indicate different metallicities relative to the solar value. The true model spectra are shown as solid lines, while the points are the simulated observations.
Download figure:
Standard image High-resolution image

The results for planet b support our selection of the planet as one with an exceptionally high potential for atmospheric characterization. We simulated a range of metallicities and cloud opacities, then calculated the S/N of the water feature at 1.5 μm. We did this by simulating the spectrum 1000 times, determining the equivalent width of the feature in each, and finding the standard deviation of the resulting distribution. Our results are displayed in Figure 12. We emphasize that this is a conservative estimate of the quality of the atmospheric spectra, but that even with middlingly optimistic assumptions about metallicity and cloudiness, an observation of a single transit with JWST can produce useful spectra of this planet, and more observations would of course produce correspondingly more precise spectra.

**Figure 12.** The S/N of the 1.5 μm water feature equivalent width of HD 63935 b, calculated after a single JWST transit observation for a variety of metallicities and cloud opacities (log(kcld)). log(kcld) is a model parameter that encapsulates cloud opacity in a single graying parameter; the range used here spans the gamut from cloud-free atmospheres in light green to opaque atmospheres in black. Optimistic but reasonable atmospheric clearness assumptions are between −30 and −29 (so S/N ∼ 10 for the 100× solar metallicity case). Obviously, an opaque atmosphere would produce no visible water features, though we describe in the text why not to expect such a scenario. Note as well that this is a very conservative estimate of the quality of atmospheric observations. There are other water features available for obtaining a water detection, the equivalent width calculation used here ignores binning as well as more advanced retrieval techniques, and the data here come from only a single JWST transit. The S/N should scale roughly with $\sqrt{{N}_{\mathrm{obs}}}$ if more transits are added.
Download figure:
Standard image High-resolution image

**Figure 12.** The S/N of the 1.5 μm water feature equivalent width of HD 63935 b, calculated after a single JWST transit observation for a variety of metallicities and cloud opacities (log(kcld)). log(kcld) is a model parameter that encapsulates cloud opacity in a single graying parameter; the range used here spans the gamut from cloud-free atmospheres in light green to opaque atmospheres in black. Optimistic but reasonable atmospheric clearness assumptions are between −30 and −29 (so S/N ∼ 10 for the 100× solar metallicity case). Obviously, an opaque atmosphere would produce no visible water features, though we describe in the text why not to expect such a scenario. Note as well that this is a very conservative estimate of the quality of atmospheric observations. There are other water features available for obtaining a water detection, the equivalent width calculation used here ignores binning as well as more advanced retrieval techniques, and the data here come from only a single JWST transit. The S/N should scale roughly with $\sqrt{{N}_{\mathrm{obs}}}$ if more transits are added.
Download figure:
Standard image High-resolution image

Additionally, there is reason to expect that the atmosphere of HD 63935 b has high-amplitude atmospheric features in its transmission spectrum. Although cloud formation in exoplanetary atmospheres is not fully understood, extensive modeling and analysis of existing transmission spectra have been done to attempt to understand which atmospheres are dominated by clouds or hazes. Most recently, Gao et al. (2020) used an aerosol microphysics model to predict the dominant opacity sources in giant-exoplanet transmission spectra. Their results suggest that the height of spectral features (specifically of the 1.4 μm water band) increases with increasing temperature until ∼950 K, with opacity in this temperature regime dominated by high-altitude photochemical hazes derived from hydrocarbons. Above ∼950 K, silicate condensate clouds become the dominant opacity source. These clouds rise higher into the atmosphere as temperature increases, resulting in reduced spectral feature amplitude with increasing temperature from 950 K until ∼1800 K. These predictions are highly relevant to this work, as HD 63935 b sits quite close to the critical equilibrium temperature of 950 K, which we expect to be the local maximum of feature amplitude.

This prediction is also in line with earlier work by Crossfield & Kreidberg (2017), who analyzed the published transmission spectra of exclusively Neptune-sized exoplanets and found a positive correlation between feature height and equilibrium temperature in their domain of 500–1100 K. Taken together, these results suggest that HD 63935 b is in a region of parameter space that is maximally likely to produce observable spectral features. Combined with its other desirable qualities as an atmospheric target (bright host star, large TSM value), it is clear that HD 63935 b has exceptional promise as a target for atmospheric characterization.

Finally, we note that HD 63935 is quite Sun-like (M_* ∼ 0.933 M_⊙, [Fe/H] = 0.07 ± 0.06), making characterization of its planetary system of additional interest to comparative planetology with our own solar system. We also emphasize that only one sub-Neptune planet around a similar star has a published atmospheric characterization study (HD 3167 c; Mikal-Evans et al. 2020). Combined with other such planets with G-type host stars being discovered by TESS, it could form part of a robust population study.

7. Conclusions

We have described the discovery of two planets around the bright G star HD 63935.

1.
HD 63935 b and c have periods of ${9.0600}_{-0.00070}^{+0.00064}$ and 21.40 ± 0.0019 days, radii of 2.99 ± 0.14 and 2.90 ± 0.13 R_⊕, and masses of 10.8 ± 1.8 and 11.1 ± 2.4 M_⊕.
2.
Both confirmed planets have a radius greater than the median sub-Neptune radius and fall on the radius cliff.
3.
Planet b is an outstanding target for transmission spectroscopy, being the best target in its parameter space niche and the second-best among targets on the radius cliff, while planet c is also amenable to atmospheric characterization. This quality for follow-up is exceptionally enticing given the compositional degeneracies that exist in this region of parameter space, which could be broken with high-quality atmospheric characterization.
4.
Atmospheric characterization of this system could provide valuable input to theories of planetary interiors, formation, and evolution, especially given the fact that the two planets similar in all observed properties except insolation flux are as close to a variable-controlled experimental setup as exoplanet astronomy typically comes.

HD 63935 b and c attest to the bright future of exoplanet astronomy and we expect this system to be an excellent test case for studying exoplanetary atmospheres in the coming years.

We thank the anonymous referee for helpful feedback that improved the quality of this work. We thank the time assignment committees of the University of California, the California Institute of Technology, NASA, and the University of Hawai'i for supporting TKS with observing time at Keck Observatory and on APF. We thank NASA for funding associated with our Key Strategic Mission Support project. We gratefully acknowledge the efforts and dedication of the Keck Observatory staff for support of HIRES and remote observing. We recognize and acknowledge the cultural role and reverence that the summit of Maunakea has within the indigenous Hawaiian community. We are deeply grateful to have had the opportunity to conduct observations from this mountain. We thank Ken and Gloria Levy, who supported the construction of the Levy Spectrometer on APF. We thank the University of California and Google for supporting Lick Observatory and the UCO staff for their dedicated work scheduling and operating the telescopes of Lick Observatory.

This paper is based on data collected by the TESS mission. Funding for the TESS mission is provided by NASA's Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS SPOC. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Resources supporting this work were provided by the NASA High-end Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center for the production of the SPOC data products. This paper includes data collected by the TESS mission that are publicly available from MAST. This research has made use of the NASA Exoplanet Archive and Exoplanet Follow-up Observation Program website, which are operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Based on observations made with the Italian TNG operated on the island of La Palma by Fundación Galileo Galilei of Istituto Nazionale di Astrofisica at the Spanish Observatorio del Roque de los Muchachos of Instituto de Astrofisica de Canarias under programs CAT19A_162 and CAT19A_96. This work has made use of data from the European Space Agency mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work is partly supported by JSPS KAKENHI grant Nos. JP17H04574 and JP18H05439, JST PRESTO grant No. JPMJPR1775, and a Grant-in-aid for JSPS Fellows, No. JP20J21872. This work makes use of observations from the LCOGT Network.

D.D. acknowledges support from the TESS Guest Investigator Program grant 80NSSC19K1727 and the NASA Exoplanet Research Program grant 18-2XRP18_2-0136. J.M.A.M. is supported by the National Science Foundation (NSF) Graduate Research Fellowship Program under grant No. DGE-1842400. J.M.A.M. acknowledges the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining grant No. 1829740, the Brinson Foundation, and the Moore Foundation; his participation in the program has benefited this work. M.R.K. is supported by an NSF Graduate Research Fellowship, grant No. DGE 1339067. P.A.D. acknowledges support from an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1903811.

Facilities: APF (Levy), HARPS-N (TNG), HIRES (Keck I), LCOGT, NIRC2 (Keck II), PHARO (Palomar), TESS.

Software: AstroImageJ (Collins et al. 2017), Astropy (Robitaille et al. 2013), batman (Kreidberg 2015), emcee (Foreman-Mackey et al. 2013), exoplanet and its dependencies (Foreman-Mackey et al. 2021; Agol et al. 2020; Kumar et al. 2019; Luger et al. 2019; Salvatier et al. 2016; Theano Development Team 2016), isoclassify (Huber et al. 2017), juliet (Espinoza et al. 2019), Jupyter (Kluyver et al. 2016), kiauhoku (Claytor et al. 2020), matplotlib (Hunter 2007), numpy (Harris et al. 2020), pandas (McKinney 2010), Radvel (Fulton et al. 2018), smint (Piaulet et al. 2021), SpecMatch-Syn (Petigura et al. 2017), Transit Least Squares (Hippke & Heller 2019).

TESS-Keck Survey. V. Twin Sub-Neptunes Transiting the Nearby G Star HD 63935

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. TKS: Atmospheric Target Selection