Tentative Evidence for Water Vapor in the Atmosphere of the Neptune-Size Exoplanet HD 106315 c

We present a transmission spectrum for the Neptune-size exoplanet HD 106315 c from optical to infrared wavelengths based on transit observations from the Hubble Space Telescope/Wide Field Camera 3, K2, and Spitzer. The spectrum shows tentative evidence for a water absorption feature in the $1.1 - 1.7\mu$m wavelength range with a small amplitude of 30 ppm (corresponding to just $0.8 \pm 0.04$ atmospheric scale heights). Based on an atmospheric retrieval analysis, the presence of water vapor is tentatively favored with a Bayes factor of 1.7 - 2.6 (depending on prior assumptions). The spectrum is most consistent with either enhanced metallicity, high altitude condensates, or both. Cloud-free solar composition atmospheres are ruled out at $>5\sigma$ confidence. We compare the spectrum to grids of cloudy and hazy forward models and find that the spectrum is fit well by models with moderate cloud lofting or haze formation efficiency, over a wide range of metallicities ($1 - 100\times$ solar). We combine the constraints on the envelope composition with an interior structure model and estimate that the core mass fraction is $\gtrsim0.3$. With a bulk composition reminiscent of that of Neptune and an orbital distance of 0.15 AU, HD 106315 c hints that planets may form out of broadly similar material and arrive at vastly different orbits later in their evolution.


INTRODUCTION
The origins of Uranus and Neptune remain mysterious. Based on current data, it is not known if they formed by core accretion or gravitational instability, whether they grew in their current locations or underwent migration, or how long it took them to form (Atreya et al. 2020, and references therein). One of the challenges in modeling these planets' origin is that their bulk composition is poorly constrained. Uranus and Neptune are so cold that many of the dominant carbon, nitrogen, and oxygen-bearing molecules have condensed out of the observable atmosphere, leaving only methane accessible by remote observation (Helled et al. 2020). There are calls for a space mission to explore one of the ice giants in situ and measure their atmospheric abundances directly with a probe; however, such a mission is over a decade away (Simon et al. 2020).
Meanwhile the search for extrasolar planets has revealed an abundance of Neptune-size worlds (e.g. Coughlin et al. 2016). Many of these have short orbital periods and correspondingly high equilibrium temperatures (up to 2000 K), meaning that major volatile species are expected to be in the gas phase in the observable part of the atmosphere (Moses et al. 2013). Atmosphere characterization of these hotter exo-Neptunes provides an opportunity to determine their chemical compositions, well in advance of in situ measurements of the Solar System ice giants.
Precise near-infrared transmission spectra are available for fewer than a dozen exoplanets in the Neptunemass range, 10 − 40 M ⊕ (Crossfield & Kreidberg 2017;Kreidberg et al. 2018b;Spake et al. 2018;Mansfield et al. 2018;Benneke et al. 2019a,b;Libby-Roberts et al. 2020;Guo et al. 2020;Chachan et al. 2019). Planets of this size are expected to have modest H/He envelopes ( 1% by mass), with a diversity of atmospheric metal enrichment (e.g. Fortney et al. 2013;Wolfgang & Lopez 2015). The transmission spectra measured to date have a wide range of properties that match the diversity expected from theoretical models. Some planets appear to have very high metallicity envelopes (e.g. the ∼ 1000× solar composition inferred for GJ 436b; Morley et al. 2017). Others have lower metallicity, more akin to Jupiter's < 10× solar composition (HAT-P-26b; Wakeford et al. 2017). The planets also have a wide range of cloud and haze properties, from cloud-free to very high altitude condensates, which complicate the interpretation of the measured spectra (Kreidberg et al. 2014;Crossfield & Kreidberg 2017). To fully explore the diversity of the exo-Neptune population and identify cloud-free atmospheres, a larger sample size is needed, which is the goal of the ongoing large HST program GO 15333 (PIs I. Crossfield and L. Kreidberg). In total this program will obtain transmission spectra for five additional Neptune-size exoplanets, including the subject of this work, HD 106315 c.
First observed by the K2 mission Rodriguez et al. 2017), HD 106315 c has a radius of 4.0 ± 0.4 R ⊕ and a mass of 15.2 ± 3.7 R ⊕ (Barros et al. 2017) 1 . The planet has a 21.06 day orbit around its F5-type host star, and an equilibrium temperature of 870 ± 20 K (assuming full heat redistribution and zero Bond albedo). Thanks to the bright host (H magnitude = 8.0), HD 106315 c is one of the most accessible candidates for atmosphere characterization with transmission spectroscopy, with a Transmission Spectroscopy Metric equal to 119 (this metric is a proxy for the expected signal-to-noise of the transmission spectrum; Kempton et al. 2018). Compared to other exo-Neptunes with precise spectra, HD 106315 c has a longer period and a more massive host star (Crossfield & Kreidberg 2017). It is also part of a multi-planet system, with an interior 2.6 R ⊕ planet orbiting the star every 9.6 days. There was also an observation on 15 June 2018 that failed due to lost guiding. Each transit observation consists of time series exposures over six continuous HST orbits. The first exposure of each orbit was a direct im-age with the F126N filter. Subsequent exposures used the G141 grism, which covers the wavelength range from 1.1 − 1.7µm. The exposures used the SPARS25, NSAMP = 8 readout mode which has a total integration time of 138.4 s. The observations used spatial scanning mode, which spreads the light in the cross-dispersion direction during the exposure, enabling longer integration times before saturation (Deming et al. 2013). The scan rate was 0.213"/sec, yielding a total scan height of 31" (238 pixels). We observed 14 exposures per orbit, for an observing efficiency of 73%.
A single transit was also observed by the Spitzer Space Telescope (Werner et al. 2004;Fazio et al. 2004) with the IRAC2 4.5 µm photometric channel on 2017 April 19 -20 as part of Program 13052 (PI: M. Werner). The observations used PCRS peak-up mode, which positions the target precisely on a pixel with minimal sensitivity variations 2 . The observation began with a 30-minute stare to allow the spacecraft to thermally settle, followed by 32168 s (8.9 hours) of science data with an exposure time of 0.4 seconds. Two transits were also observed by K2 (previously described in Crossfield & Kreidberg 2017;Rodriguez et al. 2017).

HST/WFC3
We used a custom data reduction pipeline to process the HST transit observations (described in detail in Kreidberg et al. 2018b). The starting point for our reduction was the ima data product provided by the Space Telescope Science Institute. These images have an intermediate level of processing, with corrections applied for dark current, linearity, and flat fielding. To extract spectra from the images, we used the optimal extraction routine of (Horne 1986). This algorithm iteratively masks bad pixels in the image, and provides a convenient method to reject cosmic rays from spatial scan data. To estimate the background, we identified a region of pixels that was not contaminated by flux from the target or any nearby stars (rows 10−70 and columns 400 − 500) and calculated the median count rate in this region. We subtracted the background and extracted the spectra from each up-the-ramp sample separately, and summed them to produce a final spectrum from the exposure. To account for spectral drift, we interpolated each spectrum to the wavelength scale of the first exposure of the first visit. We generated spectroscopic light curves by binning the spectra into 22 wavelength chan-nels over the wavelength range 1.125 − 1.65 µm. This binning corresponds to roughly five pixels in the spectral direction. The binning is about twice as coarse as the native resolution of the grism, and was chosen to average over variations in sensitivity between individual pixels. Figure 1 shows the band integrated light curve, the background counts, and the spectral shifts for each visit.
We fit the light curves with a joint model of the transit and the instrument systematic trends. In agreement with previous work, we found that the first orbit of every visit and the first exposure in each orbit were strongly affected by a ramp-like systematic (caused by charge traps filling up in the detector; Zhou et al. 2017). This systematic is visible in the raw data, shown in Figure 1. Following past studies, we removed the first orbit of the visit and the first exposure of the remaining orbits in our analysis (e.g. Kreidberg et al. 2014). The trimmed data had three full orbits per visit of out-of-transit baseline, which is sufficient to fit for visit-long trends.
To model the transit signal, we used the batman package (Kreidberg 2015). For the broadband light curve fit, the free parameters for the transit model were the ratio of planet to stellar radius R p /R s , the time of central transit T c , the orbital inclination i, the ratio of semimajor axis to stellar radius a/R s . We fixed the eccentricity to zero. We ran an initial fit with free parameters for a linear limb darkening coefficient, and found excellent agreement with predictions from a Kurucz ATLAS9 stellar model with T eff = 6250K, log g = 4.5 (cgs), and [Fe/H]=−0.2 (Castelli & Kurucz 2003). We therefore fixed the limb darkening on the predicted quadratic coefficients from the model for the remainder of the analysis 3 . For the spectroscopic channels, we fixed T c , a/R s , and i on the best-fit values from the broadband light curve.
To model systematic noise from the instrument, we multiplied the model transit light curve by the analytic model-ramp function, previously used for WFC3 data analysis (Equation 3; Kreidberg et al. 2018b). Briefly, this function fits an exponential ramp to each orbit, and a visit-long trend. For the HD 106315 c data, there was no significant improvement to the light curve fit for a quadratic term in the visit-long trend, so we used a linear term only. The light curves and best-fit models are shown in Figure 2. To test for correlated noise in the light curve, we binned the data in time over a range of bin sizes from 2 − 20 points and calculated the rms for  The open circles in the raw flux correspond to data points that we did not include in our light curve fits due to larger amplitude instrumental ramps. The vertical offset in the top row is due to spatial scanning, which alternates between forward and reverse directions on the detector. The total counts are higher when the detector is read out in the same direction as the spatial scan.
each bin size (see Figure 3). The rms decreases with the square root of the number of points per bin, indicating that the noise is uncorrelated in time.
We used the dynesty package to estimate parameter uncertainties for our model (Speagle 2020). The package uses dynamic nested sampling to evaluate constant likelihood contours over the full prior volume. To ensure that we did not underestimate the uncertainties, we rescaled the per point errors on the data such that the reduced χ 2 of the best fit model was unity. The error bars increased by a median (mean) of 9% (12%). The dynesty runs were halted when the remaining contribution to the evidence was estimated to be below 0.01 of the total. The resulting median and 68% credible intervals for the transit depths are listed in Table 1.

Independent Analysis of the WFC3 Data
We also carried out an independent data reduction and analysis. The data were reduced following the methodology previously described by Evans et al. (2016Evans et al. ( , 2017. Briefly, the spectra are extracted from each ima frame by summing the difference of successive up-the-ramp samples while masking crossdispersion regions away from the target to reject cosmic rays and nearby contaminating sources. A wavelengthindependent background value was subtracted from each spectrum by taking the median pixel value in a 30 × 250 pixel box away from the target. Broadband light curves were produced for each visit by summing each spectrum along the full dispersion axis. The broadband light curves were fit jointly, with the systematics and transit mid-times allowed to vary separately, and R p /R s shared across visits. Other transit parameters such as a/R s and i were fixed to the median values reported in . For the systematics, we adopted the double-exponential ramp treatment described in de Wit et al. (2018) and also allowed the white noise to vary as a free parameter, implemented as an increase above the formal photon noise value.
Following the broadband light curve fit, we produced spectroscopic light curves in 14 channels spanning the 1.122 − 1.642µm wavelength range, following the methodology described in Evans et al. (2016), which is based on an original implementation of Deming et al. (2013). This procedure effectively removes the commonmode component of the systematics in each wavelength channel, which is dominated by the ramp systematic. As such, for our spectroscopic light curve fits, a simple linear time trend and variable white noise level were adequate for modeling the systematics. We also allowed R p /R s to vary, while holding all other transit parameters fixed to the white light curve values. In both the white light curve and spectroscopic light curve fits, a quadratic limb darkening law was adopted with coefficients held fixed to values determined by fitting the limb darkened profiles of an ATLAS9 stellar model with the same parameters listed in 3.1.
The resulting transmission spectrum is compared to the first analysis in Figure 4. The two spectra agree to well within 1σ, and the uncertainties on the transit depths are consistent after accounting for the difference in wavelength bin size. The model-ramp fit has a me-     dian uncertainty on the transit depth of 23 ppm (for 0.025 µm bins), and the common-mode fit has a median uncertainty of 17 ppm (for 0.037 µm bins). Given the good agreement between the two independent analyses, we use the model-ramp results (listed in Table 1) for the remainder of the analysis.

Spitzer
In addition to the HST and K2 data, we also analyzed a single transit of HD 106315 c observed with Spitzer in the 4.5 µm bandpass. We follow a similar approach as the one described in Berardo et al. (2019), which detrends the data using the Pixel Level Decorrelation method outlined in Deming et al. (2015). We first applied a median filter to each pixel in the image and calculated a background level for each frame by taking the median of the flux in an annulus centered on the point spread function. We estimated the centroid of each frame by fitting a two dimensional Gaussian to the  image, and obtained a light curve using a fixed radius aperture.
We modeled systematics in the light curve by weighting the nine brightest pixels individually as well as fitting for a quadratic time ramp. We then chose the combination of pixel coefficients, aperture size, and timeseries binning that resulted in the smallest root mean square (rms) deviation. The optimal aperture radius was found to be 2.4 pixels. We used an MCMC sampler to estimate uncertainties, and fit the systematic model simultaneously with a transit model from batman (Kreidberg 2015). We kept the period, inclination, and distance a/R fixed to the values 21.0564 days, 88.501 • , and 26.769, respectively (based on the HST white light curve fit), and allowed the depth and transit center to vary. We also left the uncertainty of the data points as a free parameter, which we found converged to the rms scatter of the raw light curve. We also held fixed the quadratic limb darkening parameters, which were also estimated from a Kurucz ATLAS9 stellar model. The transit light curve and best fit model are shown in Figure 5, and the fit results are summarized in Table 2.

K2
To provide a broadband, optical-wavelength transit depth for comparison with our infrared observations, we reanalyzed the 30-minute-cadence K2 photometry of HD 106315. Although several analyses of these K2 data have already been published Rodriguez et al. 2017), our reanalysis takes advantage of the tighter constraints on orbital parameters (a/R s and i) provided by the higher-cadence Spitzer and HST observations. Our analysis used largely the same approach described by ), but with a few changes. First, we used a different set of K2 photometry 4 which had a substantially lower rms. Second, we fixed two key orbital parameters to the following values: a/R s = 26.769, and i = 88.501 • . Third, in contrast with the analysis of  we allowed no dilution that would potentially decrease the observed transit depth (and so bias the analysis toward larger R p /R s ). We neglected dilution because high-resolution imaging and spectroscopy show no nearby stars within 5 magnitudes of HD 106315 at distances < 0.1 arcsec (Kosiarek et al., in prep). Finally, we held the quadratic limb darkening parameters fixed to the values predicted by an ATLAS9 stellar model (u 1 = 0.365 and u 2 = 0.244). The transit parameters derived from this analysis are listed in Table 2.

Potential Impact of Star Spots
Unocculted star spots and plages can significantly contaminate exoplanet transmission spectra (Pont et al. 2013;Rackham et al. 2018Rackham et al. , 2019. In general, F-stars like HD 106315 have lower spot covering fractions than stars of later spectral type (Rackham et al. 2019). The K2 light curve for HD 106315 shows variability with amplitude 0.1% over a timescale of 75 days, a typical value for mid-F stars (Rodriguez et al. 2017). This amplitude corresponds to a covering fraction of 0.1 ± 0.1%. The expected amplitude of the stellar contamination spectrum is 0.0001 − 0.0002× the transit depth. For the transit depth of HD 106315 c (1000 ppm), the expected stellar contamination is 0.1 − 0.2 parts per million -a negligible contribution.

petitRADTRANS retrieval analysis
We used the open-source software package petitRADTRANS (pRT) (Mollière et al. 2019), which is a fast spectral synthesis tool for exoplanet atmospheres. We connected pRT to the PyMultiNest tool (Buchner et al. 2014), which is a Python wrapper of the Multi-Nest (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2013 implementation of nested sampling (Skilling 2004).
The atmosphere was modeled with the vertically constant temperature and absorber mass fractions of H 2 O, CH 4 , CO 2 , CO and N 2 as free parameters. We also included the cloud top pressure of a gray cloud deck as a free parameter. The atmospheric mean molecular weight (MMW) was calculated from the parameterized absorber abundances, assuming that the remaining mass is contributed by H 2 and He, with a H 2 :He mass ratio of 3:1. Our full model included the line opacities of H 2 O, CH 4 , CO 2 and CO, as well as the Rayleigh scattering cross-sections of these species, in addition to those of H 2 , He and N 2 . N 2 may thus be thought of as a proxy for (mostly) invisible species in the atmosphere that can increase its MMW. Instead of retrieving a reference radius  at a given pressure we retrieved a reference pressure P 0 at a given radius, where we made sure that the fixed reference radius is chosen at values appropriate for placing the retrieved reference pressure values within the atmospheric pressure domain. We placed log-uniform priors on the absorber mass fractions of H 2 O, CH 4 , CO 2 , CO and N 2 between 10 −10 and 1, requiring that the sum of all mass fractions is below unity. The temperature was allowed to vary between 400 and 1000 K. The cloud and reference pressure could be placed at any location within the atmospheric pressure domain, imposing a loguniform prior. However, we note that the posterior distribution for the water abundance (which is the only species we detect tentatively) is sensitive to the choice of prior bounds, particularly if regions are explored that do not produce any difference in the model spectrum (for example, very deep clouds). We additionally tested retrieving a scattering haze (κ haze = κ 0 [λ/λ 0 ] γ ), a cloud patchiness parameter (mixing clear and cloudy terminators), and the planet's gravity within measurement uncertainties, but none of these tests significantly changed our results.
In order to test how reliably water can be detected in our spectrum we followed the approach introduced in Benneke & Seager (2013). Our full model retrieved the abundance of all absorbers listed above. Then we iteratively removed one absorber at a time and re-ran the retrieval. Comparing the evidences between the full model and the model lacking a given species allows to assess whether the observation is in favor of that species being included in the model. The retrieved atmospheric properties for the full model are shown in Figure 6. In addition to the "full retrieval model", which includes all five molecules (H 2 O, CH 4 , CO 2 , CO and N 2 ), we ran five additional retrieval models, each with one molecular species removed at a time. This approach of removing one molecular species from full model at a time enables us to unambiguously check for each individually molecule and captures any ambiguity that may be introduced by overlapping absorption features (Benneke & Seager 2013). The resulting evidences Z are listed in Table 3. The Bayes factor is the ratio of evidences. Bayes factors greater than 100 are considered decisive, 10 − 100 is strong, 3.2 − 10 is substantial, and below 3.2 is insignificant (Kass & Raftery 1995). None of the tested species is substantially favored to be included in our model. However, the full model is slightly favored when compared to a model that removed H 2 O, with a Bayes factor of 1.7. The fact that water is the only molecule that can lead to noticeable differences in the fit is not surprising, since WFC3 spectra are predominantly sensitive to absorp- -7.5 -2.5 n CO2 -7.5 -2.5 n CO -7.5 -2.5 n CH4 -7.5 -2.5 n N2 -5.88 +3.31 −2.99 Figure 6. The retrieved atmospheric properties for the full retrieval with petitRADTRANS. The panels show the posterior distribution of parameters from the nested sampling run. Darker shading corresponds to higher posterior probability. The diagonal shows a 1-dimensional histogram for each parameter, with dotted lines denoting the median and 1σ credible interval. The molecular abundances are the logarithm (base 10) of their volume mixing ratio. For reference, a solar composition gas at 1 mbar pressure and 800 K has nH 2 O = −3.65, nCO2 = −6.32, nCO = −3.37, nCH 4 = −4.43 and nN 2 = −4.23.
tion from H 2 O. Higher precision measurements of the transmission spectrum are needed to uniquely identify absorbing species in the atmosphere of HD 106315 c.
The retrieved spectrum is shown in Figure 7. The amplitude of spectral features in the best fit model is about 30 ppm (7× smaller than that expected for a solar composition, cloud-free atmosphere). This observed peakto-trough amplitude corresponds to 0.8 ± 0.04 H 2 /He scale heights (assuming µ = 2.3 atomic mass units, T = 800 K, and g = 6.0 m/s 2 ). In general, to produce features of this amplitude, models have either (1) enhanced mean molecular weight (which decreases the atmospheric scale height and shrinks the spectral features), (2) high altitude clouds, which truncates the spectral feature at the cloud deck altitude, or (3)  . The transmission spectrum of HD 106315 c (points with 1σ uncertainties) compared to retrieved spectra from the FULL model (teal shading) with petitRADTRANS. The H2/He atmospheric scale height is indicated on the right y-axis, assuming a solar composition atmosphere at the planet's equilibrium temperature (the true scale height is likely smaller, due to enhanced metallicity and/or lower temperature). The tentative detection of water absorption is driven by the small increase in transit depth near 1.4µm.
of HD 106315 c, the retrieval prefers scenarios (2) and (3), with the highest posterior probability for moderate cloud coverage. A broad range of H 2 O abundances are consistent with the data (3 × 10 −4 − 290× solar at 1σ). The "solar" water abundance corresponds to the chemical equilibrium water volume mixing ratio for a solar composition gas at 1 mbar pressure and 800 K (2.2 × 10 −4 ). The cloud-top pressure is between P cloud = 0.04 − 130 mbar (at 1 σ). There is some degeneracy between n H2O and P cloud , because higher water abundance pushes the photosphere to lower pressures. There is also a tail of probability towards water-rich solutions with deep clouds (below the observable photosphere). Very high water abundances cannot be ruled out (n H2O < 2100× solar at 2σ confidence, n H2O < 4200× solar at 3σ).

SCARLET retrieval analysis
As an independent check of the results from pe-titRADTRANS, we also interpreted the transmission spectrum with the SCARLET atmospheric retrieval framework (e.g., Benneke & Seager 2012Kreidberg et al. 2014;Knutson et al. 2014;Benneke 2015;Benneke et al. 2019a,b;Wong et al. 2020). Employing SCARLET's free molecular composition mode we defined the mole fractions of H 2 O, CH 4 , CO 2 , CO and N 2 as free parameters and assumed a well-mixed atmosphere. The remainder of the atmosphere gas was assumed to be composed of H 2 and He in solar abundance ratio. We included a gray cloud deck using a free parameter describing the cloud top pressure, and an additional free parameter to capture our prior ignorance of the temperature in the photosphere of HD 106315 c near the terminator.
To evaluate the likelihood for a particular set of atmospheric parameters, the SCARLET forward model in free molecular composition mode computes the hydrostatic equilibrium and line-by-line radiative transfer. We consider the latest line opacities of H 2 O, CO, and CO 2 from HiTemp (Rothman et al. 1998) and CH 4 from ExoMol (Tennyson et al. 2016), as well as the collisioninduced absorption of H 2 and He. We employed loguniform priors between 10 −10 and 10 −0.5 = 31%, but required that the sum of all mass fractions is below unity. We employed log-uniform priors for the cloud top pressures 10 −3 and 10 7 Pa. We used a uniform prior on the photospheric temperature between 620K and 1150K (70-130% of the equilibrium temperature).
SCARLET then determined the posterior constraints by combining the SCARLET atmospheric forward model with nested sampling (Skilling 2004). We ran the analyses well beyond formal convergence to obtain smooth posterior distribution and capture the contours of the wide parameter space in agreement with the transmission spectrum of HD 106315 c. As in Section 4.1, we evaluated the presence of individual molecular species in the atmosphere of HD 106315 c following the strategy outlined in (Benneke & Seager 2013).
The retrieval results are shown in Figure 8. Our analysis reveals that a Bayes factor of 2.6 in favor of the pres- ence of water vapor in the atmosphere of HD 106315 c, which can be regarded as tentative evidence. No other molecular species is favored by the data. We also test for the presence of clouds by comparing the full retrieval model to a model that lacks clouds in the hypothesis space, but find no evidence in the observational data. The small differences between the evidence computed with SCARLET versus petitRADTRANS can be at-tributed in the difference in prior volume for the two analyses. We perform the final parameter estimation using the full retrieval model including the five molecules (H 2 O, CH 4 , CO 2 , CO and N 2 ) and gray clouds. The best-fitting model matches all data points within their 1-σ uncertainties. A wide range of models is consistent with the data, in agreement with the results from peti-tRADTRANS.

CLOUD AND HAZE MODELS
The retrieval analysis from the previous section showed that the muted water feature in the transmission spectrum is consistent with a low metallicity composition with high altitude condensates. To explore what condensate properties are plausible for HD 106315 c, we ran forward models with physically motivated cloud and haze opacity.
We include the condensation of Na 2 S, KCl, and ZnS, which are expected to condense at the temperature of HD 106315 c (T eq = 800 Kelvin). We calculate cloud altitude and height along the cloud-free pressuretemperature profile; the cloud properties are calculated using the methods described in Ackerman & Marley (2001); Morley et al. (2012) for each metallicity, assuming a range of sedimentation efficiencies (f sed =2, 1, 0.5, and 0.1), a parameter which controls the cloud particle size and cloud height. This model calculates the cloud optical depth, single-scattering albedo, and asymmetric parameter for each layer of the atmosphere. Example pressure-temperature profiles with cloud condensation curves are shown in Figure 9.
To calculate transmission spectra, we use the transmission spectrum model described in the appendix of Morley et al. (2017). Gas opacity from H 2 collisionally induced absorption, CO 2 , H 2 O, CH 4 , CO, NH 3 , PH 3 , H 2 S, Na, K, TiO, VO, and HCN is included. We calculate models for each metallicity and f sed combination considered. Figure 10 shows the goodness of fit for the cloudy model grid compared to the WFC3 transmission spectrum. The K2 and Spitzer data are not precise enough to significantly affect the goodness of fit. The best fits have small water absorption features with amplitude of around 30 ppm. The amplitude of features in the models is a trade-off between metallicity and cloud altitude: higher metallicity models tend to have a smaller scale 400 600 800 1000 1200 1400 1600 1800 temperature (K) 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 pressure (bar) M n S N a2 S Z n S C r K C l solar 10 × solar 100 × solar 300 × solar Figure 9. Pressure-temperature profiles (solid lines) for HD 106315 c compared to condensation curves for expected cloud species (dashed lines). The models assume are in radiative-convective and chemical equilibrium. The condensation curves are calculated for a 100× solar metallicity composition; for lower metallicities, the condensation curves shift left (by approximately 100 K per 1 dex metallicity). The shaded region marks the range of pressures sensed by transmission spectroscopy, assuming 100× solar metallicity.
height and thus smaller features. Lower sedimentation efficiency also decreases the feature amplitude. Small f sed values loft cloud particles higher in the atmosphere, obscuring spectral features. As shown in Figure 9, the cloud base is typically below the pressure level sensed by the observations, so f sed values 0.5 are required to loft the cloud into the observable atmosphere. For the HD 106315 c spectrum, the best fit models are high metallicity atmospheres (100 − 300× solar), or lower metallicity with high-altitude clouds (f sed < 0.5).

Haze Models
We also calculate transmission spectra for hazy atmospheres using the photochemistry, microphysics, and transmission spectrum models of Kawashima & Ikoma (2018) in the same way as Kawashima et al. (2019) and Kawashima & Ikoma (2019). We first perform photochemical simulations to derive the steady-state distribution of gaseous species. Then, since haze monomer production rate is uncertain for exoplanets, we assume a certain fraction of the sum of the photodissociation rates of the major hydrocarbons in our photochemical model, CH 4 , HCN, and C 2 H 2 , would result in haze monomer production. We call this fraction as haze formation efficiency f haze following Morley et al. (2013). With this assumption, we derive the steady-state distribution of haze particles by microphysical simulations. Finally, we  Figure 10. A sample of cloud forward models compared to the observed WFC3 spectrum (left) and a goodness-of-fit for the full grid of cloudy models (right). The grid cell shading indicates the reduced χ 2 of the fit to the WFC3 data. The fit has 18 degrees of freedom (21 data points, free parameters for metallicity, f sed , and reference radius).
model transmission spectra of the atmospheres with the obtained profiles of haze particles and gaseous species. For the temperature-pressure profile, we use the online-available 5 analytical model of Parmentier & Guillot (2014) assuming an internal temperature of 100 K and their correction factor of 0.25, which corresponds to the case where the irradiation is efficiently redistributed over the entire planetary surface. For the other input parameters, we use the default opacities (Valencia et al. 2013;Parmentier et al. 2015) and Bond albedo. We include convection. Solar elemental abundance ratios are taken from Lodders (2003). For the UV spectrum of HD 106315, we use that of the Sun from Segura et al. (2003) because of its similar stellar type (F5, Houk & Swift 1999). We assume a constant eddy diffusion coefficient of 10 7 cm 2 s −1 throughout the atmosphere for both photochemistry and microphysics calculations. We assume a monomer radius of 1 nm and an internal density of haze particles of 1 g cm −3 . The refractive index of haze is uncertain for exoplanets, so we consider two representative cases, tholin (Khare et al. 1984) and soot (Hess et al. 1998).
For the calculation of transmission spectra, we treat the reference radius at 10 bar pressure level as a parameter. We find the appropriate value with a grid of 0.1% of the observed transit radius which yields the minimum reduced χ 2 with 18 degrees of freedom (21 data points minus 3 free parameters of metallicity, f haze , and reference radius), for each case. We account for the transmission curve of the WFC3 G141 grism from the SVO Filter Profile Service 6 (Rodrigo et al. 2012;Rodrigo & Solano 2013).
The left panel of Figure 11 shows several representative models compared to the measured WFC3 spectrum; models for clear atmospheres of three different metallicities, as well as hazy (tholin) atmospheres with haze formation efficiency tuned to fit the WFC3 data well. The error bars for the K2 and Spitzer points are large, and therefore have a negligible effect on the goodness of fit. The right panel of Figure 11 shows the goodness of fit for the model grids. We find that modest haze formation efficiencies of 10 −3 − 10 −4 generally match the observed spectra for all the three different metallicities, for both tholins and soots. This is because smaller scale height due to increasing metallicity can be compensated out by smaller fiducial monomer production rate. Overall, these haze production efficiencies are orders of magnitude lower than the extreme values required to  Figure 11. (Left) Representative haze forward models compared to the measured WFC3 spectrum; models for the clear atmospheres of 1 (dark green line), 10 (orange line), and 100 (purple line) × Solar metallicity and those for the hazy (tholin) atmospheres of 1 (pink line), 10 (light green line), and 100 (yellow line) × Solar metallicity. The haze formation efficiency that fits the observed data well is chosen. (Right) Goodness-of-fit for the full grid of tholin and soot models. The grid cell shading indicates the reduced χ 2 of the fit to the WFC3 data. The fit has 18 degrees of freedom (21 data points, free parameters for metallicity, f haze , and reference radius).
reproduce the featureless spectrum of GJ 1214b (Morley et al. 2015;Kawashima et al. 2019). As noted above, the NUV irradiation of HD 106315 c is likely higher than that of GJ 1214b, so more haze precursors are present and lower haze production efficiency is needed.

INTERIOR STRUCTURE MODELS
Comparison between interior structure and envelope metallicity can provide additional constraints on the bulk composition of the planet (Kreidberg et al. 2018b;Thorngren & Fortney 2019). For example, given knowledge of the envelope metallicity, it is possible to put limits on the core mass, that otherwise suffers from large degeneracy for planets in the 2−5 R ⊕ radius range (Rogers & Seager 2010). We evaluate the internal structure of HD 106315 c with a model consisting of an inner core and a H/He outer envelope enriched with some various amounts of water and rock (in a 50-50 ratio), using the methods described by (Thorngren et al. 2016). We explored two limiting cases for the core composition: one is composed entirely of isothermal rock with radioactive heating, and the other is composed of convective water. Using the observed mass (with error), radius (with error), age (with error), and flux (ignoring error), we retrieved the core mass over a range of envelope metallicities. Our results are shown in Figure 12.
In the absence of any information about the envelope composition, the core mass fraction for HD 106315 c could range anywhere from 0 to 1. The higher the envelope metallicity, the lower the core mass fraction required to explain the observed mass, radius, and age of the planet. To help break this degeneracy, we compared the results from the atmospheric retrieval with Posterior Probability Figure 12. Core fraction versus envelope metallicity from interior structure modelling for a rocky core (red line with 1σ uncertainty shaded) and a water core (blue line with 1σ uncertainty). The retrieval results for the envelope metallicity are over-plotted as a histogram, with the 1σ credible interval indicated by the blue shaded region. the interior structure model (using water abundance as a proxy for envelope metallicity). The retrieval results are shown alongside the interior structure model in Figure 12. Using the retrieved abundance of H 2 O from petitRADTRANS as a proxy for the envelope metallicity (5 × 10 −4 − 290× solar at 1σ confidence), we estimate that the core mass fraction is greater than 30% regardless of the core composition (rock or ice). These conclusions generally resemble our understanding of the bulk composition of Uranus and Neptune, which are expected to have a core mass of 80 -90% (Hubbard et al. 1991;Fortney & Nettelmann 2010;Nettelmann et al. 2013). Follow-up atmosphere characterization with higher precision and broader wavelength coverage can further constrain the envelope metallicity and core mass of HD 106315 c.

DISCUSSION AND CONCLUSIONS
The number of small exoplanets with precise transmission spectra is growing, and already the population shows diversity in atmospheric properties. Some appear to have envelope metallicities below that of Neptune (e.g. HAT-P-26b; Wakeford et al. 2017), whereas others require higher metallicity (GJ 436b; Morley et al. 2017). Some planets are blanketed with thick high altitude clouds or haze (particularly GJ 1214b; Kreidberg et al. 2014), while others have deeper condensates or even cloud-free atmospheres (Benneke et al. 2019b;Madhusudhan et al. 2020). This diversity is expected from theoretical models. For example, planet population synthesis predicts a wide range of envelope enrichment for sub-Neptunes (e.g Fortney et al. 2013). Similarly, cloud and haze models indicate that condensate properties may vary widely across different planets. Condensate formation depends on many different atmospheric properties (e.g. temperature, metallicity, UV irradiation and vertical mixing) so there is no one-size-fits-all model to predict whether an atmosphere is cloudy or clear at the pressure levels sensed by transmission spectroscopy (Morley et al. 2015;Gao & Benneke 2018;He et al. 2018;Hörst et al. 2018;Kawashima et al. 2019;Ohno et al. 2020).
Where does HD 106315 c fit into this diverse picture? The small amplitude of spectral features is consistent with other sub-Neptunes, which all have feature amplitudes attenuated relative to expectations for solar composition, cloud-free atmospheres (Crossfield & Kreidberg 2017). Intriguingly, the amplitude of spectral features (0.8 ± 0.0.04 H 2 /He scale heights) agrees well with the demographic trend noted in Crossfield & Kreidberg (2017); Libby-Roberts et al. (2020), that shows an increase in the amplitude of the WFC3 water feature with planet equilibrium temperature. This is a somewhat surprising finding, because there are many factors (noted above) that affect the observed spectral feature amplitude for planets in this population. A further demographic study of water absorption in sub-Neptunes will be explored in a follow-up paper.
The tentative water detection for HD 106315 c is consistent with a wide range of abundances (3×10 −4 −290× solar at 1σ confidence), and is most comparable to that estimated for HAT-P-11b and K2-18b (Fraine et al. 2014;Benneke et al. 2019b;Chachan et al. 2019). Low metallicities (< 50× solar), akin to those GJ 3470b, HAT-P-26b, and WASP-107b (Benneke et al. 2019a;Wakeford et al. 2017;Kreidberg et al. 2018a) are possi-ble for HD 106315 c, provided it has some condensates in its atmosphere that truncate the amplitude of the water feature. The condensate properties are modest relative to extremes like GJ 1214b, which has a featureless spectrum requiring either very low sedimentation efficiency clouds and high atmospheric metallicity (f sed ≤ 0.1 and 1000× solar composition), or very efficient photochemical haze production ( 10% efficiency) for a 50× solar metallicity composition (Kreidberg et al. 2014;Morley et al. 2015). For comparison, the transmission spectrum of HD 106315 c is fit well with f sed ≤ 0.5 or haze production efficiencies of 10 −3 − 10 −4 .
The tentative detection of a small water feature in HD 106315 c makes it an intriguing candidate for followup observations to further characterize its atmosphere. Infrared observations are a particularly promising avenue -spectroscopy in the 4−5µm range is sensitive to absorption from CO 2 , a prominent feature expected in high metallicity atmospheres (Moses et al. 2013). If the atmosphere has lower metallicity but is cloudy/hazy, infrared observations are also promising because the condensates may have lower opacity at longer wavelengths (e.g GJ 3470b; Benneke et al. 2019b). Future transmission spectroscopy observations with JWST could potentially distinguish between these possibilities (Greene et al. 2016), and confirm whether HD 106315 c does indeed have a Neptune-like core mass and envelope composition. If it does, that will provide new incentive for formation models to produce ice giants on a wide range of orbits from 0.15 AU (that of HD 106315 c) to 30 AU (that of Neptune).