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

As an independent determination of the stellar parameters, we also performed an analysis of the broadband spectral energy distribution (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); 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 2MASS, the W1–W4 magnitudes from WISE, the $GG_{\rm BP}G_{\rm RP}$ magnitudes from Gaia, and the FUV and NUV magnitudes from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.15–22 $\mu$m (see Figure 2).

We performed a fit using Kurucz stellar atmosphere models, with the effective temperature ($T_{\rm eff}$), metallicity ([Fe/H]), and surface gravity ($\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 $\chi^{2}$ of 1.1 and best-fit $A_{V}=0.02\pm 0.02$. Integrating the (unreddened) model SED gives the bolometric flux at Earth, $F_{\rm bol}=1.060\pm 0.012\times 10^{-8}$ erg s^-1 cm^-23. Taking the $F_{\rm bol}$ and $T_{\rm eff}$ together with the Gaia EDR3 parallax, gives the stellar radius, $R_{\star}=0.967\pm 0.035R_{\odot}$. In addition, we can use the $R_{\star}$ together with the spectroscopic $\log g$ to obtain an empirical mass estimate of $M_{\star}=0.82\pm 0.20M_{\odot}$, which is consistent with that obtained via empirical relations of Torres et al. (2010), $M_{\star}=1.02\pm 0.06M_{\odot}$. These parameters are also consistent with those we derive from isoclassify. Finally, the $R_{\star}$ and $M_{\star}$ together yield a mean stellar density of $\rho_{\star}=1.59\pm 0.20$ g cm^-3.

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 radial velocity signal at 21 days. Although the 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 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_{\mathrm{rot}}=32.4\pm 6.5$ days, while we predicted $P_{\mathrm{rot}}=31.1\pm 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 non-detection of a rotation period in the TESS data.

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 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 33208 data points with integration times of two minutes each. The Science Processing Operations Center (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 (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 lightcurve (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)³³3https://archive.stsci.edu/tess/. The full light curve is plotted in Figure 3 and the phase-folded transits in Figure 4.

As part of our standard process for validating transiting exoplanets to assess the 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 Apr 18 UT in a standard 5-point quincunx dither pattern with steps of 5″ in the narrow-band $Br-\gamma$ filter $(\lambda_{o}=2.1686;\Delta\lambda=0.0326~\mu$m). Each dither position was observed three times, offset in position from each other by 0.5″ for a total of 15 frames; with an integration time of 1.4 seconds per frame, the total on-source time was 21 seconds on target. PHARO has a pixel scale of $0.025\arcsec$ per pixel for a total field of view of $\sim 25\arcsec$. 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-Mar-25 UT in the standard 3-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\arcsec$ and was repeated twice, with each dither offset from the previous dither by $0.5\arcsec$. NIRC2 was used in the narrow-angle mode with a full field of view of $\sim 10\arcsec$ and a pixel scale of approximately $0.0099442\arcsec$ per pixel. The Keck observations were made in the narrow-band $Br-\gamma$ filter $(\lambda_{o}=2.1686;\Delta\lambda=0.0326~\mu$m) with an integration time of 0.5 second for a total of 4.5 seconds on target. Observations were taken in narrow camera mode with a 1024” x 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 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 15 dithered science frames; each science image was then sky-subtracted and flat-fielded. The reduced science frames were combined into a single combined image using a intra-pixel interpolation that conserves flux, shifts the individual dithered frames by the appropriate fractional pixels, and median-coadds the frames (Figure 6). The final resolution of the combined dithers was determined from the full-width half-maximum of the point spread function; 0.094″ and 0.050″ for the Palomar and Keck observations respectively.

The TESS pixel scale is $\sim 21\arcsec$ pixel^-1, and photometric apertures typically extend out to roughly 1 arcminute, 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 (TFOP)⁴⁴4https://tess.mit.edu/followup 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.

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 modelling used the functions based on batman (Kreidberg, 2015) for transit fitting and 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 two quadratic limb darkening parameters (following Kipping (2013), the out-of-transit flux, and a jitter term from TESS. We keep all other parameters fixed, including the mean stellar density, for which we fix the value to that derived in Section 3. We impose 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 constrain 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 ($\sim 18$ days, the one initially reported by SPOC, and $\sim$9 days, which we ultimately selected based on RV measurements and was later confirmed by TESS S34) 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 radial velocity 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 2 transits of planet c at almost exactly the period predicted by our radial velocity data (21.40 days compared to our predicted value of 21.35 days). Both planets have high SNR transits (18.6 and 23.0 for planets b and c, respectively). The results are of these fits are reported in Table 3.

To 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^{\circ}$ at separations of integer multiples of the central source’s FWHM (Furlan et al., 2017; Lund, Submitted). The brightness of each injected source was scaled until standard aperture photometry detected it with $5\sigma$ significance. The resulting brightness of the injected sources relative to TOI 509 set the contrast limits at that injection location. The final $5\sigma$ 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).

We used the RadVel⁷⁷7https://radvel.readthedocs.io/en/latest/ package (Fulton et al., 2018) to model the radial velocity measurements of HD 63935. RadVel uses the Markov chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al., 2013a) to sample the posterior space of the model’s parameters. In these fits, we fix the period and time of inferior conjunction to the values derived from photometry (previous section). We enforce circular orbits in our fits. Allowing eccentricity to vary produced a fit which was not preferred by the information criterion analysis. Varying $e_{b}$ or $e_{c}$ only were “somewhat disfavored” with $\Delta$AICc$=2.28$ and $\Delta$AICc$=3.34$, respectively, while varying both eccentricities was “strongly disfavored” with $\Delta$AICc$=6.02$). If the eccentricities are allowed to vary, we find 1$\sigma$ 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 $\sim$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 radial velocity fits are displayed in Figure 8.

The photometry provides clear evidence for two planet candidates, whose presence we confirm with radial velocity followup. 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 suggested that there was something still unaccounted for in our model, which could be a third planet. To test this, we generated a Lomb-Scargle periodogram of the radial velocity residuals from our two-planet fit (Figure 9). There are two primary visible peaks at longer periods, at $\sim$59 and $\sim$102 days. Because the peaks in the periodograms of the radial velocity 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 $\log$R’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.

As a method of examining the significance of these peaks, we perform a bootstrap analysis. In this analysis we repeatedly resample our entire RV dataset with replacement, calculate the power at the locations in period space where the peak is highest in our true dataset, and repeat this 10000 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-$\sigma$ for both models. Further radial velocity monitoring or photometric followup could provide clarity as to the source of this additional signal.

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⁸⁸8https://github.com/cpiaulet/smint (Piaulet et al., 2021), developed by Caroline Piaulet, which utilizes the models of Lopez & Fortney (2014) and Zeng et al. (2016) to do a Markov Chain Monte Carlo exploration of the posterior space of planetary interior compositions based on their mass, radius, age, and insolation flux. We use 25 MCMC walkers and 10000 steps for each in this analysis. Our results from smint are that HD 63935 b and c have $3.6\pm 0.8$% and $3.4\pm 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 that the 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.

Weiss et al. (2018) identified the phenomenon of Kepler planets in multi-planet systems being more likely to be similar to each other than drawn from a random distribution. They refer 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 one sigma. Higher precision measurements have the potential to distinguish differences between the two planets. In particular, if planet c is discovered to be higher density (as nominally appears to be the case, though the difference is not presently statistically significant), this would be interesting, because Weiss & Marcy (2014) note that the larger and lower density planet tends to be the exterior in their sample. They suggest 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 calculate the $\Lambda$ parameter described by Fossati et al. (2017), which estimates whether atmospheric erosion is relevant for a given planet. $\Lambda<25$ at T${}_{\textrm{eq}}=1000K$ and $\Lambda<35$ at T${}_{\textrm{eq}}=500K$ (Fig 4 in Fossati et al. (2017)) are the relevant regions of atmospheric erosion for G star hosts. The respective $\Lambda$ 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) propose that denser outer planets could be explained via impacts of ice-rich planetesimals (see also Marcus et al., 2009).

Also potentially of interest is the differing equilibrium temperatures of these planets ($\sim$911K and $\sim$684K for b and c, respectively). Given their otherwise similar properties, they could serve as an experimental testing ground for the role of insolation on planetary composition and nature. As we describe in more detail in the following section, planets cooler than 1000K 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 of particular significance because the planets orbit the same star (eliminating a possible confounding variable).

As described above in Section 2, we identified HD 63935 b as the best target for atmospheric characterization followup in its region of parameter space (between 2.6 and 4 Earth radii, between 10 and 100 times Earth insolation flux, and stellar effective temperatures between 5200 and 6500K). The TSM value for HD 63935 b, incorporating our measured mass, is 108.8$\pm$ 35.8. In addition to its uniqueness within our algorithm’s defined parameter space bins, planet b is also of interest as a tool to probe the radius cliff, which refers to the apparent drop-off in planet occurrence rate above $\sim 2.5R_{\oplus}$. 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. (submitted)) 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 $\pi$ Mensae c) have host star magnitudes that saturate JWST , meaning they are not suitable targets for transmission spectroscopy with JWST .

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.

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 SNR of the water feature at 1.5$\mu$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 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.

Additionally, there is reason to expect that the atmosphere of HD 63935 b will have high-amplitude atmospheric features in its transmission spectrum. Although cloud formation in exoplanetary atmospheres is not fully understood, extensive modelling and analysis of existing transmission spectra have been done to attempt to understand which atmospheres will be 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 the 1.4 $\mu m$ water band) increases with increasing temperature until $\sim$950 K, with opacity in this temperature regime dominated by high-altitude photochemical hazes derived from hydrocarbons. Above $\sim$950K, 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 until $\sim$1800K. These predictions are highly relevant to this work, as HD 63935 b sits quite close to the critical equilibrium temperature of 950K which we expect to be a local maximum of feature amplitude.

This prediction is also in line with earlier work by Crossfield & Kreidberg (2017), which 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-1100K. 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${}_{*}~\sim 0.933\textrm{M}_{\odot}$, [Fe/H] = 0.07 $\pm$ 0.06), making characterization of its planetary system of additional interest for comparative planetology with our own solar system. We also emphasize that only one sub-Neptune planet around such 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.