TESS Hunt for Young and Maturing Exoplanets (THYME). III. A Two-planet System in the 400 Myr Ursa Major Group

The TESS mission (Ricker et al. 2014) observed TIC 130181866 (TOI 1726, HD 63433, HIP 38228) between 2019 December 24 and 2020 January 21 (Sector 20) using Camera 1. The target was proposed by three guest investigator programs (G022032, T. Metcalfe; G022038, R. Roettenbacher; G022203, J. Ge) and hence has has 2-minute cadence data. The abstracts for these GI programs suggest that HD 63433 was targeted because it is bright (V ≃ 7 mag) and/or active.

For our analysis, we used the Presearch Data Conditioning simple aperture photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2014) TESS light curve produced by the Science Process Operations center (SPOC; Jenkins et al. 2016) and available through the Mikulski Archive for Space Telescopes (MAST).⁴² We only included data points with DQUALITY = 0, i.e., those with no flags from the SPOC pipeline. No obvious flares were present in the data, so we did no further data processing.

We obtained time series photometry during three predicted transits of HD 63433 b using ground-based facilities in order to rule out nearby eclipsing binaries (NEBs) that could be the source of the transit signal. We observed an egress on 2020 February 22 UT and a full transit on 2020 February 29, both in r with the 0.6 m World Wide Variable Star Hunters (WWVSH) telescope in Abu Dhabi, United Arab Emirates. Both transits were observed under clear (near-photometric) conditions. The telescope is equipped with a Finger Lakes Instrumentation FLI 16803 camera, giving a pixel scale of 047 pixel⁻¹. We obtained 118 and 350 exposures with an exposure length of 180 and 90 s for the two observations, respectively. We also observed a full transit of the 2020 February 29 event with one of the 1 m telescopes at the Las Cumbres Observatory (LCOGT; Brown et al. 2013) node at the South African Astronomical Observatory, South Africa, under clear conditions. We observed in the z_s band using a Sinistro camera, giving a pixel scale of 0389 pixel⁻¹. We obtained 161 60 s exposures, which were reduced using the BANZAI pipeline (McCully et al. 2018). In all cases, we deliberately saturated the target star in order to search for faint NEBs.

We performed aperture photometry on all three data sets using the AstroImageJ package (AIJ; Collins et al. 2017). TIC 130181879 (TESS magnitude 13.1) and TIC 130181877 (TESS magnitude 14.6) are the only two stars within 25 of HD 63433 that are bright enough to cause the TESS detection, so we extracted light curves for both. We used the aperture size and set of comparison stars that yielded the best precision for each observation and star of interest. For both WWVSH observations, we used a 7 pixel (33) radius circular aperture to extract the source and an annulus with a 12 pixel (56) inner radius and a 17 pixel (8'') outer radius for the sky. For the LCO photometry, we used a 6 pixel (23) radius circular aperture for the source and an annulus with a 16 pixel (62) inner radius and a 23 pixel (9'') outer radius for the sky. For all observations, we centered the apertures on the source, weighted all pixels within the aperture equally, and used the same aperture setup for TIC 130181879 and TIC 130181877. The extracted light curves (including for additional nearby faint sources), field overlays, and further information on this follow-up can be found on ExoFOP-TESS.⁴³

To reproduce the observed transit depths, eclipses around either nearby star would need to be large (≃20%), but no binary was detected down to ≲1%. Thus, the observations ruled out any NEB scenario consistent with the observed transit.

We utilized new and archival high-resolution spectra and radial velocity measurements of HD 63433 in our analysis. We list all of the radial velocity data in Tables 1 and 2.

2.3.1. LCOGT/NRES

2.3.2. Tillinghast/TRES

2.3.3. Goodman/SOAR

2.3.4. HARPS-N/TNG

Between 1997 March and 2016 April, HD 63433 was observed 15 times from the 1.93 m telescope at the Haute-Provence Observatory located in France. The first three were taken with the ELODIE high-resolution spectrograph (Baranne et al. 1996), and the next 12 were taken by ELODIE's replacement, SOPHIE (Perruchot et al. 2008). We retrieved the spectra and barycentric radial velocities given on the SOPHIE/ELODIE archives (Moultaka et al. 2004).⁴⁶ To correct for differences in the zero-point between ELODIE and SOPHIE, we apply an offset of 87 ± 23 m s⁻¹ to ELODIE velocities as described in Boisse et al. (2012). SOPHIE velocities were all taken after the upgrade to SOPHIE+ (Bouchy et al. 2013) and have formal uncertainties of 1–3 m s⁻¹, not including stellar jitter or long-term drift in the instrument (≃5 m s⁻¹; Courcol et al. 2015).

As part of the Lick planet search program, HD 63433 was observed 11 times between 1998 January and 2011 December using the Hamilton Spectrograph and iodine cell (Vogt 1987) at Lick Observatory in California, USA. We utilize the velocities and errors reported in Fischer et al. (2014). Velocity errors from the Lick planet search include instrument stability but do not account for stellar jitter. These are relative velocities (the star compared to itself) and hence cannot be directly compared to other measurements without modeling an offset (Díaz et al. 2016).

The Ursa Major Group (UMaG) has long been proposed as a kinematically similar grouping of stars (e.g., Proctor 1869; Rasmuson 1921; Eggen 1965) centered on several of the stars comprising the Ursa Major constellation. While UMaG has a clear core of members that are homogeneous in kinematics and color–magnitude diagram position (Soderblom & Mayor 1993; King et al. 2003), many associations with large spatial distributions have turned out to be larger star formation events with multiple ages (e.g., Sco-Cen and Taurus-Auriga; Rizzuto et al. 2011; Kraus et al. 2017). Further, the spatial spread of UMaG members outside the core leads to a large number of interloping stars with similar Galactic orbits but different ages. Tabernero et al. (2017) found that $\simeq 2/3$ of UMaG members have similar chemical compositions, suggesting either multiple stellar populations or a large fraction of contaminants in the membership list. Thus, to be useful for age-dating HD 63433, we need to establish its association to the core members of UMaG using both kinematics and independent metrics (e.g., rotation and abundances).

UMaG has recent age estimates ranging from 390 to 530 Myr (e.g., Brandt & Huang 2015; David & Hillenbrand 2015). Direct radius measurements from long-baseline interferometry for the A stars in UMaG point toward a common age of τ = 414 ± 23 Myr (Jones et al. 2015). We adopt this measurement as the cluster age for analyses in this paper.

HD 63433 was first identified as a possible member of UMaG by Gaidos (1998) based on its kinematics and X-ray luminosity and has since been included as a candidate or likely member in multiple analyses (e.g., King et al. 2003; Fuhrmann 2008 ;Vereshchagin et al. 2018). The spatial and kinematic definition of UMaG was most recently updated by Gagné et al. (2018), who identified central values of $(X,Y,Z)=(-7.5,+9.9,+21.9)$ pc and $(U,V,W)=(+14.8,+1.8,-10.2)$ km s⁻¹, along with full covariance matrices for these parameters. These are marginally consistent with the central values from Mamajek et al. (2010) of $(U,V,W)=(+15.0,2.8,-8.1)$ ± (0.4, 0.7, 1.0) km s⁻¹. The Gaia DR2 proper motion, parallax, and our radial velocity for HD 63433 give $(U,V,W)=(+13.66,+2.42,-7.75)$ km s⁻¹. Following the method from Gagné et al. (2018),⁴⁷ we calculate a membership probability of P_mem = 99.98% using the Gagné et al. (2018) space velocity for UMaG and P_mem = 95% using the space velocity from Mamajek et al. (2010).

HD 63433's Galactic position of (X, Y, Z) = (−19.89, −4.70, 9.16) pc is 23 pc from the core of UMaG at (X, Y, Z) = (−7.5, 9.9, 21.9) pc. While the core members of UMaG are within ≃4 pc of the core, more than half of known members are >20 pc away (Madsen et al. 2002). The large-scale velocity dispersion is ≃1 km s⁻¹ (≃1 pc Myr⁻¹), and the age is τ ∼ 400 Myr; members could easily be spread over >100 pc in Y (where orbits can freely diverge) and a still substantial distance even in X and Z (where epicyclic motion prevents the spatial distribution from broadening to the same degree). Mamajek et al. (2010) and Schlieder et al. (2016) argue that the measured dispersion of ≃1 km s⁻¹ is mostly an artifact of including spectroscopic binaries in the sample and that the true dispersion is smaller. However, given the age of the cluster, HD 63433 only needs to have a velocity difference of ≃0.06 km s⁻¹ to explain a 23 pc separation. It is more likely that the velocity difference is larger and HD 63433 was evaporated from the cluster core in the past 100 Myr.

The photospheric lithium abundance provides an age and membership diagnostic that is independent of the 6D position–velocity phase space of HD 63433. Li is destroyed at temperatures common to the cores of stars (∼2.5 × 10⁶ K), which slowly depletes surface Li at a rate that depends on the core–surface transport efficiency (e.g., convection). While there is significant scatter within a single age group (e.g., Somers & Pinsonneault 2015), there is still a shift in the average A(Li)–${T}_{\mathrm{eff}}$ sequence with age.

We compared the A(Li) abundance of HD 63433 from Ramírez et al. (2012) to A(Li) for members of Pleiades (125 Myr; Dahm 2015) and Hyades (700 Myr; Martín et al. 2018). A(Li) measurements for the Pleiades were taken from Bouvier et al. (2018), and values for the Hyades were taken from Boesgaard et al. (2016). We also considered UMaG members that were confirmed using kinematics and chromospheric activity by King et al. (2003). We retrieved A(Li) from King & Schuler (2005), Ramírez et al. (2012), Aguilera-Gómez et al. (2018), and the Hypatia catalog (Hinkel et al. 2014). HD 63433 has A(Li) between that of Hyades and Pleiades stars of similar ${T}_{\mathrm{eff}}$ and is consistent with the core members of UMaG, furthering the case for membership (Figure 1).

Zoom In Zoom Out Reset image size

Stellar rotation provides an additional check on the age and membership of HD 63433. Once on the main sequence, young stars lose angular momentum with time, decreasing their rotation periods. After 100–600 Myr, Sun-like stars eventually converge to a sequence (van Saders et al. 2016; Douglas et al. 2016). Angus et al. (2015) used this information to provide a calibration between B − V, P_rot, and age, which predicts a rotation period of 6.9 ± 0.4 days for HD 63433 if it is a member of UMaG.

We estimated the rotation period of HD 63433 and other likely UMaG members from their TESS or K2 light curves using a combination of the Lomb–Scargle periodogram following Horne & Baliunas (1986) and the autocorrelation function as described in McQuillan et al. (2013). For this, we used the simple aperture photometry (SAP) light curves, as PDCSAP curves tend to have long-term signals removed or suppressed (e.g., Smith et al. 2012; Van Cleve et al. 2016). For HD 63433, this yielded a period estimate of 6.45 days (Figure 2). The Lomb–Scargle power is relatively broad, although the power is high and bootstrap resampling of the light curve suggest errors on this period of ≲3%. Our derived period is also consistent with the literature estimate of 6.46 ± 0.01 days from Gaidos et al. (2000). Rotation periods for other likely members are listed in Table 4. We note that because of the narrow window provided by TESS photometry for many UMaG members (<30 days), our estimates are subject to aliasing (i.e., true periods may be double or half the assigned value) separate from the smaller formal errors (≃3%).

Zoom In Zoom Out Reset image size

As with A(Li), we compare our UMaG rotation periods to those from the older Praesepe cluster (from Douglas et al. 2019) and younger Pleiades cluster (from Rebull et al. 2016). We show the results in Figure 3. There is significant scatter in each sequence (e.g., from binarity and nonmember interlopers). The UMaG sequence lands just below the Praesepe sequence in period space, and most members are above (longer periods than) the Pleiades sequence. As expected, HD 63433 follows the overall trend for UMaG. The rotation sample considered here is small, and a more complete accounting of UMaG membership is needed to explore the overlap between Praesepe and UMaG rotation sequences.

Zoom In Zoom Out Reset image size

Altogether, the available evidence confirms the age and membership of HD 63433. For all analyses in the rest of the paper, we adopt the cluster age (414 ± 23 Myr) as the age of HD 63433.

We fit HD 63433's spectral energy distribution using available photometry, our Goodman optical spectrum (Section 2.3.3), and spectral templates of nearby stars (e.g., Rayner et al. 2009; Falcón-Barroso et al. 2011). To this end, we followed the basic methodology of Mann et al. (2015). The procedure gives precise (1%–5%) estimates of ${F}_{\mathrm{bol}}$ from the integral of the absolutely calibrated spectrum, L_* from ${F}_{\mathrm{bol}}$ and the Gaia distance, and ${T}_{\mathrm{eff}}$ from comparing the calibrated spectrum to atmospheric models. This method reproduces angular diameter measurements from long-baseline optical interferometry (e.g., von Braun et al. 2012). As a check, the same procedure also provides an estimate of R_* from the infrared-flux method (Blackwell & Shallis 1977), i.e., the ratio of the absolutely calibrated spectrum to the model spectrum is ${R}_{* }^{2}/{D}^{2}$ (also see Equation (1) of Cushing et al. 2008).

We first combined our template and observed spectra with Phoenix BT-SETTL models (Allard et al. 2011) to cover wavelength gaps and regions of high telluric contamination and assumed a Rayleigh–Jeans law redward of where the models end (25–30 μm). To absolutely calibrate the combined spectra, we used literature optical and near-IR photometry from the Two-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006), the Wide-field Infrared Survey Explorer (WISE; Cutri et al. 2014), Gaia data release 2 (DR2; Evans et al. 2018; Lindegren et al. 2018), AAVSO All-Sky Photometric Survey (Henden et al. 2016), Tycho-2 (Høg et al. 2000), Hipparcos (van Leeuwen et al. 1997), and the General Catalog of Photometric Data (Mermilliod et al. 1997). To account for variability of the source, we assumed that all photometry had an addition error of 0.02 mag (for optical) or 0.01 mag (for near-infrared), and filter zero-points are assumed to have errors of 0.015 mag unless a value is given in the source. We then compared the literature photometry to synthetic magnitudes derived from combined spectrum using the relevant filter profiles and zero-points (e.g., Mann & von Braun 2015; Maíz Apellániz & Weiler 2018). We assumed no reddening, as HD 63433 lands within the Local Bubble (Sfeir et al. 1999). In addition to the overall scale of the spectrum, there are four free parameters of the fit that account for imperfect (relative) flux calibration of the spectra and both the model and template spectra used.

We show the best-fit spectrum and photometry in Figure 4. There is no significant near-IR excess seen out to W4, consistent with most stars at this age (Cieza et al. 2008). Our joint fitting procedure yielded ${T}_{\mathrm{eff}}$ = 5640 ± 71 K, ${F}_{\mathrm{bol}}$ = $(4.823\pm 0.12)\,\times {10}^{-8}$ erg cm⁻² s⁻¹, ${L}_{* }=0.753\pm 0.026\,{L}_{\odot }$, and R_* = 0.912 ± 0.034 R_⊙. Our derived ${T}_{\mathrm{eff}}$ is <1σ consistent with literature determinations using high-resolution spectra (5600–5700 K; Baumann et al. 2010; Ramírez et al. 2012; Luck 2017), and all parameters match our model interpolation below.

Zoom In Zoom Out Reset image size

Version 8 of the TESS Input Catalog (TICv8; Stassun et al. 2018, 2019) lists stellar parameters of ${T}_{\mathrm{eff}}$ = 5693 ± 153 K, [Fe/H] = 0.017 ± 0.017, R_* = 0.903 ± 0.055 R_⊙, and L_* = 0.772 ± 0.020 L_⊙. These are all in ≃1σ agreement with our SED-based parameters.

We derived spectral parameters from the TRES spectra of HD 63433 using the Spectral Parameter Classification (SPC) tool (Buchhave et al. 2012). SPC cross-correlates the observed spectrum against a grid of synthetic spectra based on Kurucz atmospheric models (Kurucz 1993). ${T}_{\mathrm{eff}}$, log g, bulk metallicity ([M/H]), and $v\sin {i}_{* }$ are allowed to vary as free parameters. This yielded ${T}_{\mathrm{eff}}$ = 5705 ± 50 K, $\mathrm{log}\,g$  = 4.59 ± 0.10, and [M/H] = −0.09 ± 0.08.

We ran a similar analysis using the HARPS-N stacked spectrum with ARES/MOOG following Sousa et al. (2015). Including empirical corrections from Sousa et al. (2011) and Mortier et al. (2014) yielded ${T}_{\mathrm{eff}}$ = 5764 ± 73 K, $\mathrm{log}\,g$  = 4.65 ± 0.12, and [Fe/H] = 0.04 ± 0.05.

Both methods were consistent with our ${T}_{\mathrm{eff}}$ derived using the SED, and the metallicity estimate agrees with the established value for UMaG (−0.03 ± 0.05; Ammler-von Eiff & Guenther 2009).

To determine the mass of HD 63433, we used Mesa Isochrones and Stellar Tracks (MIST; Choi et al. 2016). We compared all available photometry to the model-predicted values, accounting for errors in both the photometric zero-points and stellar variability as in Section 3.2. We restrict the comparison to 300–600 Myr and solar metallicity based on the properties of the cluster. We assumed Gaussian errors on the magnitudes but included a free parameter to describe underestimated uncertainties in the models or data. The best-fit parameters from the MIST models were M_* = 0.991 ± 0.027 M_⊙, R_* = 0.895 ± 0.021 R_⊙, ${T}_{\mathrm{eff}}$ = 5690 ± 61 K, and L_* = 0.784 ± 0.031 L_⊙. These are consistent with our other determinations, but we adopt our empirical L_*, ${T}_{\mathrm{eff}}$, and R_* estimates from the SED and only utilize the M_* value from the evolutionary models in our analysis.

Using the combination of projected rotation velocity ($v\sin {i}_{* }$), P_rot, and R_*, we can estimate the stellar inclination (i_*) and hence test whether the stellar spin and planetary orbit are consistent with alignment. In principle, this is done by estimating the V term in $v\sin {i}_{* }$ using V = 2π R_*/P_rot, although in practice it requires additional statistical corrections, including the fact that we can only measure alignment projected onto the sky. To this end, we follow the formalism from Masuda & Winn (2020), which handles the hard barrier at i_* > 90° by rewriting the relation in terms of $\cos (i)$.

We used our P_rot from Section 3.1 estimated from the TESS light curve and our R_* derived in Section 3.2. HD 63433 has $v\sin {i}_{* }$ measurements from a range of literature sources, with estimates from 7.0 km s⁻¹ (Marsden et al. 2014) to 7.7 km s⁻¹ (Luck 2017). Our fit of the TRES spectra yielded a consistent estimate of $v\sin {i}_{* }$  = 7.3 ± 0.3 km s⁻¹ with a macroturbulent velocity of 4.2 ± 1.2 km s⁻¹, and the SpecMatch run on the NRES spectra yielded 7.1 ± 0.3 km s⁻¹. We adopted 7.3 ± 0.3 km s⁻¹, which encompassed all estimates.

The combined parameters yielded an equatorial velocity (V) of 7.16 ± 0.29 km s⁻¹ and a lower limit for the inclination of i_* > 74° at 68% confidence and i > 58° at 95% confidence. This is consistent with the stellar rotation being aligned with the planetary orbits (i ≃ 90°).

HD 63433 has adaptive optics or interferometric data spanning almost a decade, from 1999 (Mason et al. 2001) to 2008 (Raghavan et al. 2012). The deepest extant high-resolution imaging reported in the literature for HD 63433 was obtained with the NaCo instrument at the Very Large Telescope on 2004 January 16 UT (Program 072.C-0485(A); PI Ammler). The observation consisted of a series of individual 35 s exposures, totaling 980 s in all, taken with the K_s filter and with the central star behind a 07 opaque Lyot coronagraph. The results of these observations were reported by Ammler-von Eiff et al. (2016), who found no candidate companions within ρ < 9''. The detection limits were reported in a figure in that work and achieved contrasts of ΔK_s ∼ 7 mag at 05, ΔK_s ∼ 9 mag at 1'', ΔK_s ∼ 11 mag at 2'', and ΔK_s ∼ 13 mag at ≥3''. Given an age of τ ∼ 400 Myr, the evolutionary models of Baraffe et al. (2015) would imply corresponding physical limits of M < 65M_Jup at ρ ∼ 11 au and M < 50M_Jup at ρ > 22 au.

Mason et al. (2001) and Raghavan et al. (2012) reported null detections at higher spatial resolution using speckle and long-baseline interferometry. These observations are consistent with the limits set by the lack of Gaia excess noise as indicated by the Renormalized Unit Weight Error (Lindegren et al. 2018). HD 63433 has RUWE = 0.98, consistent with the distribution of values seen for single stars. Based on a calibration of the companion parameter space that does induce excess noise, this corresponds to contrast limits of ΔG ∼ 0 mag at ρ = 30 mas, ΔG ∼ 4 mag at ρ = 80 mas, and ΔG ∼ 5 mag at ρ ≥ 200 mas. The same evolutionary models would imply corresponding physical limits for equal-mass companions at ρ ∼ 0.7 au, M < 0.4 M_⊙ at ρ ∼ 1.8 au, and M < 0.3 M_⊙ at ρ > 4.4 au.

Finally, the Gaia DR2 catalog (Lindegren et al. 2018) did not report any comoving, codistant companions within <1° of HD 63433. Oh et al. (2017) reported a comoving companion based on Gaia DR1 astrometry, but the DR2 parameters for the two stars are inconsistent with each other, and the claimed companion is >3° away from HD 63433. Given the Gaia catalog's completeness limit of G ∼ 20.5 mag at moderately high galactic latitudes and sensitivity at ρ > 3'' (Ziegler et al. 2018; Brandeker & Cataldi 2019), the absence of wide companions corresponds to a physical limit of M < 0.05 M_⊙ at ρ > 66 au to ρ  ∼ 80,000 au.

We identified a 7-day period planet candidate during a visual survey of the 2-minute cadence TESS Candidate Target List data (Stassun et al. 2018) via the light-curve-examining software LcTools⁴⁸ (Schmitt et al. 2019), which initially was introduced by Kipping et al. (2015). Further inspection revealed two additional transits of similar depth and duration located at 1844.057760 and 1864.606371 TBJD (BJD −2,457,000), indicating the presence of a second planet candidate with a period of approximately 21 days. The first candidate was released as a TESS object of interest (TOI) from analysis of the SPOC light curve, the first on 2020 February 19 and the second on 2020 February 20, nearly simultaneous with our visual search. The SPOC data validation reports (Twicken et al. 2018; Li et al. 2019) note a significant centroid offset for the outer planet, but the offset is not consistent with any nearby source, and such offsets are common for young variable stars (Newton et al. 2019). This offset was likely due to saturation (expected at T = 6.27), which degrades the centroids in the row direction.

We searched for additional planets using the notch filter, as described in Rizzuto et al. (2017). We recovered the two planets identified above but found no additional planet signals that passed all checks. Instead, we set limits on the existence of additional planets using an injection/recovery test, again following Rizzuto et al. (2017). To briefly summarize, we generated planets using the BATMAN package following a uniform distribution in period, b, and orbital phase. Half of the planet radii are drawn from a uniform distribution and the other half from a β distribution with coefficients α = 2 and β = 6. We used this mixed distribution to ensure higher sampling around smaller and more common planets. We then detrended the light curve using the notch filter and searched for planets in the detrended curve using a box least squares algorithm (Kovács et al. 2002). The results are summarized in Figure 5. We found that our search would be sensitive to R_P ≃  1 R_⊕ planets at periods of <5 days and ${R}_{P}\simeq 2\,{R}_{\oplus }$ out to 15 days. We required at least two transits to consider a signal recovered, but the light curve covers <30 days, so most injected planets with >15-day periods were not recovered.

Zoom In Zoom Out Reset image size

We fit the TESS photometry simultaneously with the HARPS-N velocities during transit (of the RM effect) using the misttborn (MCMC Interface for Synthesis of Transits, Tomography, Binaries, and Others of a Relevant Nature) fitting code⁴⁹ first described in Mann et al. (2016a) and expanded on in Johnson et al. (2018). misttborn uses BATMAN (Kreidberg 2015) to generate model light curves and emcee (Foreman-Mackey et al. 2013) to explore the transit parameter space using an affine-invariant Markov Chain Monte Carlo (MCMC) algorithm. We did not include any of the other radial velocities in this analysis, especially given the complication of stellar activity, as they are not precise enough to detect the reflex motion due to these small planets.

The standard implementation of misttborn fits for six parameters for each transiting planet—time of periastron (T₀), orbital period of the planet (P), planet-to-star radius ratio (R_p/R_⋆), impact parameter (b), and two parameters ($\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$) to characterize the orbital eccentricity (e) and argument of periastron (ω)—as well as three parameters related to the star: the stellar density (ρ_⋆) and the linear and quadratic limb-darkening coefficients (q₁, q₂) following the triangular sampling prescription of Kipping (2013).

To model stellar variations, misttborn includes a Gaussian process (GP) regression module, utilizing the celerite code (Foreman-Mackey et al. 2017). For the GP kernel, we mostly followed Foreman-Mackey et al. (2017) and used a mixture of two stochastically driven damped simple harmonic oscillators (SHOs), at periods P_GP (primary) and 0.5P_GP (secondary). In addition to the stellar rotation period (ln (P_GP)), the light-curve kernel is characterized by a variability amplitude of the fundamental signal (ln A₁), the decay timescale for the secondary signal (ln Q₂, the quality factor), the difference between the quality factor of the first and second signal ($\mathrm{ln}{\rm{\Delta }}Q={Q}_{1}-{Q}_{2}$), and a mix parameter (Mix) that describes the relative contribution of the two SHOs (where ${A}_{1}/{A}_{2}=1+{e}^{-\mathrm{Mix}}$).

For the RM data, we used misttborn to fit for eight additional parameters for planet b. The primary parameters were the sky-projected spin–orbit misalignment (λ), the stellar rotation broadening ($v\sin {i}_{* }$), and the intrinsic width of the Gaussian line profile of individual surface elements (v_int), which approximates the combined effects of thermal, microturbulent, and macroturbulent broadening. We also included a quadratic polynomial fit to the out-of-transit variations in the radial velocity data (γ, $\dot{\gamma }$, and $\ddot{\gamma }$). We used a generic polynomial because the overall slope in the radial velocity curve is likely dominated by stellar activity, rather than a predictable sinusoidal curve induced by the planets. The last two parameters are the two limb-darkening coefficients (${q}_{1,\mathrm{RM}}$ and ${q}_{2,\mathrm{RM}}$). These RM limb-darkening parameters were fit separately from those for the photometry because of differences in the HARPS-N and TESS wavelength coverage. From these parameters, we produced an analytic RM model following the methodology of Hirano et al. (2011) and Addison et al. (2013). This model consists of an analytic function of v sin i and v_int, multiplied by the flux drop due to the transiting planet; we calculated the flux decrement needed for this model with BATMAN, following the same methodology as for the photometric light curves. We note that the GP described above is not used for the RM data, only the photometric curve.

We ran two separate MCMC chains, the first as described above, and the second with e and ω locked at 0. For both chains, we ran the MCMC using 100 walkers for 250,000 steps, including a burn-in of 20,000 steps. The autocorrelation time indicated that this was sufficient for convergence. We also applied Gaussian priors on the limb-darkening coefficients (for both TESS and the HARPS-N data) based on the values in Claret & Bloemen (2011) and Parviainen & Aigrain (2015), with errors accounting for the difference between these two estimates (which differ by 0.05–0.07). For $v\sin {i}_{* }$ and v_int, we used Gaussian priors of 7.3 ± 0.3 km s⁻¹ and 4.2 ± 1.2 km s⁻¹ based on analysis from Section 3.5 and the investigation of Doyle et al. (2014). For the fit with e = 0, we applied Gaussian priors for the stellar density taken from our derived stellar parameters derived in Section 3.2. All other parameters were sampled uniformly with physically motivated boundaries: $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$ were restricted to (−1, 1), $| b| \lt 1+{R}_{p}/{R}_{s}$, and T₀ to the time period sampled by the data. The GP mix parameter was restricted to be between −10 and 10. We let the linear and quadratic terms of the radial velocity curve in the RM data float, as this is produced by some combination of actual reflex motion of the star due to the two planets and stellar activity of this young, relatively rapidly rotating star. The full list of fit parameters, priors, and imposed limits is given in Table 5.

Table 5.  Parameters and Priors

Parameter	Prior
	Planet b	Planet c
T0 (TJD)a	${ \mathcal U }[1916.4,1916.5]$	${ \mathcal U }[1842,1860]$
P (days)	${ \mathcal U }[0,15]$	${ \mathcal U }[15,30]$
${R}_{P}/{R}_{\star }$	${ \mathcal U }[0,1]$	${ \mathcal U }[0,1]$
b	${ \mathcal U }[| b| \lt 1+{R}_{P}/{R}_{* }]$	${ \mathcal U }[| b| \lt 1+{R}_{P}/{R}_{* }]$
${\rho }_{\star }$ (${\rho }_{\odot }$)	${ \mathcal N }[1.30,0.15]$
${q}_{\mathrm{1,1}}$	${ \mathcal N }[0.30,0.06]$
${q}_{\mathrm{2,1}}$	${ \mathcal N }[0.37,0.05]$
$\sqrt{e}\sin \omega$	${ \mathcal U }[-1,1]$	${ \mathcal U }[-1,1]$
$\sqrt{e}\cos \omega$	${ \mathcal U }[-1,1]$	${ \mathcal U }[-1,1]$
$v\sin {i}_{* }$(km s−1)	${ \mathcal N }[7.3,0.3]$	...
λ (deg)	${ \mathcal U }[-180,180]$	...
${q}_{1,\mathrm{RM}}$	${ \mathcal N }[0.53,0.08]$	...
${q}_{2,\mathrm{RM}}$	${ \mathcal N }[0.39,0.06]$	...
${v}_{\mathrm{int}}$ (km s−1)	${ \mathcal N }[4.2\,1.2]$	...
γ1 (km s−1)	${ \mathcal U }[-17,-14]$	...
$\dot{\gamma }$ (km s−1)	${ \mathcal U }[-1,1]$	...
$\ddot{\gamma }$ (km s−1)	${ \mathcal U }[-1,1]$	...
$\mathrm{ln}{P}_{\mathrm{GP}}$	${ \mathcal U }[1,2]$
$\mathrm{ln}{A}_{1}$	${ \mathcal U }[-\infty ,1]$
$\mathrm{ln}({Q}_{0})$	${ \mathcal U }[0.5,\infty ]$
$\mathrm{ln}{\rm{\Delta }}Q$	${ \mathcal U }[0,\infty ]$
Mix	${ \mathcal U }[-10,10]$

Notes. ${ \mathcal U }[X,Y]$ denotes a uniform prior limited to between X and Y, and ${ \mathcal N }[X,Y]$ denotes a Gaussian prior with mean X and standard deviation Y.

^aIt is standard to report T₀ as the midtransit point for the first transit. However, for computational reasons in the RM fit, we restrict T₀ around the RM observations for planet b.

Download table as:  ASCIITypeset image

The resulting fit light curve is shown in Figure 6 with the velocity curve in Figure 7. The best-fitting model and derived parameters, along with 68% credible intervals, are listed in Table 6. Figures 8(a) and b show partial posteriors for the MCMC fit parameters.⁵⁰

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  Transit Fit Parameters

Parameter	Planet b	Planet c	Planet b	Planet c
	e, ω Fixed		e, ω Free
Transit Fit Parameters
T0 (TJD)	${1916.4526}_{-0.0027}^{+0.0032}$	${1844.05799}_{-0.00087}^{+0.00084}$	${1916.4533}_{-0.0027}^{+0.0037}$	${1844.05791}_{-0.00076}^{+0.0008}$
P (days)	${7.10793}_{-0.00034}^{+0.0004}$	${20.5453}_{-0.0013}^{+0.0012}$	${7.10801}_{-0.00034}^{+0.00046}$	20.5455 ± 0.0011
${R}_{P}/{R}_{\star }$	0.02161 ± 0.00055	0.02687 ± 0.0007	${0.02168}_{-0.00058}^{+0.00065}$	${0.02637}_{-0.00074}^{+0.00077}$
b	${0.18}_{-0.13}^{+0.17}$	${0.512}_{-0.033}^{+0.063}$	${0.26}_{-0.17}^{+0.19}$	${0.29}_{-0.2}^{+0.21}$
${\rho }_{\star }$ (${\rho }_{\odot }$)	${1.293}_{-0.2}^{+0.079}$		1.3 ± 0.14
${q}_{\mathrm{1,1}}$	${0.302}_{-0.057}^{+0.058}$		${0.299}_{-0.055}^{+0.06}$
${q}_{\mathrm{2,1}}$	0.369 ± 0.048		0.372 ± 0.048
$\sqrt{e}\sin \omega$	0 (fixed)	0 (fixed)	$-{0.08}_{-0.13}^{+0.11}$	${0.09}_{-0.18}^{+0.13}$
$\sqrt{e}\cos \omega$	0 (fixed)	0 (fixed)	${0.01}_{-0.37}^{+0.35}$	${0.09}_{-0.43}^{+0.39}$
RM Parameters
$v\sin {i}_{\star }$ (km s−1)	7.28 ± 0.29	⋯	${7.3}_{-0.3}^{+0.29}$	⋯
λ (deg)	${1.0}_{-43.0}^{+41.0}$	0 (fixed)	${8.0}_{-45.0}^{+33.0}$	0 (fixed)
${q}_{1,\mathrm{RM}}$	${0.534}_{-0.081}^{+0.08}$		${0.524}_{-0.078}^{+0.084}$
${q}_{2,\mathrm{RM}}$	${0.388}_{-0.059}^{+0.056}$		${0.388}_{-0.059}^{+0.057}$
${v}_{\mathrm{int}}$ (km s−1)	4.3 ± 1.1	⋯	4.2 ± 1.2	⋯
${\gamma }_{1}$ (km s−1)	$-{15.74559}_{-0.00064}^{+0.0007}$	⋯	$-{15.74542}_{-0.0007}^{+0.00075}$	⋯
$\dot{\gamma }$ (km s−1 day−1)	${0.0144}_{-0.0048}^{+0.0049}$	⋯	${0.014}_{-0.0049}^{+0.0048}$	⋯
$\ddot{\gamma }$ (km s−1 day−2)	−0.184 ± 0.078	⋯	$-{0.198}_{-0.083}^{+0.082}$	⋯
Gaussian Process Parameters
$\mathrm{log}({P}_{\mathrm{GP}})$	${1.872}_{-0.017}^{+0.018}$		1.872 ± 0.017
$\mathrm{log}({A}_{1})$	$-{8.5}_{-1.2}^{+1.7}$		$-{8.6}_{-1.1}^{+1.6}$
$\mathrm{log}({\rm{\Delta }}Q)$	${1.6}_{-1.1}^{+1.9}$		${1.6}_{-1.2}^{+1.9}$
$\mathrm{log}({Q}_{0})$	${3.07}_{-0.85}^{+1.42}$		${3.03}_{-0.82}^{+1.3}$
Mix	$-{3.5}_{-1.9}^{+1.3}$		$-{3.3}_{-1.9}^{+1.3}$
Derived Parameters
$a/{R}_{\star }$	${16.95}_{-0.82}^{+0.34}$	${34.38}_{-2.0}^{+0.69}$	${16.75}_{-0.74}^{+0.47}$	${36.1}_{-1.7}^{+1.1}$
i (deg)	${89.38}_{-0.64}^{+0.43}$	${89.147}_{-0.2}^{+0.069}$	${89.1}_{-0.69}^{+0.59}$	${89.51}_{-0.35}^{+0.34}$
δ (%)	${0.0467}_{-0.0023}^{+0.0024}$	${0.0722}_{-0.0037}^{+0.0038}$	${0.047}_{-0.0025}^{+0.0029}$	${0.0696}_{-0.0038}^{+0.0041}$
T14 (days)	${0.134}_{-0.0013}^{+0.0014}$	${0.1695}_{-0.0013}^{+0.0015}$	${0.133}_{-0.02}^{+0.018}$	0.17 ± 0.031
T23 (days)	${0.1279}_{-0.0012}^{+0.0013}$	${0.1573}_{-0.0015}^{+0.0014}$	${0.127}_{-0.02}^{+0.016}$	0.159 ± 0.029
${g}_{\mathrm{1,1}}$	${0.402}_{-0.06}^{+0.063}$		${0.403}_{-0.059}^{+0.062}$
${g}_{\mathrm{2,1}}$	${0.143}_{-0.054}^{+0.058}$		${0.139}_{-0.054}^{+0.058}$
${g}_{1,\mathrm{RM}}$	${0.562}_{-0.093}^{+0.094}$		${0.558}_{-0.094}^{+0.091}$
${g}_{2,\mathrm{RM}}$	${0.162}_{-0.081}^{+0.088}$		${0.161}_{-0.082}^{+0.087}$
e	0 (fixed)	0 (fixed)	${0.085}_{-0.067}^{+0.179}$	${0.114}_{-0.068}^{+0.204}$
Rp (${R}_{\oplus }$)	2.15 ± 0.10	2.67 ± 0.12	2.15 ± 0.10	2.64 ± 0.12
a (au)a	${0.0719}_{-0.0044}^{+0.0031}$	${0.1458}_{-0.0101}^{+0.0062}$	${0.0710}_{-0.0041}^{+0.0033}$	${0.1531}_{-0.0092}^{+0.0074}$

Note.

^aR_p derived using the R_* value from Table 3.

Download table as:  ASCIITypeset image

The GP fit accurately reproduced the overall out-of-transit variability (Figure 6). Similarly, the resulting period matched the rotation period from Section 3.1 (6.54 days vs. 6.45 days) and agrees with the predicted value for the star's mass and membership to UMaG (6.9 ± 0.4 days).

The transit duration suggests a low eccentricity for both planets (e < 0.2), as is common for multitransiting systems (Van Eylen & Albrecht 2015). Our analysis did not include any correction for biases in the eccentricity distribution of planets (Kipping 2014), but accounting for this would only drive the resulting eccentricity values down. Consistent with this, the ρ_* derived from the transit assuming e = 0 and a uniform density prior yields ${\rho }_{* }={1.331}_{-0.1}^{+0.056}{\rho }_{\odot }$, in excellent agreement with our derived value from Section 3 (${\rho }_{* }=1.3\pm 0.15{\rho }_{\odot }$). Thus, if we were to assume that the eccentricities are ≃0, we can consider this an additional verification of our adopted stellar parameters.

Due to the relatively low amplitude of the RM effect for HD 63433 b and the complication of stellar activity, our posterior for the sky-projected spin–orbit misalignment λ is broad. However, we clearly demonstrated that the planetary orbit is prograde; retrograde orbits would yield $| \lambda | \gt 90^\circ$, which is completely ruled out (Figure 8). We discuss the implications of this measurement in more detail in Section 7. Furthermore, we clearly detected the RM effect due to the transit of HD 63433 b (the signal is inconsistent with no RM signal), confirming that the planet is real.

As an additional test on our RM fit, we ran an MCMC chain using a linear trend in velocity rather than a second-order polynomial (i.e., $\ddot{\gamma }$ fixed at 0). The resulting λ posterior was not significantly different (λ = 70 ± 35°). The second-order polynomial was preferred statistically (ΔBIC = 5), so we use it for all results reported in Table 6.

We compared the wealth of radial velocity data (Table 1) to the predicted velocity curve of a planet or binary at the orbital period of the outer planet (20.5 days). We assumed a low eccentricity (e < 0.1) and sampled over the whole range of ω and mass ratios. We included a zero-point correction term between instruments, which takes the value preferred by the predicted/model velocity curve. Including this term meant that two to three epochs each from TRES and NRES provided little information and were not included. To account for stellar jitter and instrumental drift, we inflated errors in the velocities based on the scatter between points for each instrument (37 m s⁻¹ for SOPHIE/ELODIE and 24 m s⁻¹ for Lick).

The velocities are not precise enough to detect either planet but easily rule out (at 99.7%) any stellar or planetary companion down to ≃Jupiter mass at the orbital period of planet c (Figure 9).

Zoom In Zoom Out Reset image size

As detailed in Vanderburg et al. (2019), if the observed transits are due to blends from a background eclipsing/transiting system, the true radius ratio can be determined from the ratio of the ingress time (T₁₂ or T₃₄) to the time between first and third contact (T₁₃). This provides a constraint on the brightest possible background source that could produce the observed transit depth: ${\rm{\Delta }}{m}_{\mathrm{TESS}}\leqslant 2.5{\mathrm{log}}_{10}({T}_{12}^{2}/{T}_{13}^{2}/\delta )$, where δ is the transit depth. Using our results from Section 4, we find Δm < 2.4 mag at 99.7% confidence. The combination of AO imaging and Gaia DR2 (Section 3.6) rules out any such background star down to <80 mas and out to 45.

We also rule out a background eclipsing binary behind HD 63433. The high-resolution imaging spans 9 yr (1999.16–2008.28) and is sensitive to companions brighter than our magnitude threshold down to 80 mas. Due to its proper motion, HD 63433 has moved more than 100 mas over the same time period; thus, any foreground or background star not visible in the earliest data set would be visible in the final observation.

To explore the range of possible stellar companions, we used a Monte Carlo simulation of 5 × 10⁶ binaries, comparing each generated system to the velocities, high-resolution imaging, and limits from Gaia imaging and astrometry. Companions were generated following a lognormal distribution in period following Raghavan et al. (2010), but uniform in other orbital parameters. The radial velocity data listed in Table 1 span more than 22 yr, which overlaps in parameter space with the high-contrast imaging data (down to ≃1 au) and Gaia constraints. A negligible fraction (<0.01%) of generated companions are consistent with the external constraints, resolved in TESS, and reproduce the observed transit depth, statistically ruling out this scenario.

Thanks to HD 63433's brightness (V = 6.9, K = 5.3; see Figure 10), this system is ideal for a variety of follow-up observations to characterize the planets and the system as a whole. Observations over the coming years will allow us to determine the system's 3D architecture, measure mass loss from the planets and study their atmospheres, and potentially measure the masses of the planets.

The HARPS-N observations of the b transit demonstrate that HD 63433 is well suited for additional RM observations. Repeat observations of the planet b would enable more detailed accounting of stellar variability (Montet et al. 2020; Zhou et al. 2020) and provide more robust constraints on λ. Observations of the c transit would both confirm the planet and allow a measurement of the mutual inclination between the orbits of the two planets.

One of the first discoveries from the young planet population has been that young planets are statistically larger than their older counterparts (e.g., Mann et al. 2018; Rizzuto et al. 2018). This offset could be explained by thermal contraction of an H/He-dominated atmosphere (Lopez & Fortney 2014), atmospheric mass loss from interactions with the (still active) host star (Murray-Clay et al. 2009), or photochemical hazes making the atmosphere larger and puffier (Gao & Zhang 2020). Right now, the difference is only an offset in the planet radius distribution with age, making it difficult to distinguish between these scenarios. Instead, planet masses (and hence densities) are needed. While challenging, both of HD 63433's planets may be within reach of existing PRV spectrographs. Assuming masses of 5.5 and 7.3 M_⊕ for planets b and c, respectively (Chen & Kipping 2017), the predicted radial velocity amplitudes are ≃2 m s⁻¹. This signal is within the reach of existing instruments, but still much smaller than the estimated stellar jitter (20–30 m s⁻¹; Figure 9 and Table 1). Planet b is especially challenging given the similarity of its orbital period to the stellar rotation period (7.11 days vs. 6.45 days). However, a focused campaign designed to separate planetary and stellar signals, as was done for the young system K2-100 (Barragán et al. 2019), will likely yield a mass constraint for planet c.

Wang & Dai (2019) and Gao & Zhang (2020) argue that young planets are likely to have flat transmission spectra owing to either dust or photochemical hazes. There is some evidence to support this from transmission spectroscopy follow-up of young systems (Libby-Roberts et al. 2020; Thao et al. 2020). However, a wider set of observations are required to explore under what conditions young atmospheres are dominated by hazes, dust, and/or clouds. Because the host is bright (H ≃ 5), both planets are well within reach of transmission spectroscopy with the Hubble Space Telescope or James Webb Space Telescope.

Both of the planets lie on the large-radius, gas-rich side of the radius valley (e.g., Owen & Wu 2013; Fulton & Petigura 2018). Given the young age of the system, it is likely that both planets are actively losing their atmospheres through photoevaporation (e.g., Owen & Jackson 2012; Lopez & Fortney 2013) or core-powered mass loss (Ginzburg et al. 2018). Given the X-ray flux of HD 63433 observed by XMM-Newton as a part of its slew catalog (94 erg s⁻¹ cm²) and the energy-limited mass-loss relation (e.g., Owen 2019), we estimate mass-loss rates of $\eta \approx 2.79\times {10}^{11}$ g s⁻¹ and $\eta \approx 7.07\times {10}^{10}$ g s⁻¹ for planets b and c, respectively, where η describes the heating efficiency of the atmospheres. This is higher than many other planets of similar size, including Gl 436b and GJ 3470b, both of which have detected exospheres (Ehrenreich et al. 2015; Ninan et al. 2020).