# Hazy with a chance of star spots:
constraining the atmosphere of the young planet, K2-33b

Pa Chia Thao NSF Graduate Research Fellow  Jack Kent Cooke Foundation Graduate
Scholar  Department of Physics and Astronomy, The University of North Carolina
at Chapel Hill, Chapel Hill, NC 27599, USA  Andrew W. Mann Department of
Physics and Astronomy, The University of North Carolina at Chapel Hill, Chapel
Hill, NC 27599, USA  Peter Gao Earth and Planets Laboratory, Carnegie
Institution for Science, 5241 Broad Branch Road, NW, Washington, DC 20015, USA
Dylan A. Owens Department of Physics and Astronomy, The University of North
Carolina at Chapel Hill, Chapel Hill, NC 27599, USA  Gemini Observatory, NSF’s
NOIRLab, 670 N. A’ohoku Place Hilo, Hawaii, 96720, USA  Andrew Vanderburg
Department of Physics and Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Elisabeth R.
Newton Department of Physics and Astronomy, Dartmouth College, Hanover, NH
03755, USA  Yao Tang  University of California, Santa Cruz, Santa Cruz, CA
95064, USA  Matthew J. Fields Department of Physics and Astronomy, The
University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA
Trevor J. David Center for Computational Astrophysics, Flatiron Institute, New
York, NY 10010, USA  American Museum of Natural History, Central Park West,
New York, NY 10024, USA  Jonathan M. Irwin  Institute of Astronomy, University
of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  Tim-Oliver Husser
Institut für Astrophysik, Georg-August-Universität Göttingen, Friedrich-Hund-
Platz 1, 37077 Göttingen, Germany  David Charbonneau Harvard-Smithsonian
Center for Astrophysics, Harvard University, Cambridge, MA 02138, USA  Sarah
Ballard Bryant Space Science Center, Department of Astronomy, University of
Florida, Gainesville, FL 32611, USA

###### Abstract

Although all-sky surveys have led to the discovery of dozens of young planets,
little is known about their atmospheres. Here, we present multi-wavelength
transit data for the super Neptune-sized exoplanet, K2-33b – the youngest
($\sim$10 Myr) transiting exoplanet to-date. We combined photometric
observations of K2-33 covering a total of 33 transits spanning $>$2 years,
taken from K2, MEarth, HST, and Spitzer. The transit photometry spanned from
the optical to the near-infrared (0.6-4.5µm), enabling us to construct a
transmission spectrum of the planet. We find that the optical transit depths
are nearly a factor of two deeper than those from the near-infrared. This
difference holds across multiple datasets taken over years, ruling out issues
of data analysis and unconstrained systematics. Surface inhomogeneities on the
young star can reproduce some of the difference, but required spot coverage
fractions ($>$60%) are ruled out by the observed stellar spectrum ($<$20%). We
find a better fit to the transmission spectrum using photochemical hazes,
which were predicted to be strong in young, moderate-temperature, and large-
radius planets like K2-33b. A tholin haze with CO as the dominant gaseous
carbon carrier in the atmosphere can reasonably reproduce the data with small
or no stellar surface inhomogeneities, consistent with the stellar spectrum.
The HST data quality is insufficient for the detection of any molecular
features. More observations would be required to fully characterize the hazes
and spot properties and confirm the presence of CO suggested by current data.

Exoplanet atmospheres; Transit photometry; Open star clusters; Exoplanet
evolution; Starspots; M dwarf stars; Markov chain Monte Carlo; Light curves;
Starspots

††software: LDTK (Parviainen & Aigrain, 2015), Iraclis , POET \+ BLISS
(Stevenson et al., 2012), MISTTBORN , emcee (Foreman-Mackey et al., 2013) ,
corner.py (Foreman-Mackey, 2016), celerite (Foreman-Mackey et al., 2017),
matplotlib (Hunter, 2007), batman (Kreidberg, 2015), Fleck (Morris, 2020),
Astropy (Astropy Collaboration et al., 2013, 2018), numpy (Harris et al.,
2020)††facilities: Spitzer (IRAC), K2 , MEarth , Hubble (WFC3)

##  1 Introduction

The study of exoplanet atmospheres has progressed at a tremendous pace in the
last decade due to the high-sensitivity spectroscopic observations covering
secondary eclipse spectroscopy, high-resolution Doppler spectroscopy, and
direct imaging spectroscopy. Transmission spectroscopy, a measurement of the
planet radius as a function of wavelength, has enabled the detection of a
range of atomic and molecular species (e.g., Seager & Sasselov, 2000; Brown,
2001; Hubbard et al., 2001; Huitson et al., 2013; Sing et al., 2015, 2016).
Detecting such molecules can provide information about the planet’s formation
and migration history, although complexities in the planet and disk properties
can make such detections difficult to interpret. For example, planets formed
beyond the snow line will have an atmosphere that is carbon rich (an enhanced
C/O) compared to planets formed closer-in (Öberg et al., 2011; Booth et al.,
2017; Booth & Ilee, 2019).

Measurements of young exoplanet atmospheres offer a more direct correlation to
the disk, which is unadulterated by later-stage evolution. However, the
majority of planets that have had their atmosphere characterized are mature
($>$1 Gyr), or have unconstrained ages. This is because young planetary
systems with precise ages are rare (usually limited to those in young
associations, e.g., Mann et al., 2016a; Newton et al., 2019; Bouma et al.,
2020; David et al., 2019; Obermeier et al., 2016) and are harder to detect due
to stellar variability. Only in the last decade has there been a significant
increase in the number of known young planetary systems. This has been driven
by improvements in near-infrared radial velocities (e.g., Johns-Krull et al.,
2016), mitigation of stellar variability (e.g., Donati et al., 2017; Rizzuto
et al., 2017), as well as, NASA’s K2 (Howard et al., 2014) and TESS
(Transiting Exoplanet Survey Satellite; Ricker, 2014) missions surveying
nearby young clusters and associations.

Only three young ($<$1 Gyr) exoplanets have had their atmospheres
characterized through transmission spectroscopy: 111While there are several
young planets ($\beta$ Pictoris b, 51 Eridani b, HR 8799 b; Chilcote et al.,
2017; Samland et al., 2017; Lee et al., 2013) with direct spectroscopy, these
are more massive and at wider separations ($\gg$ 1 AU) than the planets
transits survey; they likely have a different formation history. K2-25b,
Kepler 51b, and Kepler 51d – all of which are adolescent ($\approx$ 500-700
Myr) and may be in the midst of undergoing thermal contraction. Statistical
analyses of Kepler planet candidates suggest that most planets larger than 1.6
R⊕ – as is the case for the aforementioned planets – have an extended envelope
with a low mean molecular weight (e.g., of H/He) (Rogers, 2015), and therefore
were predicted to have a large scale height, resulting in strong spectral
features. However, all three planets showed a featureless transmission spectra
(Thao et al., 2020; Libby-Roberts et al., 2020), which were interpreted as
ongoing atmospheric mass loss (expected in such young systems) leading to
dusty atmospheric outflows with small dust grains (Wang & Dai, 2019) or high-
altitude photochemical hazes (Gao & Zhang, 2020). The small particle size
leads to a lower opacity (and hence shallower transit) in the NIR compared to
the optical. Grains or hazes also inflate observed radii (Lammer et al., 2016;
Cubillos et al., 2017) and weaken the spectral features. Gao & Zhang (2020)
showed that this effect is strongest on young, low-mass planets with a
moderate or low equilibrium temperature. Additional objects and broader
wavelength coverage could confirm these predictions.

A reoccurring challenge in interpreting transmission spectra is the
contamination by star spots and faculae (Rackham et al., 2017). Spots can
change the observed signal, whether or not the planet crosses them. In the
case of unocculted spots, the transit chord will be brighter than the stellar
average, leading to a deeper transit with increasing spot coverage and
contrast even for fixed planet size. Since spot contrast depends on
wavelength, such unocculted spots can introduce features and produce haze-like
variations in the transmission spectrum that are not present in the planet.
Barclay et al. (2021) showed that spots alone can explain the observed H2O
features seen in the transmission spectrum of K2-18 b (Benneke et al., 2019;
Tsiaras et al., 2019). The imprint of the star on the planet’s spectrum is
particularly concerning for young planets, as spot coverage fractions are
statistically higher for young stars (Morris, 2020; Luger et al., 2021a).

To understand the evolution and formation of planetary atmospheres, a greater
sample size of young planets – particularly new-born planets – need their
atmospheres characterized. To this end, we explore the atmospheric
transmission spectrum of the youngest transiting exoplanet to-date, K2-33b
(David et al., 2016b; Mann et al., 2016b). This $\simeq$ 10 Myr, super
Neptune-sized exoplanet orbits a late type pre-main sequence M3.5 dwarf in the
Upper Scorpius OB Association every 5.42 days. Its youth, moderate equilibrium
temperature ($\sim$770 K) and abnormally large radius (5$R_{\oplus}$ in the
optical wavelengths; Figure 1) make it an ideal target for a strong haze
detection. We combined 14 transits from the discovery K2 data with 7 transits
taken from the MEarth survey, 10 transits taken from Spitzer, and a partial
transit taken with HST.

The paper is presented as follows: Section 2 describes our observations and
data reduction. Using the precise parallax from Gaia, we update K2-33’s
stellar parameters in Section 3. We utilize this information in our fit to the
transit light curve, as described in Section 4. Due to the planet’s youth, we
discuss the effect of surface inhomogeneities on the transmission spectrum in
Section 5. In Section 6, we investigate photochemical hazes inferred from the
transmission spectrum of K2-33b. In Section 7, we conclude with a brief
summary of our results.

Figure 1: The detected population of all planets from Kepler and K2 as a
function of planet radius ($R_{\oplus}$) and orbital period (days). Young
($<$1 Gyr), transiting planets that are members of clusters or stellar
associations are plotted in circles colored by their age. K2-33b is outlined
as a star and lands at the edge of the sub-Saturn desert. The youngest planets
(yellow) have radii between Neptune and Jupiter, while no mature planets
(purple) are located within this region. K2-33b is likely still evolving and
offers an opportunity to observe an atmosphere in transition. Planet
properties from NASA Exoplanet Science Institute (2020).

##  2 Observation and Data Reduction

We collected 33 total transits of K2-33b obtained from 2014 to 2017 taken by
K2, the MEarth Project, the Spitzer Space Telescope, and the Hubble Space
Telescope (HST). The combined data sets span from the visible to the near-
infrared (0.64-4.5$\mu$m). The details of each dataset are summarized in Table
1.

###  2.1 K2

The K2 light curve covered 14 transits from 2014 August 23 to 2014 November 13
(UT). We used the K2 data with extraction and reduction as described in one of
the the discovery paper (Mann et al., 2016a). To briefly summarize, following
Vanderburg et al. (2016) and Mann et al. (2016a), we derived a correction to
the raw light curve using least-square minimization to simultaneously fit for
low frequency variations from stellar activity, the Kepler flat field
(Vanderburg & Johnson, 2014), and the transits of K2-33b. As shown in
Grunblatt et al. (2016), this simultaneous fit was required to avoid biasing
both the flat-field correction and the resulting transit depth; the bias is
larger for highly variable stars like K2-33.

Unlike Mann et al. (2016a), we did not apply our initial least-squares stellar
variability fit before passing it to our transit analysis (Section 4). We
instead included the final stellar variability fit as part of the MCMC (Markov
chain Monte Carlo) analysis in the form of a periodic GP (Gaussian Process)
kernel (Section 4). While the initial variability fit is consistent with that
derived from our transit analysis, the stellar signal is large enough on
transit timescales that the earlier fit may yield underestimated errors on the
transit parameters.

Figure 2: K2 light curve of K2-33 after removing signals from the pointing
drift of the telescope (grey points). The pink shaded region shows the
3$\sigma$ range of fits from the GP variability fit, while the pink solid line
shows the combined GP and the transit model.

###  2.2 MEarth

Including the two transits discussed in one of the discovery paper, Mann et
al. (2015), we observed a total of nine transits (or partial transits) of
K2-33b using the MEarth-North and/or MEarth-South arrays (Nutzman &
Charbonneau, 2008; Irwin et al., 2015) spanning 2016 February 16 to 2017 June
13 (UT). Due to differences in weather and target visibility, some transits
were observed with one of the MEarth arrays, and some transits were observed
with both MEarth arrays (see Table 1).

MEarth-North includes eight 40-cm telescopes at Fred Lawrence Whipple
Observatory on Mount Hopkins, Arizona, and MEarth-South uses a nearly
identical set of telescopes located at Cerro Tololo Inter-American Observatory
(CTIO) in Chile. All telescopes were equipped with a $2048\times 2048$ pixel
CCD; MEarth-North CCDs had a pixel scale of 0.78″/ pixel and MEarth-South had
a pixel scale of 0.84″/ pixel. For all observations, we used a Schott RG715
filter (Dittmann et al., 2016). All telescopes integrated for 60 s for a
cadence of $\simeq$90 s per telescope.

MEarth data were reduced following the basic methodology from Irwin et al.
(2007) with additional steps detailed in the documentation of the fourth
MEarth data
release222https://www.cfa.harvard.edu/MEarth/DR4/processing/index.html. Since
out-of-transit baseline varied considerably between transits, we did not
perform additional systematics correction to the data (e.g., for stellar
variability), and instead fit such systematics with our GP regression
simultaneously with the transit (see Section 4).

Following reduction, we decided to remove two transits from our analysis,
resulting in only 7 transits. The two transits were removed due to: 1) the
data taken on 2016 July 17 (UT) only had 1.5hr of usable data near ingress; 2)
the data taken on 2017 May 22 (UT) showed two sudden changes in flux during
the transit, which is likely due to the instrument’s meridian flip and a later
halt in data from MEarth-South, or a poorly-timed flare. For the latter
transit, it was difficult to remove or mitigate these problems because this
transit did not contain an egress and our GP model could not simultaneously
fit for such effects.

###  2.3 Spitzer

We obtained 10 full transits of K2-25b, five in 3.6$\mu$m (Channel 1) and five
in 4.5$\mu$m (Channel 2), taken by the InfraRed Array Camera (IRAC; Fazio et
al., 2004) on the Spitzer Space Telescope (Werner et al., 2004). Observations
were executed over the period of 2016 Nov 08 to 2017 Nov 29 (Program ID:
13037, PI: A. Mann). Two additional transits were also taken in each band for
another program (Program ID: 11026, PI: M. Werner), but were not included in
this analysis. All Spitzer data analyzed here used the $32\times 32$ pixels
sub-array mode, with each image taken as a 2 s exposure. Each transit
consisted of a $\simeq$24 minute dither, a $\simeq$360 minute stare of the
full transit and out-of-transit baseline, followed by another $\simeq$7 minute
dither333https://irachpp.spitzer.caltech.edu/. The initial dither allows an
initial settling time at the new pointing position. For the long stare, we
used the peak-up pointing mode, which keeps the star stable on a
0.5$\times$0.5 pixel box region of the IRAC CCD with relatively uniform
sensitivity (the ‘sweet spot’; Ingalls et al., 2012, 2016).

We first processed the flat-fielded and dark-subtracted Basic Calibrate Data
(BCD) images produced by the Spitzer pipeline using the Photometry for Orbits,
Eccentricities, and Transits (POET; Stevenson et al. 2012; Campo et al. 2011)
444https://github.com/kevin218/POET pipeline to create systematics-corrected
light curves. The process includes masking and flagging bad pixels, and
calculating the Barycentric Julian Dates for each frame. The center position
of the star was fitted using a two-dimensional, elliptical Gaussian in a 15
pixel square window centered on the target’s peak pixel. Simple aperture
photometry was performed using a radius of 2.0 to 4.0 pixels, in increments of
0.25 pixels, with an inner sky annulus set to 7 pixels, and an outer sky
annulus set to 15 pixels. We elected to set the radius to 2.25 pixels for both
channels as it minimizes the standard deviation of the normalized residuals
(SDNR). Additional reduction was done simultaneously with fitting the transit
and is described in Section 4.

###  2.4 Hubble

HST observed a single transit of K2-33b using the Wide Field Camera 3 (WFC3)
on 2017 May 23 (ID: 14887; PI: B. Benneke), spanning six 30-min spacecraft
orbits. Observations used the G141 grating, which provided a wavelength
coverage of 1.1-1.7$\mu$m, and was taken in the 256 $\times$ 256 sub-array
mode to reduce overhead times with the NSAMP=7 and SAMP-SEQ=SPARS25 readout
settings. In the beginning of the observations, a direct (non-dispersed) image
of K2-33 was taken in the F139M Filter for calibration purposes. Afterwards, a
total of 95 spectroscopic images were taken in the forward direction of the
spatial-scanning mode, in which the stellar spectrum is spread along the
spatial direction to avoid non-linearity or saturation. All spectra were taken
with a scan rate of 0.037 arcseconds/second for an effective exposure time of
$\sim$112 seconds to produce a scan 4.38 arcseconds in length.

Initial visual inspection of the spectra indicated that the telescope slewed
while observing K2-33, resulting in some spectra that are slanted or not
visible at all (Figure 3). There are 6 total orbits. After analyzing the
header information, 32 consecutive frames made up of Orbit 5 and Orbit 6
indicated that the guide star acquisition failed and observation was taken
only with the gyro guiding. For this reason, these frames were discarded and
not used in our analysis. Following the standard procedure (Berta et al.,
2012), we excluded the first orbit from our analysis, as it exhibits a much
steeper ramp than the subsequent orbits. Of the remaining three HST orbits,
Orbit 2 provided the out-of-transit baseline, Orbit 3 captured the ingress,
and Orbit 4 captured the mid-transit.

Figure 3: Left:  A typical raw HST image from Orbits 1-4; Right:  A typical
raw image from Orbits 5-6. The telescope slewed so images from these orbits
were discarded.

The data set was reduced and calibrated with the publicly available Iraclis
pipeline (Tsiaras et al., 2016) 555https://github.com/ucl-exoplanets/Iraclis.
Here, we provide a brief description of their method. We first used the _raw
frames provided by the Hubble Legacy Archive
666https://hla.stsci.edu/hlaview.html. The reduction process is performed in
the following order: zero-read and bias-level correction, non-linearity
correction, dark current subtraction, gain conversion, sky background
subtraction, calibration, flat-field correction, and bad pixels and cosmic ray
correction. Iraclis then extracted 1D spectra from the reduced data by
applying a 128 pixel-wide aperture along the dispersion axis and summing along
the cross-dispersion axis for each wavelength. We set the vertical starting
and ending position of the extraction box to 5 pixels above and below the
spectrum in the spatial direction. There is less than a 0.2 pixel shift in
both the x- and y-direction.

The Iraclis extraction resulted in a broadband (white) light curve, covering
the whole wavelength range in the G141 grism (1.088–1.68$\mu$m) and spectral
light curves for each wavelength band. While it is common to adopt a bin width
that yields similar SNR in all bins and/or avoid regions of higher stellar
variability (Mandell et al., 2013), the partial transit yielded relatively low
SNR overall and changes to the bin width did not change any other conclusions
about the system. As a result, we adopted 14 equal-size bins, each with a bin
width of 0.0423$\mu$m.

Table 1: Observation Log  Telescope | Filter/Grism | Transit Number | UT Start Date
---|---|---|---
K2 Campaign 2 | Kepler | 1-14a,ba,bfootnotemark:  | 2014 Aug 03
MEarth North | RG715 | 99b,cb,cfootnotemark:  | 2016 Feb 16
|  | 104bb Data was already published in Mann et al. (2016b) | 2016 Mar 14
|  | 106cc Only a partial transit was observed | 2016 Mar 25
|  | 111 | 2016 Apr 21
|  | 122cc Only a partial transit was observed | 2016 Jun 20
|  | 186 | 2017 Jun 02
|  | 188cc Only a partial transit was observed | 2017 Jun 13
MEarth South | RG715 | 104bb Data was already published in Mann et al. (2016b) | 2016 Mar 14
|  | 111 | 2016 Apr 21
|  | 186 | 2017 Jun 02
|  | 188cc Only a partial transit was observed | 2017 Jun 13
|  | 191cc Only a partial transit was observed | 2017 Jun 29
Spitzer/IRAC | Channel 1 | 148 | 2016 Nov 08
|  | 150 | 2016 Nov 19
|  | 188 | 2017 Jun 13
|  | 190 | 2017 Jun 24
|  | 217 | 2017 Nov 18
Spitzer/IRAC | Channel 2 | 153 | 2016 Dec 06
|  | 186 | 2016 Jun 03
|  | 187 | 2017 Jun 08
|  | 191 | 2017 Jun 30
|  | 219 | 2017 Nov 29
HST /WFC3 | G141 | 184cc Only a partial transit was observed | 2017 May 23

##  3 Updated Stellar Parameters

Mann et al. (2016b) estimated effective temperature ($T_{\text{eff}}$),
stellar luminosity ($L_{*}$), stellar radii ($R_{*}$), and reddening ($A_{V}$)
for K2-33 using the combination of moderate resolution spectra, magnetic
models (Feiden, 2016), the transit-fit stellar density (assuming $e=0$), and
the distance to the Upper Scorpius OB association (de Zeeuw et al., 1999). The
more precise distance and photometry from Gaia Data Release Two (Gaia DR2;
Lindegren et al., 2018) enabled us to improve on these parameters.

We followed the methods in Mann et al. (2016b) for estimating bolometric flux
($F_{\mathrm{bol}}$ ) and $A_{V}$ by simultaneously comparing the flux-
calibrated spectra of K2-33 to unreddened young templates and observed
photometry of K2-33 from Gaia DR2 (Evans et al., 2018), the Two Micron All-Sky
Survey (2MASS; Skrutskie et al., 2006), the Wide-field Infrared Survey
Explorer (WISE; Cutri & et al., 2014), the Sloan Digital Sky Survey (SDSS;
Abolfathi et al., 2018), and the AAVSO All-Sky Photometric Survey (APASS;
Henden et al., 2012). This yielded $F_{\mathrm{bol}}$ of $2.40\pm 0.18\times
10^{-10}$erg cm-2 s-1 and $A_{V}=0.89\pm 0.13$. Adding the Gaia DR2 distance
yielded $L_{*}=0.146\pm 0.012L_{\odot}$, and using the $T_{\text{eff}}$ from
Mann et al. (2016b) and the Stefan-Boltzmann relation gave $R_{*}=1.017\pm
0.057R_{\odot}$. These are broadly consistent, but more precise than values
derived in Mann et al. (2016b).

To update $M_{*}$, we compared all above photometry to the magnetic model grid
from Feiden (2016) using a Gaussian prior on age of 10$\pm$3 Myr and $A_{V}$
from our spectral fit above. This yielded $M_{*}=0.571\pm 0.054M_{\odot}$.
Other model-based parameters ($R_{*}$, $L_{*}$ and $T_{\text{eff}}$) were
consistent with our spectral analysis, but we selected the empirical
parameters above. Using the $R_{*}$ above and the model-based mass gave a
stellar density ($\rho_{*}$) of $0.54\pm 0.12\rho_{\odot}$. A summary of our
adopted stellar parameters used in our transit analysis are listed in Table 2

Table 2: Updated Stellar Parameters for K2-33  Parameters | Value | Source
---|---|---
$A_{v}$ | 0.89 $\pm$ 0.13 | This Paper
$F_{\mathrm{bol}}$ ($\times 10^{-10}$ erg cm-2 s-1) | 2.40$\pm$ 0.18 | This Paper
$R_{*}$ $(R_{\odot})$ | 1.017$\pm$0.057 | This Paper
$M_{*}$ $(M_{\odot})$ | 0.571$\pm$0.054 | This Paper
$\rho_{*}$ $(\rho_{\odot})$ | 0.54$\pm$0.12 | This Paper
$L_{*}$ $(L_{\odot})$ | 0.146$\pm$0.012 | This Paper
$T_{\text{eff}}$ (K) | 3540$\pm$70 | Mann et al. (2016b)

We estimated limb-darkening parameters for all observations using the above
stellar parameters and the LDTK toolkit (Parviainen & Aigrain, 2015)777v1.5,
which added support for limb-darkening parameters past 5$\mu$m, which uses
PHOENIX models (Husser et al., 2013) and propagated uncertainties in stellar
parameters onto the output limb-darkening values. Different models give
different limb-darkening parameters at the level 0.04–0.08, so we inflated
these errors by those values to account for differences between input stellar
models. For MEarth, K2, and both Spitzer bands, we opted for a two-parameter
limb-darkening value ($g_{1}$, $g_{2}$). For HST, we used the non-linear
(4-parameter) limb-darkening formula detailed in Claret (2000) ($a_{1}$,
$a_{2}$, $a_{3}$, $a_{4}$). In each case, we used the relevant filter profile
from the instrument documentation. To simulate the HST spectral bands, we
multiplied the G141 filter profile with top-hat profiles corresponding to the
limits of each spectral band. The resulting limb-darkening estimates were fed
directly into MISTTBORN (K2, MEarth, and Spitzer) or Iraclis (HST), as
described in the next section. A summary of our adopted limb darkening
coefficients used in our transit analysis are listed in Table 5.

##  4 Transit Fitting

###  4.1 K2 and MEarth

We fit both the K2 and MEarth data individually using the MISTTBORN (MCMC
Interface for Synthesis of Transits, Tomography, Binaries, and Others of a
Relevant Nature) fitting code888https://github.com/captain-exoplanet/misttborn
described in Mann et al. (2016a) and in more detail in Johnson et al. (2018).
MISTTBORN uses BATMAN (Kreidberg, 2015) to generate model transit curves and
uses emcee (Foreman-Mackey et al., 2013) to explore the transit parameter with
a MCMC algorithm.

With MISTTBORN, we fit four parameters directly relating to the transiting
planet: time of periastron ($T_{0}$), orbital period of the planet ($P$),
planet-to-star radius ratio ($R_{p}/R_{\star}$), and impact parameter ($b$).
For the K2 data, all of these parameters evolved under uniform priors. For the
MEarth data, we placed a weak Gaussian prior on $T_{0}$ and $P$ ($\pm$0.01
days) to prevent the solution from wandering away from the data; other
parameters evolved under uniform priors. An additional parameter, stellar
density ($\rho_{\star}$), allowed us to replace the more common transit
duration and impose a Gaussian prior from our stellar parameters in Section 3.
For each of the two wavelengths, we fit two linear and quadratic limb-
darkening coefficients ($q_{1}$, $q_{2}$) following the triangular sampling
prescription of Kipping (2013). Both evolved under Gaussian priors derived
from LDTK (see Section 3). We locked the eccentricity ($e$) to zero because
gas drag and gravitational interactions are expected to dampen out
eccentricities and inclinations of extremely young planets like K2-33b (Tanaka
& Ward, 2004). In the case that eccentricity is non-zero, this does change the
overall depth, but it does not change the relative depths because all impact
parameters will move together.

An important difference between our analysis here and that of Mann et al.
(2016a) was that we fit the stellar variability simultaneously with the
transit parameters. This change was driven, in part, because of cases where
our removal of stellar variability would impact the transit, yielding a
systematically smaller transit depth than when fitting simultaneously. For
example, the transit depths determined by the two discovery papers differed
slightly despite using overlapping data ( 0.23 vs. 0.19%; David et al., 2016b;
Mann et al., 2016b).

To this end, we used the celerite package (Foreman-Mackey et al., 2017) for
Gaussian Process fitting which has now been incorporated into MISTTBORN. We
first tested the “SHOM” kernel, i.e., a mixture of two stochastically driven
damped simple harmonic oscillators (SHO) with periods $P$ and $0.5P$. However,
we found the parameters of the second SHO (at half the period) was poorly
constrained, which is likely because the primary rotation signal dominates
(Nicholson & Aigrain, 2022). Instead, we adopted a single SHO fit, which added
three free parameters: the log of period ($\ln(P_{GP})$), the log of the
variability amplitude ($\ln{A}$), and the decay timescale ($\ln{Q}$).

When fitting the K2 data, we evolved all GP parameters under uniform priors.
For the MEarth data, there is insufficient out-of-transit data to constrain
the GP, so we applied a Gaussian prior on $\ln(P_{GP})$ and $\ln{A}$ based on
the output from the K2 fit. Since the two datasets are at slightly different
wavelengths (0.6$\mu$m vs 0.8$\mu$m) the spot-variability amplitude might not
be the same (and hence so might be the GP amplitude). However, factor-of-two
changes to the amplitude prior in either direction did not change the
resulting transit depths, most likely because the GP is able to adjust to each
transit between the long data gaps between MEarth observations.

For both datasets, we ran MISTTBORN with 100 walkers for 250,000 steps after a
burn-in of 50,000 steps. A comparison to the autocorrelation time suggested
this was more than sufficient for convergence (Goodman & Weare, 2010). We
summarize the MISTTBORN output in Table 3. The fit to the K2 data is shown in
Figure 2.

Table 3: Parameters of K2-33 for Different Data Sets Description | Parameter |  K2 | MEarth |  Spitzer [3.6] |  Spitzer [4.5]
---|---|---|---|---|---
First Transit Mid-point |  $T_{0}$ (KBJD999KBJD = BJD-2454833) | $2065.6921^{+0.0028}_{-0.0027}$ | $2065.6923\pm 0.0045$ | $2065.6982^{+0.0046}_{-0.0048}$ | $2065.6954^{+0.0079}_{-0.0081}$
Orbital Period |  $P$ (days) | $5.42503\pm 0.00034$ | $5.424832^{+2.6\times 10^{-5}}_{-2.5\times 10^{-5}}$ | $5.424851^{+2.6\times 10^{-5}}_{-2.7\times 10^{-5}}$ | $5.424862^{+4.2\times 10^{-5}}_{-4.1\times 10^{-5}}$
Planet-Star Radius Ratio | $R_{P}/R_{\star}$ | $0.04735^{+0.00099}_{-0.00096}$ | $0.0489^{+0.0019}_{-0.0020}$ | $0.03545\pm 0.00086$ | $0.03746^{+0.00083}_{-0.00079}$
Impact Parameter | $b$ | $0.19^{+0.18}_{-0.13}$ | $0.27^{+0.2}_{-0.18}$ | $0.24^{+0.21}_{-0.17}$ | $0.22^{+0.2}_{-0.15}$
Stellar Density |  $\rho_{\star}$ ($\rho_{\odot}$) | $0.464^{+0.042}_{-0.066}$ | $0.527^{+0.062}_{-0.098}$ | $0.483^{+0.04}_{-0.1}$ | $0.476^{+0.038}_{-0.088}$
Limb-darkening coefficient | $q_{1}$ | $0.48^{+0.083}_{-0.082}$ | $0.401^{+0.079}_{-0.078}$ | $0.072^{+0.033}_{-0.031}$ | $0.06^{+0.027}_{-0.024}$
Limb-darkening coefficient | $q_{2}$ | $0.351\pm 0.046$ | $0.286\pm 0.048$ | $0.209^{+0.033}_{-0.034}$ | $0.204\pm 0.034$
Log Period | $\ln{P}$ | $1.838\pm 0.015$ | $1.836\pm 0.015$ | $2.6^{+2.57}_{-0.83}$ | $0.99^{+0.68}_{-0.4}$
Log Variability Amplitude | $\ln{A}$ | $-8.96^{+0.26}_{-0.22}$ | $-10.14^{+0.44}_{-0.39}$ | $-13.84^{+1.04}_{-0.86}$ | $-13.93^{+1.12}_{-0.92}$
Log Decay Timescale | $\ln{Q}$ | $0.926^{+0.13}_{-0.087}$ | $0.048^{+0.083}_{-0.036}$ | $6.2^{+2.9}_{-3.0}$ | $4.5^{+3.3}_{-2.9}$

###  4.2 Spitzer

Spitzer’s large intra-pixel sensitivity and pointing jitter can cause the
measured flux of a source to vary up to 8%, depending on where it falls on a
pixel (Ingalls et al., 2012). However, years of high-precision observations
with Spitzer have provided a number of methods to correct for model variations
in the photometric response (Ingalls et al., 2016). We tested multiple
methods, which includes using a high-resolution pixel variation gain map
(PMAP; Ingalls et al., 2012) and nearest neighbors (NNBR; Lewis et al., 2013),
but found the most consistent results using the BiLinearly-Interpolated
Subpixel Sensitivity (BLISS; Stevenson et al., 2012) mapping technique, which
we briefly summarize below. We note that changing to NNBR did not change any
of the conclusions, only the overall fit quality.

The BLISS Mapping Technique is provided as an optional part of the POET
pipeline. It corrects for both the position-dependent (intra-pixel effect) and
time-dependent (ramp effect) Spitzer systematics. The parameters in the BLISS
mapping include the grid size of the sub-pixel sensitivity map, the
astrophysical light curve model (e.g., transit, secondary eclipse), and the
ramp model.

For the Channel 1 data, to avoid overfitting, we followed POET’s
recommendation of choosing a grid size in which a nearest-neighbor
interpolation would not outperform a BLISS interpolation (Stevenson et al.,
2012). We tested various grid sizes and selected 0.007 pixels. For the Channel
2 data, we used the most recent fixed intra-pixel sensitivity map from May &
Stevenson (2020).

The transit was modeled using the Mandel & Agol (2002) transit model and
several different ramp parameterizations were tested: linear, quadratic,
rising exponential, falling exponential, quartic-log + quadratic polynomial,
log + linear ramp, and a no-ramp model. For further information about each
individual ramp, see Stevenson et al. (2012). For each light curve, we
determine the best ramp model based on three metrics: 1) the difference in
predicted t0 value compared to the expected t0, 2) the difference in the
predicted transit duration compared to the expected transit duration, and 3)
the overall minimal red noise levels in the fit residual, assessed by
considering the root-mean-squared (RMS) binned residuals as a function of
different bin sizes with the theoretical uncorrelated white noise. A small
$T_{0}$ difference, small duration difference, and low RMS were all favored.

The time-dependent component of the transit model consisted of the mid-transit
time ($T_{0}$), planet-to-star radius ratio ($R_{p}/R_{s}$), orbital
inclination ($\cos{i}$), semi-major axis ratio ($a/R_{*}$), and system flux
($\mu Jy$), as well as, parameters associated with the ramps (e.g., ramp
phase, ramp amplitude, ramp constant offset, etc.). These parameters were
explored with an MCMC process, using 4 walkers with 200,000 steps and a burn-
in region of 30,000 steps. The period was fixed to 5.4248 days, based on our
analysis of the K2, MEarth, and an initial reduction of the Spitzer data
(which provided a sub-second precision period). As recommended by May &
Stevenson (2020), temporal binning was not used for Channel 2. The resulting
POET fits for each light curve are summarized in Table 4.

Table 4: Best Fit Parameters for Spitzer Data Using POET AOR | Transit Num | BIC Value | Best Ramp | T0 [BJD] | cosi | a/R∗ | System Flux [$\mu Jy$] | Rp/R∗
---|---|---|---|---|---|---|---|---
3.6 $\mu$m
60655360 | 148 | 10773 | quadratic | 2457701.57 | 0.07 | 8.46 | 31325 | 0.037 $\pm$ 0.004
60658432 | 150 | 10765 | quadratic | 2457712.42 | 0.07 | 8.79 | 31172 | 0.063 $\pm$ 0.002
60661504 | 188 | 14433 | quartic-log + quartic polynomial | 2457918.57 | 0.00 | 10.52 | 30291 | 0.035 $\pm$ 0.003
60664576 | 190 | 14396 | linear | 2457929.41 | 0.03 | 10.00 | 29864 | 0.036 $\pm$ 0.002
60664576 | 217 | 14407 | quartic | 2458075.89 | 0.07 | 8.63 | 30618 | 0.035 $\pm$ 0.002
4.5 $\mu$m
60656128 | 153 | 10754 | quadratic | 2457728.70 | 0.00 | 10.50 | 19369 | 0.050 $\pm$ 0.003
60659200 | 186 | 14406 | log | 2457907.72 | 0.07 | 8.61 | 19550 | 0.037 $\pm$ 0.004
60662272 | 187 | 14390 | quadratic | 2457913.15 | 0.00 | 10.27 | 19662 | 0.038 $\pm$ 0.003
60665344 | 191 | 14390 | rising exponential | 2457934.85 | 0.07 | 8.54 | 19677 | 0.037 $\pm$ 0.002
60668416 | 219 | 14388 | rising exponential | 2458086.74 | 0.07 | 8.80 | 19664 | 0.041 $\pm$ 0.001
Figure 4: Planet-to-star radius ratio ($R_{p}/R_{*}$) posteriors from the POET
MCMC fits for Spitzer Channel 1 (left) and Channel 2 (right) with a bin size
of 30. Each color represents a different transit (AOR), where gray is used to
highlight the outlier transit in each channel (AOR 60658432 in Channel 1 had a
flare and AOR 60656128 in Channel 2 had poor pointing). The green dashed
vertical line is the planet-to-star radius ratio for K2, which does not
overlap with any of the Spitzer datapoints, except for the outlier transit in
Channel 2.  Figure 5: Top:  Stacked phase-folded light curve of K2-33 observed
in K2 (pink), MEarth (blue), Spitzer/Channel 1 (purple), Spitzer/Channel 2
(orange), and HST (green). The HST data here is the white light curve. Data
points correspond to the light curve binned in phase using a bin size of 10,
200, 1500, and 1500 for K2, MEarth, Channel 1, and Channel 2, respectively.
Solid line corresponds to the best-fit (highest likelihood) model from our
MCMC fit. Typical error bars are derived from scatter in the out-of-transit
data points. The transit depth for the optical wavelengths is $\sim$ 2 times
deeper than the transit depth for the near-infrared wavelengths. Bottom:
Residuals using the binned points.

A drawback of the POET results is that we could not fit all transits
simultaneously. This is particularly important since several of the Spitzer
transits yield discrepant depths (Figure 4). To this end, we re-fit the output
light curves from POET with MISTTBORN, as described in Section 4.1. We fit the
two channels separately, including all but one transit in each band. The
excluded observations were AOR 60658432 from Channel 1 (which had a strong
flare) and AOR 60656128 from Channel 2 (which had extremely poor pointing).
These two transits also yielded outlier transit depths from POET (see Figure
4).

Unlike with POET, MISTTBORN cannot simultaneously fit for instrumental effects
with the transit, so we attempt to account for that by including the GP. We
tried several GP kernels, including Matern-3/2 and a simple white noise term,
but found the SHO we used for stellar variability worked as well or better.
Since the periodic signal in this case is instrumental, not stellar, we allow
all GP parameters to evolve under uniform priors (instead of forcing it to
match the rotation period, for example). As above, the only transit parameters
we place priors on were limb-darkening and stellar density.

The resulting MISTTBORN fits are summarized in Table 3. Other than depth, the
transit-specific parameters were in broad agreement with those from K2 and
MEarth. Excluding the two problematic transits, the transit depths from our
MISTTBORN fits are in excellent agreement with those from POET. The MISTTBORN
fits gave similar or smaller uncertainties than individual POET fits, but
still larger than expected from a simple weighted mean of the individual POET
fits. We consider the MISTTBORN fits to be more realistic.

Importantly, both the MISTTBORN and POET fits to the data yielded a transit
depth $\simeq$half that from the optical data. This difference was not
explainable through the fitting, as it is clear even in the extracted light
curves (Figure 5). Further, as we will show in Section 4.3, the effect is seen
in the HST data as well.

###  4.3 HST

####  4.3.1 White Light Curve Fitting

In order to fit the extracted light curves, we must first correct the time-
dependent systematics (ramps) introduced by the WFC: 1) a long term ramp that
occurs throughout the visit and typically has a linear behavior and 2) a short
term ramp that occurs during each orbit and typically has an exponential
behavior (Tsiaras et al., 2016).

The light curve is fitted using a transit model from PyLightcurve
101010https://github.com/ucl-exoplanets/pylightcurve, and with an instrumental
systematics function $R(t)$ (Tsiaras et al., 2016) :

$R_{w}(t)=n_{w}(1-r_{a1}(t-T_{0}))(1-r_{b1}e^{-r_{b2}(t-t{{}_{0}})})$ (1)

where $n_{w}$ is a normalization factor, $t$ is time of the data, $T_{0}$ is
the model mid-transit time, $t_{0}$ is the starting time of each HST orbit,
$r_{a1}$ is the slope of the linear systematic trend, $r_{b1}$ and $r_{b2}$
are the coefficients of the exponential systematic trend. Since the first
orbit ramp (Orbit 2) is steeper compared to the ramps of the subsequent orbits
(Figure 6), a different set of the short-term exponential coefficients
($for_{b1}$ and $for_{b2}$) were used to fit this ramp. Removing these
coefficients will increase the root-mean-square (RMS) of the residuals.

Due to only having half a transit and missing both the egress and out-of-
transit baseline, we limit our ability to constrain the orbital parameters.
Therefore, we locked the inclination($i$), the semi-major axis-to-star radius
ratio ($a/R_{*}$), the period ($P$), the argument of periastron ($\omega$),
and the eccentricity ($e$) to the values from a combined fit of data from K2,
MEarth, and Spitzer. Limb-darkening coefficients ($a$) were fixed to the
model-derived values (Table 5).

The only free parameters were coefficients for Hubble systematics, the
normalization factor ($n_{w}$), the planet-to-star radius ratio
($R_{p}/R_{*}$), and the mid-transit ($T_{0}$). These parameters were explored
with an MCMC process, using 200 walkers with 40,000,000 steps, and a burn-in
region of 6,000,000 steps. This was sufficient for convergence based on the
autocorrelation time. The final $R_{p}/R_{*}$ result for the broadband (white)
light curve can be found in Table 5.

The fit white light curve is shown in Figure 6. As can be seen in the bottom
panel, the residuals are within 300 ppm for both Orbit 2 and Orbit 3, but not
for Orbit 4, due to the remaining uncorrected systematics.

The autocorrelation of the residuals was 0.37. This value is slightly above
the threshold used by Tsiaras et al. (2018), who defined a fit to be
successful if the autocorrelation was $<$ 0.3. In that work, they applied
Iraclis to 30 gaseous planets taken by HST/WFC3. However, their study was
focused on full transits and targets significantly less variable than K2-33.
Our standard deviation of the residuals relative to the photon noise
($\bar{\sigma}$) is 1.36 and this is within the reasonable limits from Tsiaras
et al. (2018) (1.10 $<\bar{\sigma}<$ 2.75). We considered this fit to be
reasonable, but given the partial transit loss, more data would be required to
further investigate the fit quality.

Figure 6: White light curve of K2-33 observed with HST/WFC3. The top panel
shows the raw normalized data (green circles) overplotted with the systematic
model (gray line). The middle panel shows the best fit model (gray line) from
our MCMC fit overplotted with the systematic-corrected data (green circles).
The bottom panel shows residuals from the light curve.

####  4.3.2 Spectral Light Curve Fitting

To extract the planetary spectrum from the spectral light curves, we must
first correct for the systematics. Iraclis fits the spectra light curve using
the divide white method, first introduced by Kreidberg et al. (2014), in which
a transit mode is multiplied with an instrumental systematics function
$R_{\lambda}(t)$:

$R_{\lambda}(t)=n_{\lambda}[1-r_{a,\lambda}(t-T_{0})]\frac{LC_{W}}{M_{W}},$
(2)

where $n_{\lambda}$ is the normalization factor for each wavelength bin,
$r_{a1}$ is the linear long term ramp, $LC_{W}$ is the white light curve, and
$M_{W}$ is the model of the white light curve.

The transit model was fixed to the same orbital parameters as the white light
curve. We did not fit for the limb-darkening coefficients. The only free
parameters are the linear long term ramp, the normalization factor, and the
planet-to-star radius ratio ($R_{p}/R_{*}$). The parameter space was explored
with an MCMC process using 100 iteration walkers, 10,000,000 steps and a burn-
in region of 2,000,000 steps. The final $R_{p}/R_{*}$ results for each
spectral bin and their uncertainties can be found in Table 5. All the spectral
light curves for each wavelength bins are plotted in Figure 7. The standard
deviation of the residuals is on average 1.12 times the expected photon noise
($\bar{\sigma}$), and the residuals autocorrelation ($R^{2}$) is 0.18. Both
these values fall within the range of the spectral light curve fits from
Tsiaras et al. (2018) ($\bar{\sigma}<1.17$ and $R^{2}<0.20$).

Figure 7: Left: Spectral light curve annotated with the mean wavelength for
each bins. Right: Residuals of the fit, along with the standard deviation
relative to the photon noise ($\bar{\sigma}$), the reduced chi-square
($\bar{\chi}^{2}$), and the square of the correlation ($R^{2}$)  Table 5:
Limb-darkening Coefficients For Each Data Set

Telescope |  $\lambda$ range ($\micron$) | $g_{1}$ | $g_{2}$ |  |  |
---|---|---|---|---|---|---
K2 | 0.42–0.90 | $0.479\pm 0.050$ | $0.211\pm 0.036$ |  |  |
MEarth | 0.69–1.00 | $0.317\pm 0.048$ | $0.274\pm 0.035$ |  |  |
Spitzer Channel 1 | 3.13–3.96 | $0.112\pm 0.015$ | $0.155\pm 0.019$ |  |  |
Spitzer Channel 2 | 3.92–5.06 | $0.085\pm 0.011$ | $0.124\pm 0.018$ |  |  |
Telescope |  $\lambda$ range ($\micron$) | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $R_{p}/R_{*}$
HST | 1.0880–1.1303 | -0.133292 | 1.848572 | -1.917848 | 0.665771 | $0.043544^{+0.004269}_{-0.004879}$
| 1.1303–1.1726 | -0.121481 | 1.783907 | -1.858905 | 0.646929 | $0.040560^{+0.004145}_{-0.004737}$
| 1.1726–1.2149 | -0.115123 | 1.742967 | -1.820934 | 0.634099 | $0.035911^{+0.004172}_{-0.004768}$
| 1.2149–1.2571 | -0.107202 | 1.694365 | -1.767156 | 0.614537 | $0.036135^{+0.00303}_{-0.00303}$
| 1.2571–1.2994 | -0.100966 | 1.657434 | -1.729879 | 0.601442 | $0.032079^{+0.003439}_{-0.003439}$
| 1.2994–1.3417 | -0.106114 | 1.672739 | -1.746672 | 0.606971 | $0.034728^{+0.003642}_{-0.003642}$
| 1.3417–1.3840 | -0.137048 | 1.853932 | -1.936428 | 0.671591 | $0.037669^{+0.002757}_{-0.00315}$
| 1.3840–1.4263 | -0.150512 | 1.907895 | -1.981589 | 0.684045 | $0.037647^{+0.003729}_{-0.003262}$
| 1.4263–1.4686 | -0.158651 | 1.977145 | -2.061429 | 0.710897 | $0.035890^{+0.003073}_{-0.003512}$
| 1.4686–1.5109 | -0.154233 | 2.000803 | -2.104788 | 0.727856 | $0.029557^{+0.003642}_{-0.003187}$
| 1.5109–1.5531 | -0.149891 | 2.013220 | -2.127071 | 0.737913 | $0.040000^{+0.00318}_{-0.00318}$
| 1.5531–1.5954 | -0.123405 | 1.967341 | -2.120444 | 0.742160 | $0.035862^{+0.003098}_{-0.002711}$
| 1.5954–1.6377 | -0.092558 | 1.879567 | -2.055896 | 0.726144 | $0.039918^{+0.00286}_{-0.00286}$
| 1.6377–1.6800 | -0.057187 | 1.743134 | -1.937097 | 0.689630 | $0.035785^{+0.005314}_{-0.005314}$
HST White Light Curve | 1.0880–1.6800 | -0.127289 | 1.860330 | -1.965358 | 0.685232 | $0.036702^{+0.001698}_{-0.001485}$

Note. — Limb darkening values for all dataset were calculated using the LDTK
toolkit (Parviainen & Aigrain, 2015). For the K2 , MEarth , and Spitzer data
set, limb darkening priors are provided as the traditional linear and
quadratic terms, but were fit using triangular sampling terms.

Overall, the HST data was too imprecise to confidently detect features,
particularly given the loss of egress. However, even the partial transits
confirm the large difference between the optical and NIR transit depths; HST
and Spitzer both yielded transit depths almost half as deep as the depths from
K2 and MEarth. This rules out that the depth difference is due to
unconstrained systematics in the datasets. We explored two possible reasons
for this difference: spots on the stellar surface (Section 5) and hazes in the
planet’s atmosphere (Section 6).

##  5 The effect of surface inhomogeneities on the transmission spectrum

Unocculted spots on the stellar surface can make a transit appear deeper. If
the transit chord has fewer spots, the planet will block a statistically
warmer (brighter) region of the star. The effect is reversed for unocculted
plages. Since surface inhomogeneities vary in intensity with wavelength, the
effect is wavelength-dependent and can mimic or complicate transmission
spectroscopy of the planet (e.g., Kreidberg et al., 2015; Rackham et al.,
2017). Similarly, occulted spots or plages can make the transit shallower or
deeper, respectively. We show below that the large difference between the
optical and NIR transit depths could be explained if more than half ($>$60%)
of the stellar surface is covered in cool spots and the transit chord is
clear, or if the transit cord is heavily populated by plages while the rest of
the star is clear.

###  5.1 Planet crossing over spots and plages

If the planet were crossing a significant number of spots or plages during
transit, we would expect to see two effects. First, we would see distortions
in the light curve during transit (e.g., as seen for HAT-P-11 in Morris et
al., 2017). While this would be a challenge to see with long-cadence (cadence
= 1765.5 s) data from K2, it should be readily visible in the high-cadence
(cadence = 60 s) MEarth data. The second effect should be variations in the
transit depth between transits due to differences between the stellar and
orbital periods; each transit crosses a slightly different region of the star,
changing the properties of the occulted surface. Changes in the underlying
spot pattern are also visible in the the K2 light curve (Figure 2), likely due
to differential rotation (as noted in David et al., 2016a). However,
individual transit depths are consistent with each other within measurement
errors, the scatter in the light curve during transit is the same as outside
of the transit, and there are no clear morphological changes in the MEarth or
K2 data during transit. Planet-crossing plages would also need to be long-
lived to explain both the K2 and MEarth data, which spans 2.7 years and has
more than 150 full stellar rotations.

Due to the reasons above, we assumed that the transit chord is pristine for
the rest of our analysis and focused on the effects of unocculted spots on the
rest of the star. We note that consistency between transit depths with time
disfavors significant unocculted spots, since it is unlikely that any given
transit would have the same spot coverage fraction. However, this is easier to
get around with polar or symmetric spot distributions and small variations in
the spot coverage fraction may get lost in the noise.

###  5.2 Unocculted spots

####  5.2.1 What kind of spot coverage can reproduce the observed transit
depths?

Faculae have shown to have less of an effect in the stellar variability in the
Kepler light curves, justifying our reasons for not including them in our
analysis (Johnson et al., 2021). To estimate what kind of spots are required
to reproduce the observed transit depths, we followed a modified procedure
from Rackham et al. (2017), Thao et al. (2020), and Libby-Roberts et al.
(2021). To summarize, we assumed the surface can be described by a simple two-
temperature model of three parameters: the surface temperature
($T_{\rm{surf}}$), the spot temperature ($T_{\rm{spot}}$), and the fraction of
the star with $T_{\text{eff}}$$=$$T_{\rm{spot}}$ ($f_{S}$).

As explained by Libby-Roberts et al. (2021), $f_{S}$ is actually time variable
for these two reasons: 1) it represents the surface seen from Earth, therefore
$f_{S}$ changes as the star rotates, and 2) the overall surface spot pattern
can change with time. We ignore this effect for now. As we will show in
Section 5.2.2, the overall spot coverage needs to be relatively stable over
the rotational period, or the light curve would show more significant
variability in the K2 data.

We built our synthetic star using surface brightness estimates from
PHOENIX/BT-SETTL 111111https://phoenix.ens-lyon.fr/Grids/BT-
Settl/CIFIST2011_2015/ atmosphere models (Allard et al., 2013, 2012). For both
the spots and surface, we assumed Solar abundances and adopted simple linear
interpolation between grid points in $T_{\text{eff}}$. Surface gravity
($\log(g)$) has negligible effects on the spectrum at the resolution
considered here, so we fixed it to 4.5. We adopted a single transit depth
($D$), which we would measure for an unspotted star (i.e., a featureless
spectrum), although this also can work as a normalization constant for a model
spectrum (as we do in Section 6). We blocked light only from the region
described by $T_{\rm{surf}}$, which we then combined with the spotted region
to compute an observed transit depth after normalizing to the combined
spectrum with no light blocked (the out-of-transit spectrum). This yielded an
observed transit depth at any wavelength, which we convolved with the relevant
filter profile (K2, MEarth, HST, and Spitzer) to compare to the observed data.
We eliminated one free parameter by fixing $T_{\rm{surf}}$ to the temperature
assigned from the optical spectrum (3540 K; Section 3). This assumption was
imperfect, as it is possible to reproduce the observed spectrum with
significant spot coverage and/or cooler spots while increasing
$T_{\rm{surf}}$. However, allowing $T_{\rm{surf}}$ to go to higher values
while forcing it to reproduce the observed spectrum did not significantly
change the answer, as we explore further in Section 5.2.3.

We compared our model to the observed transit depths data within an MCMC
framework using emcee (Foreman-Mackey et al., 2013). Each of the three free
parameters evolved under uniform priors, with physical limits ($0<f_{S}<1$ and
$0<D<1$), and those imposed by our model grid (1500 K$<$$T_{\rm{spot}}$$<$4000
K). We ran the chain with 12 walkers for 100,000 steps, following a burn-in of
10,000 steps. This was sufficient for convergence based on the autocorrelation
time.

Figure 8: The best-fit model of a planet with a flat transmission spectrum
crossing a heavily spotted star along a pristine transit chord (blue). The
observed transit depths are plotted in red circles and the labeled bars at the
bottom indicate which instruments took those observations. The yellow lines
show 100 random samples in $T_{\rm{spot}}$ $f_{S}$, and $T_{\rm{surf}}$ taken
from our fit to the IGRINS spectra (see Section 5.2.3). None of the results
drawn from the spectroscopic constraints are consistent with the data; the
only samples that come closest to the optical data fail to reproduce the
Spitzer data. Results are similar using the posterior from our fit to the
SNIFS+SpeX data.  Figure 9: Corner plot of $T_{\rm{spot}}$ and $f_{S}$ from
fitting a flat planetary transmission spectrum to the observed transit depths
and assuming a pristine transit chord (blue) compared to the same parameters
constrained by fitting the stellar spectrum from IGRINS (black) or SNIFS+SpeX
(red) spectra with a two-temperature model. Some solutions with low $f_{S}$
and low $T_{\rm{spot}}$ values were trimmed for clarity.

We show the best-fit transmission spectrum in Figure 8, and a subset of the
posteriors in Figure 9. The results suggest we can explain the deeper optical
transits if $71\%^{+14}_{-6}$ of the star is covered in spots with
$T_{\rm{spot}}$=$2750^{+200}_{-250}$ K. Such large spot coverage fractions are
rare, even for young stars (Morris, 2020), but they have been observed on
stars in star-forming regions (e.g., Gully-Santiago et al., 2017) so we could
not dismiss this option on statistical grounds alone. We explore constraints
provided by the transit light curve and stellar spectrum in the next two
sections.

The star is likely to harbor spots with a range of temperatures. However,
adding an additional spot with a different temperature and coverage fraction
did not alter our conclusions: the preferred solution is 1) to have both spots
of similar temperatures, or 2) to have one of the spots be equivalent to the
result above, and the other being either too cold or possessing too small a
coverage fraction to impact the final result. In either case, the combined
effect is not significantly different from the single spot temperature model.
The conclusion is still that the spot coverage fraction must be $>60$% to
explain the observed transit depths.

####  5.2.2 Spots from the light curve

Surface inhomogeneities imprint on the stellar light curve as the star
rotates, with spots moving into darker regions of the stellar disk before
disappearing from view. All else being equal, larger/more spots and plages
will lead to higher-amplitude variation. As discussed in Luger et al. (2021b),
it is difficult to use this relation in reverse to infer spot coverage
fractions from light curves because many different spot coverage fractions and
distributions can produce similar light curves. Rather than trying to use the
light curve to constrain the spot coverage fraction, we instead attempted to
explore what kind of surface patterns could reproduce both the light curve and
the observed transit depths.

To this end, we used Fleck121212https://github.com/bmorris3/fleck (Morris,
2020), which produced light curves from simple limb-darkened spotted stars. We
compared the model-generated light curves to the light curve within a MCMC
framework using emcee. We restricted our analysis to the K2 curve; the HST and
Spitzer curves contain too little out-of-transit data, and the MEarth out-of-
transit monitoring is not precise enough to add significant constraints above
the K2 data. The K2 data likely has long-term systematics, which we remove by
fitting the light curve with a low-order polynomial before comparing it to the
Fleck-generated model.

We fixed the rotation period to that from Mann et al. (2016b) and restricted
the spot contrast to be consistent with the results from Section 5.2.1.
Specifically, we convolved the K2 filter profile with BT-SETTL models
(described in Section 5.2.1) using the assumed $T_{\rm{surf}}$ and posterior
$T_{\rm{spot}}$ values. This yielded a spot contrast for every value of
$f_{S}$ (also from the posterior), which we included in our likelihood. The
result is that the contrast was fixed based on the model $f_{S}$ value derived
from the spot radii. This also effectively limited $f_{S}$, as values
$\lesssim 0.6$ had no corresponding contrast that could reproduce the
transmission spectrum. The goal was to determine if the high spot coverage
fraction implied by the transmission spectrum could reproduce the light curve
(and what would the resulting star look like).

There were three fit parameters for each spot: latitude, longitude, and
radius. Each parameter evolved under uniform priors with physical limits
(e.g., latitude from -90 to 90 degrees). We tested using three to seven spots;
the results did not change in any way relevant to our analysis.

We highlight some representative results in Figure 10. In order to
simultaneously reproduce the low variability ($<2\%$) in the light curve and
the inferred large spot coverage, highly symmetric polar spots that cover a
large portion of the star are required. Previous studies have suggested the
presence of polar spots in M-dwarfs (Doyle et al., 2018; Roettenbacher & Vida,
2018) so we can not rule out this possibility. Most often the model prefers
either a slight asymmetry in the polar spots and/or 1-3 smaller equatorial
spots to explain the perturbations offset from the major sinusoidal pattern
seen in the light curve (likely spot clusters at different longitudes). This
is consistent with earlier findings that large spot coverage fractions can
still produce relatively low variability Johnson et al. (2021).

Figure 10: Left:  Example results from Fleck, showing a portion of the light
curve (black points) and best-fit models (purple curves). The model only
captures the mean variability, failing to capture the peaks and troughs
because it does not include long-term spot evolution. Middle: Posteriors of
$f_{S}$. Right: Three example stellar surfaces randomly pulled from our fit
posteriors. The colors correspond to 2 (blue), 3 (red), and 4 (grey) spots. To
explain the low variability in the light curve and the inferred spot coverage,
the spots must cover a large portion of the star and be located at the poles.

As expected, it was possible to recreate a star that reproduced the light
curve and could explain the transit depths using a featureless planetary
spectrum; a single-band light curve was not uniquely constraining. However,
the resulting stellar surface had to be finely tuned to contain mostly polar
spots. In Section 5.2.2, we assumed that $f_{S}$ was not a time variable;
while this assumption was not totally correct, the light curve results
suggested that it is reasonable nonetheless. The only way to reproduce the low
photometric variability was with mostly polar and/or symmetric spots (Figure
10), which, by definition, produce small temporal changes in $f_{S}$. Or,
stated another way, if $f_{S}$ were changing at $\gtrsim 5\%$, this would be
evident in the light curve, validating our simplification in Section 5.2.2.

####  5.2.3 Spot properties from the observed stellar spectrum

Large cool spots change both the spectral shape and the strength of molecular
features, as we show in Figure 11. In extreme cases, cool spots can generate
molecular features that should otherwise not be present for a homogeneous
surface. Numerous earlier studies have taken advantage of this to study spot
properties of young stars (e.g., Petrov et al., 1994; Fang et al., 2016),
likely providing more accurate pictures of stellar surfaces than light curve
based metrics alone (e.g., Gully-Santiago et al., 2017).

Figure 11: Illustration of the effect of spots on the observed stellar
spectrum using BT-SETTL atmosphere models for a star with
$T_{\rm{surf}}$$=3500$ K and $T_{\rm{spot}}$$=2800$ K and spot coverage
fractions of 80% (green), 50% (blue), and 20% (green). The left panel shows a
moderate-resolution optical–NIR spectrum (matching the wavelength coverage and
resolution as our observational data). Since spectra only generally provide
relative flux measurements, each model is normalized to yield the same total
luminosity. The change over the three models is small in the optical, and
could largely be fit out by assuming a cooler surface temperature or higher
reddening. However, as large spots have a much larger impact on the NIR
spectrum, the combined dataset can provide strong constraints on both
$T_{\rm{surf}}$ and $T_{\rm{spot}}$.

Since the effect of spots on the observed spectrum is strongly wavelength
dependent (e.g., the cooler region has a smaller overall contribution in the
optical), it can be detected in spectra with broad wavelength coverage. The
optical data mostly probes $T_{\rm{surf}}$ and reddening, while the NIR data
is more sensitive to changes in $T_{\rm{spot}}$ and $f_{S}$. We used the
moderate resolution ($R\simeq 1000-2000$) spectra from Mann et al. (2016b) to
explore the spot properties for K2-33. The spectra were built from two
overlapping datasets, one from the SNIFS spectrograph (0.32–0.95$\mu$m; Lantz
et al., 2004) and one from the SpeX spectrograph (0.7–2.4$\mu$m; Rayner et
al., 2003). Both spectrographs provide excellent relative flux calibration
(2-4%; Rayner et al., 2009; Mann et al., 2013) ideal for this work. We refer
readers to the discovery paper (Mann et al., 2016b) for more details on the
observations and reduction of these data.

We fit the spectra with a two-temperature model (e.g., Gully-Santiago et al.,
2017; Morris et al., 2019) derived from BT-SETTL atmosphere models. In total,
we had 13 free parameters for the fit. The first three, $T_{\rm{surf}}$,
$T_{\rm{spot}}$, and $f_{S}$ were as described in Section 4, but without any
limit on $T_{\rm{surf}}$ and an additional requirement that $T_{\rm{spot}}$ be
at least 50 K less than $T_{\rm{surf}}$. Six additional parameters described
the normalization and slope of the blue-optical, red-optical, and NIR data
($a1$, $a2$, $b1$, $b2$, $c2$, and $c2$). These were needed because the SNIFS
blue end, SNIFS red end, and SpeX data were calibrated separately (SNIFS arms
are separated with a dichroic at 0.54$\mu$m); as a result, the three regions
are often offset at the 2-5% level even after calibration (Mann & von Braun,
2015). Two additional parameters handle small wavelength offsets between the
data and models, one for the optical/SNIFS data ($\lambda_{\rm{off,opt}}$) and
one for the NIR/SpeX data ($\lambda_{\rm{off,NIR}}$). The remaining two
parameters were extinction ($A_{V}$) and a parameter to handle underestimated
uncertainties or other systematics in the spectra ($\sigma_{f}$), which we
treated as a fractional error added in quadrature with measurement errors.

Following Mann et al. (2013), we masked specific regions where the models
poorly reproduced observed spectra of stars with empirically-determined
temperatures (from long-baseline interferometry). No other direct constraints
were put on the data.

All parameters were fit using emcee. We placed Gaussian priors on the slope
parameters ($a2$, $b2$, $c2$) based on estimates of the flux calibration from
Mann et al. (2013) and Rayner et al. (2009). All other parameters evolved
under uniform priors. We ran the MCMC chain with 100 walkers for 500,000
steps, which was more than sufficient for convergence based on the
autocorrelation time.

We added the resulting $f_{S}$ and $T_{\rm{spot}}$ distributions in Figure 9.
As expected, spot temperatures were unconstrained when the spot fraction was
low, and the spot fraction was unconstrained when $T_{\rm{spot}}$ was close to
$T_{\rm{surf}}$. The general result was that either spots must be relatively
warm $\gtrsim 3000$ K or cover a small fraction ($<10$%) of the star.

We performed a similar fit using the high-resolution IGRINS spectra from Mann
et al. (2016b). IGRINS covers most of the $H$\- and $K$-bands with a
resolution of $\simeq$45,000. Nearly-identical data were used effectively in
Gully-Santiago et al. (2017) to study the spot properties of a young star. Of
the 54 IGRINS orders, we focused on the four with the highest SNR, lowest
telluric contamination, and highest sensitivity to spots (e.g., right panel of
Figure 11). The fit parameters were similar to what we used for the moderate-
resolution data, but each order had separate flux calibration terms, a common
wavelength offset, and we did not fit for $A_{V}$ (which has a minor impact on
the NIR data). We add an additional parameter ($b$), which accounts for both
rotational and instrumental broadening as a single Gaussian width. All
parameters were fit under uniform priors, with emcee using 100 walkers and
100,000 steps.

The results from the IGRINS orders were broadly consistent with that from our
SNIFS+SpeX data; the fit strongly prefers either
$T_{\rm{surf}}$$\simeq$$T_{\rm{spot}}$ or $f_{S}<20\%$. The major differences
from our SNIFS+SpeX results were that the IGRINS fit preferred a warmer
$T_{\rm{surf}}$ (and hence warmer $T_{\rm{spot}}$) and that the fit to IGRINS
data allowed for lower temperature ($<2500$ K) spots with higher surface
coverage fractions (up to $\simeq$20%) than the fit to the SNIFS+SpeX data.

A comparison between the posteriors from the fit to the stellar spectrum and
the fit to the transmission spectrum (Figure 9) had no overlap. The fit to the
stellar spectrum rules out all possible $T_{\rm{spot}}$ and $f_{S}$
combinations that could explain the observed transmission spectrum. This was
true regardless of the dataset used.

While the SNIFS+SpeX results are more precise, we used the fit from IGRINS for
further analysis. The main reason was that tests on similar spectra (same
instrument and setup) for other stars suggested that the SNIFS+SpeX spectra
were sensitive to the relative flux calibration. While both SNIFS and SpeX
have excellent flux calibration (Mann et al., 2015), K2-33 was observed at an
air mass of $>1.4$, around which the blue end of the (SNIFS) data become less
reliable (Buton et al., 2013). Further, IGRINS data was taken over many epochs
(yielding similar results regardless of which IGRINS spectrum we used) and the
data was contemporaneous with the MEarth and Spitzer transits. Overall, we
consider the IGRINS fit to be both more accurate and more conservative.

We show a more fair comparison of the spectroscopic results to the transit
data in Figure 8. For this, we randomly drew $T_{\rm{surf}}$ and
$T_{\rm{spot}}$ values from the (stellar spectrum fit) posterior and computed
the predicted transit depth as a function of wavelength, normalized to the
depth at 4.5$\mu$m. As can be seen, even the most extreme samplings yielded
transit depths inconsistent with the MEarth and K2 data.

For our fit to the transmission spectrum, we had one fewer free parameter, as
we fixed $T_{\rm{surf}}$ to 3540 K; for our fit to the stellar spectrum, both
$T_{\rm{surf}}$ and $T_{\rm{spot}}$ were free parameters. Fixing
$T_{\rm{surf}}$ in the spectroscopic fit or shifting the posterior on
$T_{\rm{surf}}$ (for the IGRINS fit) did not change the conclusion: the
spectroscopic fit was still inconsistent with the spot coverage required to
explain the observed transit depths.

The uncertainties from both fits were likely underestimated due to imperfect
models, but this is also unlikely to impact our results. For example, we tried
using the PHOENIX model atmospheres from Husser et al. (2013) which yielded
$T_{\rm{surf}}$ and $T_{\rm{spot}}$ values $\simeq 60$ K cooler and $A_{V}$
values 0.1 magnitudes lower. These changes were larger than the errors within
a fit, and are likely due to differences in the handling of strong molecular
lines in the spectra of cool stars (Mann et al., 2013; Passegger et al.,
2016). However, the conclusions on $f_{S}$ were largely unchanged;
$T_{\rm{spot}}$ must either be within $\simeq 150$ K of $T_{\rm{surf}}$ or
$f_{S}$ must be below $\simeq$20%. This was expected: these model systematics
impact both the surface and spot spectra in similar ways, keeping the
difference in temperature similar. We concluded that even in the presence of
significant systematics in the model of M dwarfs, the spectrum rules out any
spot pattern that can explain the transit depths.

It is likely that the stellar surface contains more than a single spot, and
hence more than a single spot temperature. However, adding additional spots
with different $T_{\text{eff}}$ values did not meaningfully change the
conclusion. The additional spot simply followed the output distribution in
$T_{\rm{spot}}$ and $f_{S}$ of the single spot assumed above, but with larger
uncertainties. It is possible that allowing an arbitrary number of spots each
with unconstrained temperatures can reproduce the observed transmission
spectrum. However, our tests suggest that such a reconciliation would come
primarily by increasing the uncertainties until the contours had marginal
overlap; in this case, the fit does not improve meaningfully with more spots.

##  6 Photochemical Hazes

In addition to spots, the existence of submicron aerosol particles in K2-33b’s
atmosphere could also lead to enhanced optical transit depths. One possible
formation pathway for such particles is photochemistry resulting in the
breakup of small molecules at low pressures ($\sim$1 $\mu$bar) and subsequent
polymerization of photochemical products, forming high altitude hazes (e.g.
Morley et al., 2015; Kawashima & Ikoma, 2019; Gao et al., 2021). The small
size of these particles causes a rapid decrease in opacity towards longer
wavelengths and thus a decrease in transit depth. This effect is strengthened
in young planets like K2-33b that may be experiencing atmospheric outflows,
which are capable of pushing haze particles to greater altitudes and
increasing the difference in transit depth between the optical and the NIR
(Wang & Dai, 2019; Gao & Zhang, 2020; Ohno & Tanaka, 2021).

We explore the impact of hazes on K2-33b’s transmission spectra following the
method described in Gao & Zhang (2020). Briefly, we construct model
atmospheres defined by user-selected core masses and atmospheric mass
fractions composed of a convective zone at depth and a radiative region at
lower pressures separated by the radiative-convection boundary, where the
temperature is set to the equilibrium temperature of the planet assuming zero
albedo and full heat redistribution; we use the Rosseland mean opacity from
Freedman et al. (2014) to find the radiative-convective boundary and to
calculate the temperature-pressure profile in the radiative zone. The core is
assumed to have Earth-like composition, with the mass-radius relationship from
Zeng et al. (2019). As an update to Gao & Zhang (2020), we use a planetary
evolution model (Tang et al., 2022; Lopez & Fortney, 2014) to derive an
intrinsic temperature as a function of core mass, yielding values
$\sim$100-150 K for core masses of 3-10$M_{\oplus}$, atmospheric mass
fractions $<$10%, and model planet ages matching that of K2-33b.

We consider the effect of atmospheric loss pushing haze particles upwards. We
compute atmospheric loss rates assuming energy limited escape with an energy-
deposition pressure level of 1 nbar and a mass-loss efficiency of 10% (e.g.
Lopez, 2017). We use the saturated stellar XUV flux for M dwarfs (Wright et
al., 2011, 2018; Pineda et al., 2021) to drive the atmospheric loss.

The elemental composition of the atmosphere is simplified to just H, He, O,
and C, with solar abundances of O and C incorporated into water vapor and
methane, respectively, and the rest of the atmosphere composed of H/He. We
initially choose methane as the dominant carbon carrier (“CH4 models”) due to
the relatively low temperature of K2-33b (Lodders & Fegley, 2002), but we also
explore the impact of upward mixed CO (“CO models”) on the transmission
spectra (see Fortney et al., 2020, and references therein). Given this
atmospheric composition and assumed fully mixed abundance profiles, we equate
the haze production rate profile to the methane photolysis rate profile
multiplied by a haze production efficiency factor, which is a free parameter.
The methane photolysis rate profile is generated using a simplified
photochemical scheme where methane is photolyzed by Lyman $\alpha$ radiation
from the host star, with an assumed Lyman $\alpha$ flux of a saturated M dwarf
(Pineda et al., 2021).

Once produced, haze particles are transported to the deep atmosphere via
sedimentation and mixing, the latter of which is parameterized using eddy
diffusion. The eddy diffusion coefficient is a free parameter and assumed to
be constant with altitude. The haze particles can grow through coagulation
during transport, starting from a fixed minimum radius of 10 nm, and we assume
that all haze particles are spherical. The haze evolution is simulated using
the Community Aerosol and Radiation Model for Atmospheres (CARMA; Turco et
al., 1979; Toon et al., 1989; Ackerman et al., 1995; Gao et al., 2018), which
has been applied to exoplanet and solar system hazes on numerous occasions
(Gao et al., 2017; Adams et al., 2019; Gao & Zhang, 2020; Gao et al., 2020).
We use Mie scattering to calculate the haze optical properties and consider
two haze compositions: soot and tholin, with the corresponding complex
refractive indices (Lee & Tien, 1981; Khare et al., 1984; Chang &
Charalampopoulos, 1990; Morley et al., 2015; Gavilan et al., 2016; Lavvas &
Koskinen, 2017). Finally, we combine the model atmosphere and haze optical
properties to compute the model transmission spectra.

###  6.1 CH4 models

Since the mass of the planet is unknown, we considered several planet core
masses. From the planet mass-radius distribution, we calculated a predicted
mass for K2-33b of $16.59_{-11.73}^{+7.08}M_{\oplus}$ using the python
package, forecastor (Chen & Kipping, 2017). Due to the youth of the planet, we
expect the true value to be smaller because younger planets are larger and
less dense than their evolved counterparts (Owen, 2020). As such, we tested
core masses of 3, 5, and 8$M_{\oplus}$. We also examined eddy diffusion
coefficients of 109 and 1011 cm2 s-1 and haze production efficiencies of 1, 2,
4, and 8% (Lavvas & Koskinen, 2017). For each core mass, we varied the
atmospheric mass fraction to roughly match the observed planet radius,
resulting in values of 1–2.5%.

To compare our models to the data, we convolved the model spectra with the
relevant filter profiles to create synthetic transit depths corresponding to
each effective wavelength (photon weighted mean wavelength). The effective
wavelength factors in the widths of the broadband filters were calculated
using K2-33b’s spectrum and each filter’s bandpass. The results of this
calculation yielded effective wavelengths of 0.72 $\mu$m (K2), 0.84 $\mu$m
(MEarth), 3.46 $\mu$m (Channel 1), 4.43 $\mu$m (Channel 2). We added a free
parameter, a normalization constant, to allow each model spectra to shift in
median depth, and varied it to minimize the $\chi^{2}$ when compared to our
data. Although the Gao & Zhang (2020) model computes the radius of the model
planets, we allow this normalization factor to shift the planet radius to
account for model uncertainties in e.g., the intrinsic temperature,
atmospheric composition, etc.

Our results disfavored a featureless model and all haze models where the
carbon carrier in the background atmosphere is CH4 ($>4\sigma$ confidence)
regardless of the haze composition, core mass, atmospheric mass fraction, eddy
diffusion coefficient, and haze production efficiency. The lowest $\chi^{2}$
for these models that had a normalization factor within 20% was $\chi^{2}=$ 58
for tholin hazes and $\chi^{2}=$ 79 for soot hazes. These high values are
driven by the disagreement between the models and the Spitzer points. The
normalization parameters and the corresponding $\chi^{2}$ values are listed in
Table 6. Only the top three lowest $\chi^{2}$ values are listed for each haze
composition.

###  6.2 CO models

Under the assumption of thermochemical equilibrium and solar metallicity,
methane should be the dominant carbon carrier in the atmospheres of planets
with K2-33b’s equilibrium temperature (Lodders & Fegley, 2002). However, no
methane has been detected in the atmospheres of exoplanets that are similar in
temperature, but older (Kreidberg et al., 2018; Benneke et al., 2019). One
possible reason for this is the disequilibrium process of mixing bringing CO-
rich gas from depth to the observable part of the atmosphere (Fortney et al.,
2020). To test this scenario and better fit the Spitzer points, we replace CH4
in our previous models with CO and reduce the water vapor abundance
accordingly. The same core masses, atmospheric mass fractions, eddy diffusion
coefficients, and haze production efficiencies were considered in generating
the CO model grid. We use the same haze production rate even though it is
computed assuming methane photolysis so as to isolate the impact of changing
atmospheric composition on the transmission spectrum. Given the unknowns in
the haze production pathway and experimental evidence pointing to CO as a
possible haze precursor (Hörst et al., 2018; He et al., 2018; Fleury et al.,
2019), this simplification is justified.

Table 6 lists the top 3 lowest $\chi^{2}$ and normalization factors for each
haze composition. Our results disfavor hazes composed of soot to $>5\sigma$
confidence, as they are unable to fit the optical data. In contrast, tholin
hazes ($\chi^{2}=29$, dof=16) give reasonable fits to our measured
transmission spectrum ($>2\sigma$ confidence).

To summarize, consideration of photochemical hazes allowed us to put the
following constraints on our model parameters:

  * •

planet core mass: Lower values (3 or 5$M_{\oplus}$) preferred.

  * •

atmospheric mass fraction: 1–2.5% preferred.

  * •

eddy diffusion coefficient: Higher values (1011 cm2 s-1) preferred for tholin
hazes, while lower values (109 cm2 s-1) preferred for soot hazes.

  * •

haze efficiencies: Soot models prefer lower values (1%), while tholin models
prefer higher values (4%)

  * •

haze composition: Tholin strongly preferred over soot.

  * •

primary carbon carrier: CO strongly preferred over CH4.

The best-fit model from our grid has a 5$M_{\oplus}$ core and an atmospheric
mass fraction of 2.2–2.3%. The model atmosphere possesses a high (1011 cm2
s-1) eddy diffusion coefficient, a haze production efficiency of 4% (about
midway between that of Jupiter and Titan; Lavvas & Koskinen, 2017), and tholin
hazes, with CO as the dominant gaseous carbon carrier.

Table 6: Normalization Factor for Atmospheric Models

Haze | Planet Core Mass | Atmospheric Mass | Eddy Diffusion | Haze Production | Normalization |  $\chi^{2}$ aa There were 14 total consecutive transits taken by K2
---|---|---|---|---|---|---
Composition | [$M_{\oplus}$] | Fraction | Coefficient [cm2 s-1] | Efficiency [%] | Factor |
_Haze Fit with No Spots_
| 5.0 | 0.024 | 1011 | 4.0 | 0.87 | 58
CH4 Tholin | 5.0 | 0.025 | 1011 | 4.0 | 0.81 | 58
| 5.0 | 0.023 | 1011 | 4.0 | 0.93 | 59
| 3.0 | 0.013 | 109 | 1.0 | 0.84 | 79
CH4 Soot | 3.0 | 0.012 | 109 | 1.0 | 0.80 | 80
| 3.0 | 0.013 | 1011 | 1.0 | 0.86 | 84
| 5.0 | 0.022 | 1011 | 4.0 | 1.07 | 29
CO Tholin | 5.0 | 0.023 | 1011 | 4.0 | 1.00 | 29
| 5.0 | 0.021 | 1011 | 4.0 | 1.14 | 30
| 3.0 | 0.013 | 109 | 1.0 | 0.88 | 65
CO Soot | 3.0 | 0.014 | 109 | 1.0 | 0.83 | 68
| 3.0 | 0.015 | 109 | 1.0 | 0.80 | 68
Flat Line | $\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$ | 122
_Haze Fit with Spots_
| 5.0 | 0.020 | 1011 | 4.0 | 1.01 | 28
CO Tholin | 5.0 | 0.021 | 1011 | 4.0 | 0.98 | 28
| 5.0 | 0.019 | 1011 | 4.0 | 1.04 | 28
| 3.0 | 0.010 | 1011 | 1.0 | 1.12 | 56
CO Soot | 3.0 | 0.011 | 1011 | 1.0 | 1.00 | 59
| 3.0 | 0.012 | 1011 | 2.0 | 0.80 | 61

aafootnotetext: The spot-free haze fits has a dof= 16 and the combined haze
and spot fits has a dof=14.

Note. — Only the top 3 (lowest $\chi^{2}$ values) for each haze composition is
listed here. In addition to this, for the combined haze and spots fits, the
fit also needed to yield spot properties consistent with the fit to the
stellar spectrum.

Figure 12: Transmission spectrum of K2-33b from our data (gray circles)
compared to haze models composed of 1) CO Tholin (green), 2) CO Soot (pink),
3) a combined spot model and CO Tholin (purple), and 4) a combined spot model
and CO Soot (orange). Only the models that produced the lowest chi-square
values are shown. All models were normalized to give the best fit to the data.
The combined haze and spot models gave a lower chi-square value. All
calculations were done with high-resolution models.  Figure 13:  Posteriors of
spot properties from the stellar spectrum (black/grey), a fit to the
transmission spectrum using CO Tholin models (purple), and a fit to the
transmission spectrum using CO soot models (orange). The two transmission
spectrum fits shown here are those using the Tholin or Soot model with the
lowest $\chi^{2}$ that also yielded spot properties consistent with the fit to
the stellar spectrum.

###  6.3 Combining spots and hazes

As we showed in Figure 8, spots alone did a relatively poor job fitting the
overall transit depths, and the best fits required spot properties
inconsistent with the stellar spectrum. Photochemical hazes provided a better
fit to the data, but there is still some tension. Here, we endeavored to
improve the fit by combining hazes and spots.

Our methodology was the same as described in Section 5.2.1, except we replaced
the flat model with the theoretical haze models from Section 6 and $D$ was
changed to a dimensionless scale factor applied to the models (denoted
$D_{\rm{mod}}$). The scale factor can be an arbitrary number, but changing the
planet size $>$ 20 % invalidates the model assumption, so we limited the scale
factor to $0.8<D_{\rm{mod}}<1.2$; more extreme scaling would imply a different
surface gravity and hence invalidate the model atmosphere. The two other
parameters ($T_{\rm{spot}}$ and $f_{S}$) evolved under uniform priors. For
each model, we ran with 50 walkers and as many steps as required to pass 50
times the autocorrelation time (checking every 1000 steps). The number of
steps required varied from 4,000 to more than 15,000, depending on the model.

As with the fits ignoring effects from stellar spots, the models with CO as
the primary carbon carrier and tholin hazes gave the best fit. This is because
the CO models are able to match the shallower transit at 3.6$\mu$m compared to
that at 4.5$\mu$m. Spots alone cannot explain this difference, as their effect
is weak past 3$\mu$m and even strong spots tend to predict a deeper transit at
3.6$\mu$m instead of a shallower one.

Overall, adding spots provided only marginal improvement when fitting the
transmission spectrum. For example, the best fit haze model (both with and
without spots) was the CO model with a tholin haze (which we show in Figure
12), with an improvement of $\Delta\chi^{2}\simeq 1$ over a spot-free haze
model. Given that there were two additional free parameters in the spot and
haze fit (compared to haze-only), this improvement was not significant for
most of the tholin models (see Table 6 for a summary of the best fits).
However, it is worth highlighting that the joint fit of the spots and CO
Tholin model makes everything consistent: the fit to the stellar spectrum, the
fit to the light curve, and the fit to the transmission spectrum. The required
small spot coverage fractions and spot properties were consistent with the
stellar spectrum (Figure 13).

Including the effect of spots had a larger improvement for fits using the soot
models, with the best cases going from $\chi^{2}\simeq 65$ to
$\chi^{2}\simeq$35\. However, these generally required spot coverage fractions
and temperatures inconsistent with those inferred from the stellar spectrum.
If we consider just cases where the required spot properties were consistent
with the stellar spectroscopic data, the best fit soot haze and spot
combination did a poor job reproducing the full transmission spectrum
($\chi^{2}=56$).

##  7 Summary and Discussion

To explore the transmission spectrum of the $\simeq$10 Myr, K2-33b, we
combined 33 transit observations taken by K2, MEarth, HST, and Spitzer
spanning over $>$ 2.5 years. We also updated the stellar parameters based on
the parallax from Gaia.

The most striking result from the multi-wavelength data was that the optical
transits from K2 and MEarth are almost a factor of two deeper than the NIR
transits from HST and Spitzer (0.24% vs. 0.13%). We explore whether this depth
difference is due to unconstrained systematics in the datasets or if it is
astrophysical in nature. We rule out the first scenario primarily because the
difference holds across multiple datasets, with roughly consistent (but
shallower) depths from Spitzer and HST and consistent (but deeper) depths from
K2 and MEarth. Agreement between instruments rules out issues related to the
data source, data quality, or analysis method. The Spitzer, MEarth, and K2
data each cover at least five transits and span numerous rotational cycles,
ruling out single events like flares. Further, three transits were observed
simultaneously by MEarth and Spitzer (transit numbers 186, 188, and 191; Table
1). Even considering just these three transits, the MEarth and Spitzer depths
were still clearly inconsistent (Figure 14). This rules out long-term
evolutionary effects.

Figure 14: Light curves from MEarth and Spitzer for the three transits with
overlapping data in both (transits 186, 188, and 191). Data from both
facilities were folded and (median) binned for clarity. One of the MEarth
transits was partial (transit 188) and conditions varied both between and
during MEarth transits leading to inconsistent data quality. However, the
MEarth light curve still clearly favors a deeper transit depth.

In our analysis, we excluded two Spitzer and two MEarth transits. The two
MEarth transits were removed primarily because of poor coverage (insufficient
out-of-transit baseline). Including them had negligible effect on the transit
depth because the GP was too poorly constrained. As we show in Figure 4, the
two excluded Spitzer transits favored deeper transits, although Channel 1 data
this was due to a visible flare before the transit. As a test, we repeated our
fit using MISTTBORN including all transits but masking out the data right
around the flare (in AOR 60658432). Otherwise, the analysis was identical to
that described in Section 4.2. For channel 1, this gave $R_{P}/R_{*}=0.0358\pm
0.0010$ and for channel 2, $R_{P}/R_{*}=0.0383\pm 0.0012$. Neither of the two
fits were inconsistent with our original analysis and both were still
$>5\sigma$ inconsistent with the depth from K2.

We conclude that the depth difference is astrophysical. We consider two
causes: (dark) spots on the stellar surface, or hazes in the planet’s
atmosphere – both of which predict deeper transits in the optical. In the case
of star spots, their effect is largest at the optical wavelengths because the
contrast ratio of stellar spots decreases with increasing wavelength. As for
the hazes, their small particle size results in a decrease in opacity with
increasing wavelength, and thus a decrease in transit depth.

We tested three scenarios; 1) a spotted star with a (perfectly) flat planetary
transmission spectrum, 2) a mix of photochemical hazes with a pristine star,
and 3) a mix of both spots and hazes. Here is a summary of each fit:

  * •

Spots:  The deeper optical transits can be explained if $>$60 % of the star is
covered in cold ($<$3200 K) spots. However, such large spot coverage fractions
are ruled out by the stellar spectrum, which indicate the spots must either
cover $<20\%$ of the surface or be similar in temperature to the surface
($>$3,300 K).

  * •

Photochemical Hazes:  A planet with a low core mass ($\simeq 5M_{\oplus}$) and
a tholin haze provides the best fit to the data ($\chi^{2}=29$ with 16 degrees
of freedom). All soot models provided a poorer fit than tholin models. This is
because tholins provide more absorption in the blue optical and UV than in the
NIR, resulting in a steeper optical slope, whereas soots are more effective at
absorbing at all wavelengths leading to a more gentle slope (Figure 12). For
any haze, models with CO as the carbon carrier are required to explain the
difference between the two Spitzer bands. The low inferred planet core mass is
required to explain the large difference between the optical and NIR transit
depths, as higher mass cores would result in reduced atmospheric scale
heights.

  * •

Spots + Photochemical Hazes:  A combined fit offers only marginal improvement
from the spot-free haze fits. For tholin hazes, the improvement is not
significant, and for soot hazes, the spot fraction required is inconsistent
with the stellar spectrum and the fit to the transmission spectrum is still
poor (high $\chi^{2}$) compared to the tholin-only fit. It is, however,
encouraging that the resulting spot properties when using the tholin model are
consistent with the two-temperature fit to the stellar spectrum. Adding spots
does not change any of our conclusions with respect to the CO vs CH4 models or
low planet mass.

The complication to our current analysis is that spot coverage fractions may
change with time. The K2 and spectroscopy data were used to constrain the
spots, but K2 data predate the other sets of data by $\sim$ 1.5 years.
However, removing the Kepler data did not change the conclusions and the
improvement to the fit was marginal. The IGRINS spectra were taken
contemporaneously with the MEarth data and multiple spectra taken over several
months (2016 Jan 25 to 2016 Mar 28) yield consistent results. There is also no
evidence of extreme spot evolution in the long-term MEarth monitoring. Thus,
these kinds of long-term effects are unlikely to change our findings, but may
explain some of the tension in the tholin-haze fit.

A possibility we did not explore was a transiting ringed system (Mamajek et
al., 2012). Ohno & Fortney (2022) highlight that even close-in planets can
retain thick rings at young ($<100$ Myr) ages, and rings would help explain
the larger radii observed for young planets. Rings on K2-33 b would need to be
dominated by $\simeq$micron-sized dust and extend to at least twice the
planetary radius to explain the deeper optical transits. Depending on the
orientation and size distribution, rings would dominate over any signal from
the planet’s atmosphere (Alam et al., 2022), which would explain the lack of
spectral features seen here. A companion paper, Ohno (in press), investigates
this possibility further.

We conclude from our analysis that K2-33b likely hosts tholin-like hazes, with
modest ($<10\%$) spot coverage for its host star. This scenario can explain
all the data, including the light curve, stellar spectrum, and planetary
transmission spectrum. More detailed information, including better constraints
on the spot properties and haze production rate, will require more data. The
current HST data is particularly limiting, since it was only half a transit. A
full transit of HST/WFC3 combined with contemporaneous ground-based transit
and spectroscopic monitoring would let us fully disentangle the effects of
spot surface inhomogeneities from haze properties in the observed transmission
spectrum. These observations will break this degeneracy because stellar spots
generally increase the strength of the NIR H2O feature, while hazes generally
weaken it. Further, the former changes the observed out-of-transit spectrum
while the latter does not, and each scenario will have a significantly
different impact on the transit depths at bluer wavelengths depending on the
spot and haze properties. Contemporaneous data would also resolve issues of
spot properties changing with time and enable us to do a simultaneous fit of
all constraints instead of the one-at-a-time fit we did here.

Our haze results are in line with those of Ohno & Kawashima (2020), which
suggested that moderate haze production rates combined with high values of the
eddy diffusion coefficient lead to super-Rayleigh optical slopes in
transmission spectra. While both of these parameters are loosely constrained
in exoplanet atmospheres, the high eddy diffusion coefficient may be
reasonable given the young age and therefore high luminosity of K2-33b, and
potentially tying into the significant upwelling of CO-rich gas suggested by
the Spitzer data. Our results also fit those of Gao & Zhang (2020), which
predicted significant enlargement of young planet transit radii due to high
altitude hazes entrained in atmospheric outflows, though here the effect seems
to impact the optical wavelengths much more than the NIR. Both the core mass
and atmospheric mass fraction implied by our results are similar to those of
typical sub-Neptunes (Owen & Wu, 2017), suggesting that K2-33b will eventually
join that population after thermal evolution and contraction.

The need for CO to be the carbon carrier in our best-fit models (regardless of
spot levels) is extremely compelling, but it relies entirely on the Spitzer
bands. While the difference in depth holds even when we try different analysis
methods, the significance of the difference changes depending on how we handle
combining multiple transits (the difference is not significant with a single
transit).

Upward mixing of CO is only one way in which it can supersede CH4 as the
dominant carbon carrier in an atmosphere at the temperature of K2-33b.
Alternatively, K2-33b’s atmosphere can possess a low C/O and an enhanced
metallicity, which would provide clues to its formation and early evolution
(e.g., Öberg et al., 2011; Madhusudhan, 2012). Due to the large scale heights
required to match the observed transmission spectrum, the metallicity cannot
be too high (e.g. $>$100 $\times$ solar). Detecting a water and/or CO2 feature
would help differentiate between these scenarios. We eagerly await
confirmation of this detection with HST and/or JWST, which would enable a more
precise measurement of the carbon-to-oxygen ratio.

We thank the anonymous referee for their careful reading and thoughtful
comments on the manuscript. The authors would like to thank Laura Kreidberg
and Alex Teachey for helpful conversations regarding the HST reduction. The
authors also wish to acknowledge Wally, Penny, and Bandit for their input,
sound counsel, and their constant quest to discover the unknown. P.C.T was
supported by NSF Graduate Research Fellowship (DGE-1650116), the NC Space
Grant Graduate Research Fellowship, the Zonta International Amelia Earhart
Fellowship, and the Jack Kent Cooke Foundation Graduate Scholarship. This work
was made possible by grants to A.W.M. from the K2 Guest Observer Program
(80NSSC19K0097) and NASA’s Astrophysics Data Analysis Program (80NSSC19K0583).
M.J.F was supported by the NC Space Grant Graduate Research program and a
grant from NASA’s Exoplanet Research Program (80NSSC21K0393). This paper
includes data collected by the K2 mission. Funding for the K2 mission is
provided by the NASA Science Mission directorate. This work is based on
observations made with the Spitzer Space Telescope, which is operated by the
Jet Propulsion Laboratory, California Institute of Technology under a contract
with NASA. Support for this work was provided by NASA through an award issued
by JPL/Caltech. This research is based on observations made with the NASA/ESA
Hubble Space Telescope obtained from the Space Telescope Science Institute,
which is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5–26555. These observations are
associated with program ID: 14887 (PI: Benneke). The MEarth Team gratefully
acknowledges funding from the David and Lucile Packard Fellowship for Science
and Engineering (awarded to D.C.). This material is based upon work supported
by the National Science Foundation under grant AST-1616624, and work supported
by the National Aeronautics and Space Administration under Grant No.
80NSSC18K0476 issued through the XRP Program.

## References

  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42
  * Ackerman et al. (1995) Ackerman, A. S., Toon, O. B., & Hobbs, P. V. 1995, J. Geophys. Res., 100, 7121
  * Adams et al. (2019) Adams, D., Gao, P., de Pater, I., & Morley, C. V. 2019, ApJ, 874, 61
  * Alam et al. (2022) Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Royal Society of London Philosophical Transactions Series A, 370, 2765
  * Allard et al. (2013) Allard, F., Homeier, D., Freytag, B., et al. 2013, Memorie della Societa Astronomica Italiana Supplementi, 24, 128
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Barclay et al. (2021) Barclay, T., Kostov, V. B., Colón, K. D., et al. 2021, arXiv e-prints, arXiv:2109.14608
  * Benneke et al. (2019) Benneke, B., Knutson, H. A., Lothringer, J., et al. 2019, Nature Astronomy, 3, 813
  * Benneke et al. (2019) Benneke, B., Wong, I., Piaulet, C., et al. 2019, The Astrophysical Journal Letters, 887, L14
  * Berta et al. (2012) Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, ApJ, 747, 35
  * Booth et al. (2017) Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, Monthly Notices of the Royal Astronomical Society, 469, 3994
  * Booth & Ilee (2019) Booth, R. A., & Ilee, J. D. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3998
  * Bouma et al. (2020) Bouma, L. G., Hartman, J. D., Brahm, R., et al. 2020, AJ, 160, 239
  * Brown (2001) Brown, T. M. 2001, The Astrophysical Journal, 553, 1006
  * Buton et al. (2013) Buton, C., Copin, Y., Aldering, G., et al. 2013, A&A, 549, A8
  * Campo et al. (2011) Campo, C. J., Harrington, J., Hardy, R. A., et al. 2011, ApJ, 727, 125
  * Chang & Charalampopoulos (1990) Chang, H., & Charalampopoulos, T. T. 1990, Proceedings of the Royal Society of London Series A, 430, 577
  * Chen & Kipping (2017) Chen, J., & Kipping, D. 2017, ApJ, 834, 17
  * Chilcote et al. (2017) Chilcote, J., Pueyo, L., De Rosa, R. J., et al. 2017, AJ, 153, 182
  * Claret (2000) Claret, A. 2000, A&A, 363, 1081
  * Cubillos et al. (2017) Cubillos, P., Erkaev, N. V., Juvan, I., et al. 2017, Monthly Notices of the Royal Astronomical Society, 466, 1868
  * Cutri & et al. (2014) Cutri, R. M., & et al. 2014, VizieR Online Data Catalog, II/328
  * David et al. (2019) David, T. J., Petigura, E. A., Luger, R., et al. 2019, ApJ, 885, L12
  * David et al. (2016a) David, T. J., Hillenbrand, L. A., Petigura, E. A., et al. 2016a, Nature, 534, 658
  * David et al. (2016b) David, T. J., Conroy, K. E., Hillenbrand, L. A., et al. 2016b, AJ, 151, 112
  * de Zeeuw et al. (1999) de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354
  * Dittmann et al. (2016) Dittmann, J. A., Irwin, J. M., Charbonneau, D., & Newton, E. R. 2016, ApJ, 818, 153
  * Donati et al. (2017) Donati, J.-F., Yu, L., Moutou, C., et al. 2017, MNRAS, 465, 3343
  * Doyle et al. (2018) Doyle, L., Ramsay, G., Doyle, J. G., Wu, K., & Scullion, E. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2153
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4
  * Fang et al. (2016) Fang, X.-S., Zhao, G., Zhao, J.-K., Chen, Y.-Q., & Bharat Kumar, Y. 2016, MNRAS, 463, 2494
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, The Astrophysical Journal Supplement Series, 154, 10
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99
  * Fleury et al. (2019) Fleury, B., Gudipati, M. S., Henderson, B. L., & Swain, M. 2019, ApJ, 871, 158
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fortney et al. (2020) Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288
  * Freedman et al. (2014) Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25
  * Gao et al. (2018) Gao, P., Marley, M. S., & Ackerman, A. S. 2018, ApJ, 855, 86
  * Gao et al. (2021) Gao, P., Wakeford, H. R., Moran, S. E., & Parmentier, V. 2021, Journal of Geophysical Research (Planets), 126, e06655
  * Gao & Zhang (2020) Gao, P., & Zhang, X. 2020, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 52, AAS/Division for Planetary Sciences Meeting Abstracts, 213.03
  * Gao et al. (2017) Gao, P., Fan, S., Wong, M. L., et al. 2017, Icarus, 287, 116
  * Gao et al. (2020) Gao, P., Thorngren, D. P., Lee, E. K. H., et al. 2020, Nature Astronomy, 4, 951
  * Gavilan et al. (2016) Gavilan, L., Alata, I., Le, K. C., et al. 2016, A&A, 586, A106
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, Vol. 5, No. 1, p. 65-80, 2010, 5, 65
  * Grunblatt et al. (2016) Grunblatt, S. K., Huber, D., Gaidos, E. J., et al. 2016, AJ, 152, 185
  * Gully-Santiago et al. (2017) Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836, 200
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * He et al. (2018) He, C., Hörst, S. M., Lewis, N. K., et al. 2018, AJ, 156, 38
  * Henden et al. (2012) Henden, A. A., Levine, S. E., Terrell, D., Smith, T. C., & Welch, D. 2012, Journal of the American Association of Variable Star Observers (JAAVSO), 40, 430
  * Hörst et al. (2018) Hörst, S. M., He, C., Lewis, N. K., et al. 2018, Nature Astronomy, 2, 303
  * Howard et al. (2014) Howard, A. W., Marcy, G. W., Fischer, D. A., et al. 2014, ApJ, 794, 51
  * Hubbard et al. (2001) Hubbard, W., Fortney, J., Lunine, J., et al. 2001, The Astrophysical Journal, 560, 413
  * Huitson et al. (2013) Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Ingalls et al. (2012) Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2012, in Proc. SPIE, Vol. 8442, Space Telescopes and Instrumentation 2012: Optical, Infrared, and Millimeter Wave, 84421Y
  * Ingalls et al. (2016) Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2016, AJ, 152, 44
  * Irwin et al. (2007) Irwin, J., Irwin, M., Aigrain, S., et al. 2007, MNRAS, 375, 1449
  * Irwin et al. (2015) Irwin, J. M., Berta-Thompson, Z. K., Charbonneau, D., et al. 2015, in Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, Vol. 18, 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. G. T. van Belle & H. C. Harris, 767–772
  * Johns-Krull et al. (2016) Johns-Krull, C. M., McLane, J. N., Prato, L., et al. 2016, ApJ, 826, 206
  * Johnson et al. (2021) Johnson, L. J., Norris, C. M., Unruh, Y. C., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 4751
  * Johnson et al. (2018) Johnson, M. C., Dai, F., Justesen, A. B., et al. 2018, MNRAS, 481, 596
  * Kawashima & Ikoma (2019) Kawashima, Y., & Ikoma, M. 2019, ApJ, 877, 109
  * Khare et al. (1984) Khare, B. N., Sagan, C., Arakawa, K. T., et al. 1984, Icarus, 60, 127
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161
  * Kreidberg et al. (2018) Kreidberg, L., Line, M. R., Thorngren, D., Morley, C. V., & Stevenson, K. B. 2018, ApJ, 858, L6
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69
  * Kreidberg et al. (2015) Kreidberg, L., Line, M. R., Bean, J. L., et al. 2015, ApJ, 814, 66
  * Lammer et al. (2016) Lammer, H., Erkaev, N., Fossati, L., et al. 2016, Monthly Notices of the Royal Astronomical Society: Letters, 461, L62
  * Lantz et al. (2004) Lantz, B., Aldering, G., Antilogus, P., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5249, Optical Design and Engineering, ed. L. Mazuray, P. J. Rogers, & R. Wartmann, 146–155
  * Lavvas & Koskinen (2017) Lavvas, P., & Koskinen, T. 2017, ApJ, 847, 32
  * Lee et al. (2013) Lee, J.-M., Heng, K., & Irwin, P. G. J. 2013, ApJ, 778, 97
  * Lee & Tien (1981) Lee, S., & Tien, C. 1981, Symposium (International) on Combustion, 18, 1159, eighteenth Symposium (International) on Combustion
  * Lewis et al. (2013) Lewis, N. K., Knutson, H. A., Showman, A. P., et al. 2013, ApJ, 766, 95
  * Libby-Roberts et al. (2020) Libby-Roberts, J. E., Berta-Thompson, Z. K., Désert, J.-M., et al. 2020, AJ, 159, 57
  * Libby-Roberts et al. (2021) Libby-Roberts, J. E., Berta-Thompson, Z. K., Diamond-Lowe, H., et al. 2021, arXiv e-prints, arXiv:2105.10487
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2
  * Lodders & Fegley (2002) Lodders, K., & Fegley, B. 2002, Icarus, 155, 393
  * Lopez (2017) Lopez, E. D. 2017, MNRAS, 472, 245
  * Lopez & Fortney (2014) Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1
  * Luger et al. (2021a) Luger, R., Foreman-Mackey, D., & Hedges, C. 2021a, AJ, 162, 124
  * Luger et al. (2021b) Luger, R., Foreman-Mackey, D., Hedges, C., & Hogg, D. W. 2021b, arXiv e-prints, arXiv:2102.00007
  * Madhusudhan (2012) Madhusudhan, N. 2012, ApJ, 758, 36
  * Mamajek et al. (2012) Mamajek, E. E., Quillen, A. C., Pecaut, M. J., et al. 2012, AJ, 143, 72
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171
  * Mandell et al. (2013) Mandell, A. M., Haynes, K., Sinukoff, E., et al. 2013, ApJ, 779, 128
  * Mann et al. (2015) Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64
  * Mann et al. (2013) Mann, A. W., Gaidos, E., & Ansdell, M. 2013, ApJ, 779, 188
  * Mann & von Braun (2015) Mann, A. W., & von Braun, K. 2015, PASP, 127, 102
  * Mann et al. (2016a) Mann, A. W., Gaidos, E., Mace, G. N., et al. 2016a, ApJ, 818, 46
  * Mann et al. (2016b) Mann, A. W., Newton, E. R., Rizzuto, A. C., et al. 2016b, AJ, 152, 61
  * May & Stevenson (2020) May, E. M., & Stevenson, K. B. 2020, AJ, 160, 140
  * Morley et al. (2015) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2015, ApJ, 815, 110
  * Morris (2020) Morris, B. M. 2020, ApJ, 893, 67
  * Morris et al. (2019) Morris, B. M., Curtis, J. L., Sakari, C., Hawley, S. L., & Agol, E. 2019, AJ, 158, 101
  * Morris et al. (2017) Morris, B. M., Hebb, L., Davenport, J. R. A., Rohn, G., & Hawley, S. L. 2017, ApJ, 846, 99
  * NASA Exoplanet Science Institute (2020) NASA Exoplanet Science Institute. 2020, Planetary Systems Table, doi:10.26133/NEA12
  * Newton et al. (2019) Newton, E. R., Mann, A. W., Tofflemire, B. M., et al. 2019, ApJ, 880, L17
  * Nicholson & Aigrain (2022) Nicholson, B. A., & Aigrain, S. 2022, Monthly Notices of the Royal Astronomical Society, 515, 5251
  * Nutzman & Charbonneau (2008) Nutzman, P., & Charbonneau, D. 2008, PASP, 120, 317
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16
  * Obermeier et al. (2016) Obermeier, C., Henning, T., Schlieder, J. E., et al. 2016, AJ, 152, 223
  * Ohno (in press) Ohno, K. in press
  * Ohno & Fortney (2022) Ohno, K., & Fortney, J. J. 2022, ApJ, 930, 50
  * Ohno & Kawashima (2020) Ohno, K., & Kawashima, Y. 2020, ApJ, 895, L47
  * Ohno & Tanaka (2021) Ohno, K., & Tanaka, Y. A. 2021, ApJ, 920, 124
  * Owen (2020) Owen, J. E. 2020, Monthly Notices of the Royal Astronomical Society, 498, 5030
  * Owen & Wu (2017) Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29
  * Parviainen & Aigrain (2015) Parviainen, H., & Aigrain, S. 2015, MNRAS, 453, 3821
  * Passegger et al. (2016) Passegger, V. M., Wende-von Berg, S., & Reiners, A. 2016, A&A, 587, A19
  * Petrov et al. (1994) Petrov, P. P., Shcherbakov, V. A., Berdyugina, S. V., et al. 1994, A&AS, 107, 9
  * Pineda et al. (2021) Pineda, J. S., Youngblood, A., & France, K. 2021, ApJ, 911, 111
  * Rackham et al. (2017) Rackham, B., Espinoza, N., Apai, D., et al. 2017, The Astrophysical Journal, 834, 151
  * Rayner et al. (2009) Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, ApJS, 185, 289
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362
  * Ricker (2014) Ricker, G. R. 2014, Journal of the American Association of Variable Star Observers (JAAVSO), 42, 234
  * Rizzuto et al. (2017) Rizzuto, A. C., Mann, A. W., Vanderburg, A., Kraus, A. L., & Covey, K. R. 2017, AJ, 154, 224
  * Roettenbacher & Vida (2018) Roettenbacher, R. M., & Vida, K. 2018, The Astrophysical Journal, 868, 3
  * Rogers (2015) Rogers, L. A. 2015, ApJ, 801, 41
  * Samland et al. (2017) Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57
  * Seager & Sasselov (2000) Seager, S., & Sasselov, D. 2000, The Astrophysical Journal, 537, 916
  * Sing et al. (2015) Sing, D. K., Wakeford, H. R., Showman, A. P., et al. 2015, MNRAS, 446, 2428
  * Sing et al. (2016) Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Stevenson et al. (2012) Stevenson, K. B., Harrington, J., Fortney, J. J., et al. 2012, ApJ, 754, 136
  * Tanaka & Ward (2004) Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388
  * Tang et al. (2022) Tang, Y., Fortney, J., Murray-Clay, R., Thorngren, D., & Broome, M. 2022, in Bulletin of the American Astronomical Society, Vol. 54, 100.03
  * Thao et al. (2020) Thao, P. C., Mann, A. W., Johnson, M. C., et al. 2020, AJ, 159, 32
  * Toon et al. (1989) Toon, O. B., Turco, R. P., Jordan, J., Goodman, J., & Ferry, G. 1989, Journal of Geophysical Research, 94, 11359
  * Tsiaras et al. (2016) Tsiaras, A., Waldmann, I. P., Rocchetto, M., et al. 2016, ApJ, 832, 202
  * Tsiaras et al. (2019) Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J., & Yurchenko, S. N. 2019, Nature Astronomy, 3, 1086
  * Tsiaras et al. (2018) Tsiaras, A., Waldmann, I., Zingales, T., et al. 2018, The Astronomical Journal, 155, 156
  * Turco et al. (1979) Turco, R. P., Hamill, P., Toon, O. B., Whitten, R. C., & Kiang, C. S. 1979, Journal of the Atmospheric Sciences, 36, 699
  * Vanderburg & Johnson (2014) Vanderburg, A., & Johnson, J. A. 2014, PASP, 126, 948
  * Vanderburg et al. (2016) Vanderburg, A., Latham, D. W., Buchhave, L. A., et al. 2016, ApJS, 222, 14
  * Wang & Dai (2019) Wang, L., & Dai, F. 2019, ApJ, 873, L1
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, The Astrophysical Journal Supplement Series, 154, 1
  * Wright et al. (2011) Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48
  * Wright et al. (2018) Wright, N. J., Newton, E. R., Williams, P. K. G., Drake, J. J., & Yadav, R. K. 2018, MNRAS, 479, 2351
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723