Correcting Exoplanet Transmission Spectra for Stellar Activity with an Optimized Retrieval Framework

To determine how extreme the contamination effects are for each spot case considered and explore how well these effects can be mitigated, we must first produce a synthetic star–planet system, for which the stellar and planetary parameters are known a priori. We consider a synthetic K dwarf host star characterized by the following parameters (T_eff = 4750 K, R_⋆ = 0.8 R_⊙, M_⋆ = 0.8 M_⊙) and displaying activity in the form of a single-unocculted spot, which we model using the methodology described in Section 2.3. This is the same methodology outlined in Cracchiolo et al. (2021b) with some subsequent, minor improvements to its efficiency. A K dwarf was chosen as stars of later spectral types are typically more likely to be active (Berdyugina 2005; Ciardi et al. 2011; Hartman et al. 2011; McQuillan et al. 2014; Rackham et al. 2018, 2019). We decided against using an M dwarf host for this preliminary study as this is the region of the main sequence in which stars transition to a fully convective regime (Reiners & Basri 2010), as such, stellar activity may manifest differently on these stars. We choose to focus on single-spot cases for this initial benchmarking study for several reasons. First, this represents the simplest spot-contaminated disk model that can be considered, which makes it invaluable as a baseline for future investigations. Encompassed with this is that it reduces the dimensionality required for the input model. Having multiple spots present requires the location of each of them on the stellar disk to be defined in order to accurately compute the interplay between their properties and the limb-darkening effect. Second, a single, larger spot also allows us to probe the extremes of this interplay in a way that multiple spots will not as, for a given filling factor, the active region is concentrated at a single latitude rather than dispersed over multiple locations on the stellar disk. Despite this, the extension to multiple-spot cases is comparatively straightforward as described in Cracchiolo et al. (2021b). The results of several preliminary multiple-spot cases are presented in Section 4.5 for completeness.

The transiting planet is a temperate sub-Neptune (R_P = 3 R_⊕ (0.273 R_Jup), M_P = 5 M_⊕ (0.0157 M_Jup), and T_P = 400 K) with a primary atmosphere containing water ($\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$ = −3) and H and He present as fill gases with a ratio of 0.172. Rayleigh scattering and collision-induced absorption (CIA) are also included and introduce wavelength-dependent contributions to the opacity. A smaller planet with a primary atmosphere is used as its scale height should result in detectable absorption features but the smaller S/N of such features means that they are more susceptible to being masked by stellar contamination. As such, the accuracy of our correction method, and any biases that may be inadvertently introduced, are comparatively far more important here than when considering an atmosphere with a much higher S/N, for example, that of a hot Jupiter. The orbital inclination is set to 88° so that the effects of a central-unocculted spot (ϕ_spot = 0°) can also be explored.

Using the stellar and planetary parameters described above and the retrieval code TauREx3 (Al-Refaie et al. 2021, 2022), we produce a forward model that is equivalent to the idealized, uncontaminated transmission spectrum that would be observed in the presence of a completely homogeneous host star. This synthetic spectrum has a wavelength coverage of 0.5–9.5 μm and a resolution of 200. This wavelength range has been chosen as it is similar to the regions that are/will be covered by the JWST and Ariel instruments but with greater coverage and resolution in the optical regime where the effects of stellar contamination are strongest. Error bars of 10 ppm are introduced for all data points regardless of wavelength. This noise level was chosen as 10 ppm is broadly consistent with the most precise observations obtained with the Hubble Space Telescope (HST) instruments (e.g., Edwards et al. 2023), albeit for larger planets further into the IR. We acknowledge that obtaining this level of precision, in reality, would be extremely challenging and heavily reliant on an accurate characterization of the host star to mitigate the astrophysical noise. We rationalize our choice of using small error bars as our goal at this stage is to be limited by the model assumptions/limitations and not the instrument performance, which will be a focus of future work. Testing our correction methods on an idealized case first allows us to quantify how effective they are and ensure that they do not introduce any unknown, intrinsic bias before using them with real observations.

Throughout this study, we consider 27 spot-contaminated cases (and an uncontaminated case as a reference frame) for the star–planet system outlined in the previous section. From this grid, we investigate the potential biases that an unocculted spot may introduce and how these may vary as a function of the three spot parameters considered (ΔT_spot, R_spot, and ϕ_spot). The contaminated stellar disks are produced with StARPA (Section 2.3) and the following values are considered for each spot parameter: ΔT_spot = [250, 500, 1000] K, R_spot = [0.1, 0.2, 0.4] R_⋆, and ϕ_spot = [0, 30, 60]°. This results in 27 unique parameter combinations in which only one spot parameter is varied at a time. The full set of spot parameters used in each case are given in Table 1. Visual representations of each case are shown in Figure 1.

Zoom In Zoom Out Reset image size

From the uncontaminated planetary transit spectrum produced with TauREx, the grid of 27 spot-contaminated spectra described above is constructed using the Cracchiolo et al. (2021b) model with some subsequent improvements made to the implementation. The key point of this methodology with respect to this study is that, through accounting for the limb-darkening effect and defining the exact position of the spot on the stellar disk, any interactions between contamination and limb darkening can be explored. As in Cracchiolo et al. (2021b), the spots are modeled as being circular but with an elliptical 2D projection on the visible stellar disk at nonzero latitudes. The out-of-transit stellar flux from the spot-contaminated star is computed using quadrature integration; the star is divided into 1000 equally spaced annuli and the fractions of each of these annuli covered by the spot are calculated. This enables us to compute both the absolute and normalized intensity profiles (Figure 2).

Zoom In Zoom Out Reset image size

To create the contaminated stellar disks we model the contribution of the quiet photosphere and the spot using the BT-Settl (Allard et al. 2012; Baraffe et al. 2015) stellar spectral model grids from the PHOENIX library (Husser et al. 2013). The spectral emission densities (SEDs) corresponding to the photosphere and the spot are governed by three fundamental stellar parameters: the stellar effective temperature (T_eff), the stellar metallicity [M/H], and the stellar surface gravity ($\mathrm{log}g$). For the purposes of this study, the metallicity and gravity are fixed at [M/H] = 0 (solar metallicity) and $\mathrm{log}g$ = 4.5, respectively, in order to isolate the effects of active regions with contrasting temperatures. These are reasonable assumptions for low-resolution spectroscopy but may need to be reconsidered at higher resolution. In particular, the treatment of $\mathrm{log}g$ may need improvement, as it has been suggested that spots may be characterized by a lower $\mathrm{log}g$ than that of the photosphere due to the localized increase in magnetic pressure (e.g., Solanki 2003). In total, we require four SEDs, one corresponding to the photospheric temperature (T_phot = 4750 K) and three corresponding to the spot temperatures considered in this work: 3750, 4250, and 4500 K, equivalent to a ΔT_spot of 1000, 500, and 250 K, respectively.

Limb darkening is a well-known phenomenon acting to reduce the flux originating from the limbs of the stellar disk with respect to its center (Claret 2000; Howarth 2011). It also varies as a function of wavelength, and as such, different bands are characterized by different intensity profiles with the strongest effects seen in the optical. To account for the limb-darkening effect within our forward model we use the ExoTETHyS package, specifically the SAIL and BOATS subpackages (Morello et al. 2020a, 2020b, 2021), to calculate the limb-darkening coefficients (LDCs) for the star using the PHOENIX_2012_13 database (Claret et al. 2012, 2013) of BT-Settl models and the Claret four-coefficient law (Claret 2000):

$$\begin{eqnarray}&&\displaystyle \frac{{I}_{\lambda }(\mu )}{{I}_{\lambda }(1)}=1-\displaystyle \sum _{n=1}^{4}{a}_{n}{,}_{\lambda }(1-{\mu }^{n/2}),\end{eqnarray} \tag{ 1 }$$

where λ is the wavelength/bandpass being considered, μ = cosθ (where θ is the angle between the line of sight and the normal at the stellar surface), I_λ (μ) is the stellar intensity profile, I_λ (1) is the intensity at the disk center (i.e., where μ = 1), and α_n,λ are the LDCs.

We model the spot using the same LDCs that have been calculated for the star which is a reasonable approximation at first order, but again, may be revised for higher-resolution observations. The absolute fluxes of the star and the spot are also calculated using ExoTETHyS from which the normalized intensity profiles of the spotted star can be calculated as a function of wavelength (Figure 2). The emission from the spot-contaminated stellar disk is calculated as in Equation (2), where S_star,λ , S_phot,λ and S_spot,λ are the spectra of the average star, the quiescent photosphere, and the spot, respectively, for a given wavelength (λ) and F_spot is the spot-filling factor. As such the resulting spectrum for the active star is essentially a combination of the photosphere and spot SEDs weighted by their relative covering fractions (Figure 3). The limb-darkening effect, which has already been defined for the 1000 annuli considered using Equation (1), is also accounted for in this stage.

$$\begin{eqnarray}&&{S}_{\mathrm{star},\lambda }=((1-{F}_{\mathrm{spot}})\times {S}_{\mathrm{phot},\lambda })+({F}_{\mathrm{spot}}\times {S}_{\mathrm{spot},\lambda }).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The resulting chromatic contamination can be described as acting as a contamination factor (ε) relative to the nominal transit depth i.e., the uncontaminated spectrum that would be observed in the case of an inactive star (Equation (3)). Where contamination is present only in the form of lower-temperature spots, as in this study, ε_λ > 1 resulting in increased transit depths at all wavelengths.

$$\begin{eqnarray}&&{\varepsilon }_{\lambda }=\frac{1}{1-{F}_{\mathrm{spot}}\left(1-\frac{{S}_{\lambda ,\mathrm{spot}}}{{S}_{\lambda ,\mathrm{phot}}}\right)}.\end{eqnarray} \tag{ 3 }$$

From the constructed, spot-contaminated stellar disks, we then use the open-source, light-curve analysis package pylightcurve (Tsiaras et al. 2016) to produce the spot-contaminated light curves (Figure 4) for all 200 wavelength intervals. These light curves are then used to construct the contaminated transmission spectrum for each spot case. We highlight that for this preliminary investigation, as we are only aiming to construct the contaminated transmission spectra for use as inputs in the retrievals, we make the assumption that the orbital parameters used in pylightcurve are well constrained/known a priori and we do not introduce any Gaussian noise at the light-curve fitting stage to avoid introducing any additional bias that may arise arbitrarily. As in Section 2.1, error bars of 10 ppm are assumed for the input spectra. A comparison of the uncontaminated transmission spectrum and the spectrum with the highest degree of stellar contamination considered in this study is shown in Figure 5. The strongly wavelength-dependent nature of the spot contamination is evident, with it imparting a steep, positive blueward slope on the spectrum at wavelengths shortward of ∼2 μm. Longward of 2 μm, the contamination acts to introduce an almost constant positive offset equivalent to an overestimate in the transit depth on the order of 5% (>50 ppm) for the worst-case scenario.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Retrievals on the grid of spot-contaminated spectra produced with the StARPA model are conducted using the fully Bayesian atmospheric retrieval code TauREx3 (Al-Refaie et al. 2021). The main benefit of using TauREx for this study is that its modular design allows for a large amount of freedom and flexibility in testing and introducing new models. Mitigating for potential contamination due to stellar activity has been conducted in previous studies in the form of combined stellar-planetary retrievals (e.g., Pinhas et al. 2018; Bixel et al. 2019; Espinoza et al. 2019), and hopefully, this consideration of the host star will progressively become the new norm within the exoplanetary retrieval community. While this is not the first study focused on retrieval biases as a function of the degree of stellar heterogeneity, parameterized by temperature contrast and spatial extent (Iyer & Line 2020), this is the first time that that bias has been explored as a function of three spot parameters, with the inclusion of their position on the stellar disk, which governs how they interact with the limb-darkening effect.

The ASteRA plugin introduces a new, heterogeneous star class to TauREx, which follows a formalism similar to that used in previous stellar activity studies (e.g., Rackham et al. 2018, 2019). Instead of modeling the host star as being homogeneous and characterized by a single temperature/SED, the active star is modeled as a combination of multiple distinct temperature components, in this case, two: the quiet photosphere and the spot. The spot and photosphere SEDs are homogeneous disk-integrated models as they are not spatially resolved. The observed disk-integrated stellar spectrum is then produced by combining the two SEDs, weighted by the filling factor of each component, as shown in Figure 3 and described in Equation (2).

Using the heterogeneous star model requires the addition of two extra fitting parameters to a normal retrieval; the spot temperature T_spot and its filling factor F_spot. For the purposes of this study, we restrict the active regions considered to spots only; however, ASteRA is easily extended to incorporate faculae, as has been done in the analysis of HST STIS observations of the hot Jupiter WASP-17b (Saba et al. 2022).

The core difference between the stellar models of StARPA and ASteRA is the treatment of limb darkening, or lack thereof. In contrast to the forward model produced with StARPA, the ASteRA stellar model is simpler in that the interplay between the spot and the limb-darkening effect is neglected, making it similar to the initial model used in Cracchiolo et al. (2021a). Conducting retrievals without the inclusion of limb darkening helps us gauge how important its inclusion is for low-resolution spectroscopy, particularly with respect to active host stars. With limb darkening neglected, the relative position of the spot on the stellar disk becomes unimportant provided that the entire spot is unocculted, reducing the dimensionality of the model by one. ASteRA requires the fundamental stellar parameters, i.e., effective temperature, metallicity, and $\mathrm{log}g$ as inputs in order to select the correct, corresponding PHOENIX spectra and these are fixed within the retrieval, as are the orbital parameters. This is done with the assumption that for real observations, these parameters will be reasonably accurately and homogeneously constrained a priori. Indeed, this is a high-priority and ongoing effort within the Ariel consortium regarding the stellar parameters (Danielski et al. 2022; Magrini et al. 2022) and the ExoClock project for the orbital parameters (Kokori et al. 2022).

For each spot case, two separate retrievals were conducted: one accounting for activity by fitting for T_spot and F_spot and one where the activity is not accounted for, despite being present. The planetary parameters R_P , T_P , and $\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$ are fit for in all retrievals and the prior bounds set for each parameter are given in Table 2. All retrievals conducted with TauREx use MultiNest (Feroz et al. 2009; Buchner et al. 2014) as the sampler with 500 live points and an evidence tolerance of 0.5. The retrievals with and without activity have a dimensionality of five and three, respectively. Conducting two retrievals allows for a comparison of the Bayesian evidence to obtain the Bayes factor, which quantifies how strongly the model with activity is favored over the one without. This provides an additional, statistically motivated way of quantifying how strong the spot contamination effects are. The retrieved values are then compared to the input/ground truth values for each parameter. From this, we determine how accurately the stellar and planetary parameters can be retrieved using ASteRA and if any bias is introduced as a result. The results of these retrievals are given in Section 3.

Table 2. Fitting Parameters Used within the Retrievals with ASteRA and TauREx and Their Prior Bounds

Fitted Parameter	Prior Bounds	Scale
RP (RJup)	0.75RP ; 1.5RP	Linear
TP (K)	100 ; 1000	Linear
Tspot (K)	3000 ; 4700	Linear
Fspot	0 ; 0.99	Linear
$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$	−12 ; −1	log10

Download table as:  ASCIITypeset image

In this section, we present the results of the retrievals for the 27 spot cases outlined in Section 2.2. To account for the spot contamination in retrieval, the spot parameters T_spot and F_spot are fitted for using the ASteRA plugin. In contrast, in Section 3.2, retrievals are conducted on the same contaminated spectra but with the incorrect assumption that the star is homogeneous. The same retrievals are also conducted on the uncontaminated transmission spectrum, which we term "Case 0". Case 0 acts as both a frame of reference for the highest precision and accuracy obtainable with TauREx for the simulated spectra used in this study, and as a verification that the use of ASteRA does not introduce any intrinsic bias. The results of these retrievals are given in Table 3 and a visual representation of the retrieval accuracy with respect to the planetary parameters is given in Figure 7. On inspection of the retrieved planetary parameters alone, the outlook is overall very positive, with the retrieved values falling very close to the ground truth in almost all of the cases considered. Intuitively, it becomes apparent that the largest errors in the planetary parameters are obtained for the cases where the spectra are most strongly contaminated. These cases are characterized by the largest, highest contrast spots considered in this study (e.g., Cases 7, 8, and 9). Despite showing the largest errors in this study, we emphasize that these errors are not substantial enough to result in a large misinterpretation of the planet. The retrieved parameters are also substantially more accurate than those obtained when the stellar contamination is neglected in the retrieval, the results of which are presented in the following section (Section 3.2). One other thing that becomes particularly evident here is that ASteRA struggles to constrain the spot-filling factor for small spots at high latitudes (Cases 3, 12, and 21) and for small, low-contrast spots (Cases 19, 20, and 21). For these scenarios, the retrieved F_spot is highly degenerate as the comparatively low levels of contamination can be reproduced by a larger number of spot configurations.

Table 3. Retrieved Planetary and Spot Parameters for Each of the 27 Spot Cases Investigated

Case No.	RP (RJup)	T (K)	$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$	Input Tspot (K)	Retrieved Tspot (K)	Input Fspot	Retrieved Fspot
GT	0.2730	400	−3.00	⋯	⋯	⋯	⋯
0	${0.2731}_{-0.0002}^{+0.0003}$	${398.31}_{-5.62}^{+4.91}$	−2.99 ± 0.05	N/A	${4670.15}_{-19.95}^{+20.12}$	0	${0.391}_{-0.378}^{+0.400}$
1	${0.2729}_{-0.0005}^{+0.0004}$	${399.74}_{-6.36}^{+5.74}$	−3.01 ± 0.05	3750	${3918.02}_{-512.02}^{+442.34}$	0.010	${0.018}_{-0.004}^{+0.012}$
2	${0.2732}_{-0.0005}^{+0.0004}$	${395.98}_{-5.92}^{+5.58}$	−2.99 ± 0.05	3750	${4017.22}_{-653.88}^{+490.05}$	0.009	${0.015}_{-0.005}^{+0.030}$
3	0.2736 ± 0.0003	${393.65}_{-6.00}^{+5.22}$	−2.99 ± 0.05	3750	${4663.47}_{-585.54}^{+24.25}$	0.005	${0.305}_{-0.299}^{+0.456}$
4	0.2723 ± 0.0005	${403.97}_{-5.48}^{+5.68}$	−3.05 ± 0.05	3750	${3779.57}_{-110.55}^{+112.28}$	0.040	${0.059}_{-0.003}^{+0.004}$
5	${0.2728}_{-0.0005}^{+0.0004}$	${398.95}_{-5.94}^{+5.62}$	−3.01 ± 0.05	3750	${3760.82}_{-139.26}^{+140.09}$	0.035	${0.045}_{-0.003}^{+0.004}$
6	0.2734 ± 0.0005	${392.42}_{-6.34}^{+6.08}$	−2.96 ± 0.05	3750	${3773.34}_{-384.27}^{+449.70}$	0.020	${0.018}_{-0.004}^{+0.007}$
7	${0.2698}_{-0.0005}^{+0.0004}$	${429.02}_{-5.45}^{+5.98}$	−3.22 ± 0.05	3750	${3787.50}_{-59.60}^{+41.23}$	0.160	${0.234}_{-0.006}^{+0.002}$
8	${0.2711}_{-0.0004}^{+0.0005}$	${411.12}_{-4.94}^{+5.23}$	−3.09 ± 0.05	3750	${3707.40}_{-39.19}^{+39.81}$	0.139	${0.176}_{-0.002}^{+0.003}$
9	0.2737 ± 0.0004	${385.34}_{-6.09}^{+5.98}$	−2.88 ± 0.05	3750	${3689.36}_{-85.94}^{+77.58}$	0.080	0.070 ± 0.003
10	${0.2731}_{-0.0005}^{+0.0004}$	${396.84}_{-5.72}^{+6.09}$	−3.01 ± 0.05	4250	${4301.79}_{-744.13}^{+308.63}$	0.010	${0.013}_{-0.006}^{+0.031}$
11	${0.2733}_{-0.0005}^{+0.0004}$	${396.01}_{-6.14}^{+5.58}$	−3.01 ± 0.05	4250	${4534.88}_{-893.32}^{+139.69}$	0.009	${0.018}_{-0.012}^{+0.491}$
12	0.2733 ± 0.0003	${395.37}_{-5.64}^{+5.45}$	−3.00 ± 0.05	4250	${4669.58}_{-48.39}^{+20.84}$	0.005	${0.401}_{-0.387}^{+0.401}$
13	0.2730 ± 0.0004	${398.71}_{-6.01}^{+5.72}$	−3.02 ± 0.05	4250	${4517.02}_{-248.64}^{+97.03}$	0.040	${0.100}_{-0.049}^{+0.092}$
14	0.2731 ± 0.0004	${397.32}_{-5.79}^{+5.64}$	−3.00 ± 0.05	4250	${4479.76}_{-277.08}^{+128.87}$	0.035	${0.070}_{-0.035}^{+0.076}$
15	0.2733 ± 0.0004	${394.98}_{-5.52}^{+5.48}$	−2.98 ± 0.05	4250	${4299.29}_{-748.47}^{+288.83}$	0.020	${0.015}_{-0.007}^{+0.031}$
16	${0.2716}_{-0.0005}^{+0.0006}$	${413.66}_{-7.66}^{+6.10}$	−3.13${}_{-0.05}^{+0.06}$	4250	${4311.13}_{-46.12}^{+72.60}$	0.160	${0.221}_{-0.005}^{+0.044}$
17	${0.2724}_{-0.0005}^{+0.0007}$	${404.22}_{-8.11}^{+6.42}$	−3.05 ± 0.06	4250	${4310.73}_{-63.42}^{+97.73}$	0.139	${0.171}_{-0.020}^{+0.037}$
18	${0.2736}_{-0.0005}^{+0.0004}$	${391.16}_{-6.59}^{+6.20}$	−2.95 ± 0.05	4250	${4341.04}_{-178.21}^{+235.08}$	0.080	${0.073}_{-0.019}^{+0.146}$
19	${0.2733}_{-0.0004}^{+0.0003}$	${396.28}_{-5.62}^{+4.78}$	−3.01${}_{-0.05}^{+0.04}$	4500	${4657.20}_{-725.55}^{+29.61}$	0.010	${0.145}_{-0.140}^{+0.592}$
20	0.2733 ± 0.0003	${396.05}_{-5.85}^{+5.23}$	−3.00 ± 0.05	4500	${4664.47}_{-452.27}^{+24.21}$	0.009	${0.283}_{-0.278}^{+0.490}$
21	0.2732 ± 0.0003	${396.10}_{-5.87}^{+5.37}$	−2.99 ± 0.05	4500	${4670.68}_{-20.66}^{+19.06}$	0.005	${0.404}_{-0.385}^{+0.402}$
22	0.2729 ± 0.0004	${400.01}_{-6.01}^{+5.34}$	−3.02 ± 0.05	4500	${4505.79}_{-276.72}^{+102.77}$	0.040	${0.055}_{-0.028}^{+0.052}$
23	0.2730 ± 0.0004	${398.58}_{-5.77}^{+5.35}$	−3.01 ± 0.05	4500	${4456.58}_{-397.22}^{+149.41}$	0.035	${0.033}_{-0.017}^{+0.047}$
24	${0.2734}_{-0.0005}^{+0.0003}$	${395.19}_{-5.78}^{+5.29}$	−3.00 ± 0.05	4500	${4632.19}_{-842.62}^{+49.95}$	0.020	${0.035}_{-0.030}^{+0.603}$
25	0.2724 ± 0.0004	${406.09}_{-5.21}^{+5.52}$	−3.07 ± 0.05	4500	${4570.62}_{-105.43}^{+57.93}$	0.160	${0.400}_{-0.192}^{+0.02}$
26	0.2728 ± 0.0004	${401.39}_{-5.64}^{+5.21}$	−3.03 ± 0.05	4500	${4562.99}_{-105.43}^{+57.93}$	0.139	${0.302}_{-0.149}^{+0.027}$
27	${0.2733}_{-0.0004}^{+0.0003}$	${398.53}_{-5.75}^{+4.93}$	−2.98 ± 0.05	4500	${4501.57}_{-294.34}^{+111.42}$	0.080	${0.065}_{-0.033}^{+0.066}$

Note. The ground truth planetary parameters are given for comparison."Case 0" shows the parameters obtained when a retrieval using ASteRA is conducted on the uncontaminated spectrum as a frame of reference and to verify that ASteRA introduces no intrinsic bias.

Download table as:  ASCIITypeset image

The retrieved parameters obtained for the 27 contaminated spectra when the spot parameters are not fit for are given in Table 4. In comparison to the retrievals presented in Section 3.1, the errors in the retrieved planetary parameters are far more significant, particularly for the highest-activity cases. For the most contaminated case (Case 7), the decision to account for stellar activity, even with the simplified method used by ASteRA, can be the difference in retrieving an approximately solar level water abundance ($\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$ = −3.22 ± 0.05) or an incorrect subsolar water abundance ($\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$ = −5.32 ± 0.06) an underestimation equivalent to >2 orders of magnitude. The planetary temperature (T_P ) is also significantly underestimated with the retrieved temperature of 275 K being 125 K cooler than the ground truth (400 K), which, in the context of a temperate planet, represents a very substantial error (>30%).

Table 4. Retrieved Planetary Parameters Obtained for the Same 27 Cases as in Table 3 but without Accounting for the Presence of Stellar Contamination by Fitting for the Spot Parameters

Case No.	RP (RJup)	T (K)	$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$
GT	0.2730	400	−3.00
0	0.2731 ± 0.0002	${398.40}_{-5.50}^{+5.08}$	−3.00 ± 0.05
1	${0.2746}_{-0.0002}^{+0.0003}$	${390.68}_{-6.30}^{+5.63}$	−3.09 ± 0.05
2	${0.2743}_{-0.0002}^{+0.0003}$	${389.63}_{-5.93}^{+6.06}$	−3.05${}_{-0.04}^{+0.05}$
3	${0.2736}_{-0.0002}^{+0.0003}$	${392.96}_{-5.76}^{+5.61}$	−3.00${}_{-0.05}^{+0.04}$
4	0.2792 ± 0.0003	${365.30}_{-7.10}^{+6.97}$	−3.41 ± 0.05
5	0.2781 ± 0.0003	${368.64}_{-6.93}^{+6.72}$	−3.27 ± 0.05
6	0.2753 ± 0.0003	${381.90}_{-6.23}^{+5.82}$	−3.05 ± 0.05
7	0.2986 ± 0.0001	${274.94}_{-0.28}^{+0.18}$	−5.32 ± 0.06
8	0.2944 ± 0.0001	${274.93}_{-0.41}^{+0.29}$	−4.44 ± 0.05
9	0.2827 ± 0.0003	${330.52}_{-7.58}^{+7.90}$	−3.28 ± 0.07
10	${0.2738}_{-0.0002}^{+0.0003}$	${393.36}_{-5.95}^{+5.17}$	−3.04 ± 0.05
11	0.2737 ± 0.0002	${393.43}_{-5.80}^{+5.35}$	−3.03 ± 0.05
12	${0.2733}_{-0.0002}^{+0.0003}$	${395.31}_{-6.02}^{+5.07}$	−3.00 ± 0.05
13	0.2758 ± 0.0003	${385.69}_{-6.23}^{+6.47}$	−3.24 ± 0.05
14	0.2753 ± 0.0003	${386.20}_{-6.26}^{+6.24}$	−3.16 ± 0.05
15	0.2741 ± 0.0001	${390.73}_{-6.18}^{+5.56}$	−3.04${}_{-0.04}^{+0.05}$
16	${0.2833}_{-0.0003}^{+0.0004}$	${374.04}_{-7.97}^{+7.45}$	−4.16${}_{-0.04}^{+0.05}$
17	0.2819 ± 0.0003	${363.24}_{-7.01}^{+7.63}$	−3.80 ± 0.05
18	0.2774 ± 0.0003	${370.91}_{-6.51}^{+6.50}$	−3.21 ± 0.05
19	0.2735 ± 0.0002	${395.05}_{-5.87}^{+5.17}$	−3.02${}_{-0.04}^{+0.05}$
20	${0.2734}_{-0.0002}^{+0.0003}$	${395.45}_{-5.87}^{+5.07}$	−3.01 ± 0.05
21	0.2732 ± 0.0002	${396.01}_{-5.61}^{+5.06}$	−2.99${}_{-0.05}^{+0.04}$
22	${0.2745}_{-0.0002}^{+0.0003}$	${395.50}_{-6.24}^{+5.09}$	−3.13 ± 0.04
23	0.2742 ± 0.0002	${392.38}_{-5.78}^{+5.42}$	−3.09 ± 0.04
24	${0.2736}_{-0.0002}^{+0.0003}$	${393.47}_{-5.78}^{+5.08}$	−3.01${}_{-0.05}^{+0.04}$
25	${0.2784}_{-0.0002}^{+0.0003}$	${384.85}_{-6.03}^{+5.39}$	−3.61${}_{-0.04}^{+0.05}$
26	0.2775 ± 0.0003	${382.73}_{-6.67}^{+7.00}$	−3.44 ± 0.05
27	0.2752 ± 0.0003	${385.92}_{-6.34}^{+6.21}$	−3.13 ± 0.05

Note. Comparison with the ground truth shows that in cases of severe activity (e.g., Cases 7, 8, and 9) the planetary parameters retrieved are highly inaccurate.

Download table as:  ASCIITypeset image

Conducting two retrievals on the same contaminated spectrum allows us to determine which model is preferred by the data through the comparison of the Bayesian evidence. The Bayes factors (lnB) show a strong preference for the ASteRA retrieval for the majority of cases considered here (Figure 6). The cases in which only a weak preference for the corrected model is seen, or a preference for the retrieval in which contamination is neglected is indicated, correspond to those characterized by small, high-latitude, and low-temperature contrast spots. As the contamination resulting from such spots is minimal, the model neglecting contamination is still able to perform reasonably well under conditions of low activity. In contrast, a penalty is applied to the model accounting for contamination when calculating the Bayesian evidence as a result of its higher dimensionality, explaining why the activity model is not conclusively preferred in the low activity regime despite a spot being present. For these low-activity cases, the planetary parameters are still accurately retrieved even when the presence of the spot is neglected entirely (Figure 7). The Bayes factor shows the strongest preference for the retrieval without the activity correction for the uncontaminated case as one would expect. This reaffirms that no intrinsic bias is introduced by ASteRA and that ASteRA does not find evidence for stellar activity where there is none.

Zoom In Zoom Out Reset image size

In order to adequately understand and interpret the retrieval results we must first consider how each of the three spot parameters will contribute to the contamination factor and therefore the observed spectrum. The size and the temperature contrast of the spot have very intuitive impacts on the contamination spectrum. The contamination factor increases as a function of both as expected. In comparison, the effect of the spot's latitude/position on the stellar disk is more subtle. It also affects only the shortest wavelengths with the largest variations seen at λ ≤ 1μm. The importance of the spot latitude in this study is twofold. First, the projected filling factor (defined as the percentage of the 2D stellar disk that is covered by the spot) decreases as a function of latitude for a spot of a constant radius due to decreases in its 3D geometric projection onto the surface. This reduces the contamination introduced by the spatial coverage of the spot alone. Second, the position of the spot dictates how it will interact with the limb-darkening effect as shown by the normalized intensity profiles (Figure 2) in Section 2.3. We subsequently term this interaction as the limb darkening–spot interplay. The limb darkening–spot interplay is an important consideration as a central spot will remove a greater amount of flux compared to a spot (with an identical filling factor) located at the limb as the intensity maximum occurs at the disk center.

In this section, we focus on how the use of the simplified spot model within ASteRA affects the retrieved planetary parameters (Figure 7) and how these errors are observed to vary as a function of the spot parameters. As already stated in Section 3.1, the largest errors are seen for the most extreme activity cases. Figure 7 allows for the visual comparison of the values retrieved when the spot is fit for (black data points) as opposed to when it is not (red data points). This reiterates the large improvement in accuracy obtained through using even a simple activity correction method over no correction at all. In the worst-case scenario not correcting for stellar activity leads to an underestimation of the water mixing ratio by over two orders of magnitude. Such a large error would significantly impact our understanding of the planet as a whole. Another concerning result highlighted by Figure 7 is that, while the bias introduced by not correcting for stellar activity strongly affects the retrieval accuracy, it does not affect the retrieval precision, as evidenced by the small error bars. As such this high precision could be very misleading for real observations where a priori knowledge of the correct planetary parameters is not known. This is consistent with previous retrieval-oriented works, e.g., Iyer & Line (2020) and reaffirms that caution should be taken when analyzing retrieval results where stellar contamination has not been included, even if the host star is thought to be less active.

Zoom In Zoom Out Reset image size

The parameter that appears to be most influential in the highest-activity cases is the spot latitude, which acts as a proxy for limb darkening (Figure 8). It is intuitive that this spot parameter should have the greatest effect on the residual bias as ASteRA has no way of fitting for the spot's position. A decreasing trend as a function of increasing ϕ_spot is observed in the retrieved planetary temperature, with this being overestimated for the lower-latitude spot cases (Cases 7 and 8) and subsequently underestimated at the highest latitude considered (Case 9). This underestimation can be attributed to the interplay of a high-latitude spot and limb darkening. The reduced flux originating from the quiet photosphere at the limbs acts to reduce the observed spot contrast and thus there is also a reduction in the degree of contamination introduced. In contrast to this, the opposite trend is observed in the retrieved H₂O mixing ratio with this being underestimated at the two lower latitudes and overestimated at the highest. Similar trends in T and $\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$ are observed for the other cases considering large (0.4 R_*) spots (Cases 16, 17, and 18 and Cases 25, 26, and 27 respectively); however, the magnitude of the bias introduced is weaker due to the lower-spot contrasts considered.

Zoom In Zoom Out Reset image size

ASteRA is less successful in recovering the spot parameters (Figure 9). This is likely due to a combination of not accounting for limb darkening and also because many combinations of the three spot parameters are degenerate at low resolution, particularly if the spot in question has a low-to-moderate temperature contrast. An example of how different spot parameter combinations can result in degenerate solutions is given in Figure 10 for clarity. The largest errors and uncertainties in the retrieved T_spot are seen for the smallest spots considered, especially when present at high latitudes. For larger spots, T_spot is generally reasonably well constrained due to the larger contamination effects they introduce. Large errors and uncertainties are also seen in the retrieved F_spot values in several cases. Large error bars point to substantial degeneracy in the small spot cases, whereas in the case of the largest spots, F_spot is often constrained to a higher precision but significantly overestimated. The results of these retrievals indicate that we should be more cautious with retrieved stellar parameters, as the degeneracies between them mean that they are constrained less confidently by the retrieval. Simultaneous, external observations, e.g., at different bandpasses could help break some of these degeneracies.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to attribute the biases seen in the planetary parameters to not accounting for the effect of limb darkening, the worst-case scenarios (Cases 7, 8, and 9) were rerun. For each case, the contaminated spectrum was regenerated with the LDCs set to zero (Figure 11) and a further retrieval was conducted. When noise is taken into account the effect of including versus excluding the limb-darkening effect is really only distinguishable at the shortest wavelengths considered (λ < 1 μm), even for these worst-case scenarios. Ariel in particular will only be able to access this wavelength region through three photometric bands (Tinetti et al. 2021), as such, extracting as much information as possible about the activity of the host star from these three data points will be of utmost importance. The retrieved parameters when the limb-darkening effect is removed are given in Table 5. With the exception of T_spot, which was already well constrained, all other parameters are retrieved more accurately, including F_spot now being constrained correctly. As the input observations are inherently different due to the removal of the limb-darkening effect, unfortunately, we cannot use the Bayes factor as a means of model comparison here as was done in previous sections (Figure 6). However, the fact that both the planetary and stellar parameters are retrieved more accurately when the limb-darkening contribution is removed from the input observation provides compelling evidence in itself that the residual bias seen originates from the limb darkening–spot interplay. For real observations ignoring the effects of limb darkening within activity correction frameworks will unfortunately not be a viable solution as this phenomenon will always be present. As such, we believe that as a community there is a need for us to push toward developing and using stellar activity models in which the limb darkening–spot interplay is accounted for, particularly when dealing with highly active host stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. The Retrieved Spot and Planetary Parameters Obtained for Cases 7, 8, and 9 When the Effect of Limb Darkening is No Longer Present in the Input Spectra i.e., the Blue Spectra in Figure 11

Case No.	Tspot	Input Fspot	Retrieved Fspot	RP (RJup)	T	$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$
GT	3750	⋯	⋯	0.273	400	−3
7	${3700.81}_{-33.72}^{+34.68}$	0.160	0.166 ± 0.002	0.2728 ± 0.0004	${392.47}_{-6.06}^{+5.69}$	−2.91 ± 0.05
8	${3704.40}_{-38.06}^{+33.81}$	0.139	0.144 ± 0.002	0.2728 ± 0.0004	${393.19}_{-5.51}^{+6.37}$	−2.92 ± 0.05
9	${3715.14}_{-58.94}^{+82.23}$	0.080	0.084 ± 0.003	0.2730 ± 0.0004	${395.72}_{-6.15}^{+5.58}$	−2.95 ± 0.05

Note. The ground truth spot and planetary parameters are given for comparison. F_spot varies on a case-by-case basis, as such, the Input F_spot denotes the ground truth value for each case.

Download table as:  ASCIITypeset image

The posterior distributions retrieved for the two input scenarios, spot contaminated with the limb-darkening contribution and spot contaminated with the limb-darkening contribution removed, are given in Appendix A. The greatest deviation between the two posteriors is seen in the case of a central spot. This is intuitive as the spot masks the disk center where a greater proportion of the stellar flux originates from and as such, results in the largest residual contamination due to not encompassing the limb darkening–spot interplay. As the spot is modeled at progressively higher latitudes, less residual contamination remains and the posteriors for the limb darkened, contaminated spectra converge toward the correct values.

In the previous sections, we have focused only on isolated, single spots in order to assess the interactions of the spot parameters and how they influence the observed, contaminated spectrum. In order to better understand the limb darkening–spot interplay (Section 4.1) and the bias it introduces in retrievals (Section 4.4) we conducted three preliminary multiple-spot cases using the same methodology to explore how this interplay may differ when more than one spot is present. For comparability, the spots have the same temperature (T_spot = 3750 K) and total filling factor (F_spot = 0.16) as the single-spot Case 7. The multiple-spot cases we consider are a two-spot case and two cases characterized by 10 smaller spots. These 10 spots are arranged in a random configuration on the stellar disk in the first scenario, and occupy a preferred active latitude centered around ϕ = 15°N in the second scenario. The only strict requirement for these cases is that all of the spots remain unocculted. The resultant contaminated spectra are shown in Figure 12 alongside the uncontaminated spectrum and the Case 7 contamination spectrum, with and without the contribution of the limb-darkening effect, as in Figure 11. These spectra show that, as the number of spots increases and they are distributed over a larger range of μ values (i.e., located at different distances from the disk center) the limb darkening–spot interplay and its contribution to the contamination spectrum are minimized. A table of the retrieved parameters is given in Appendix B, alongside the posterior distributions for the three multiple-spot cases considered.

Table 6. The Retrieved Planetary and Spot Parameters for the Multiple-spot Cases Considered in Section 4.5 Compared to the Ground Truth Values and Case 7, a Single Spot of the Same Temperature and Filling Factor

	R (RJup)	T (K)	$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$	Tspot (K)	Fspot
Ground truth	0.273	400	−3	3750	0.16
Case 7	${0.2698}_{-0.0005}^{+0.0004}$	${429.02}_{-5.45}^{+5.98}$	−3.22 ±0.05	${3787.50}_{-59.60}^{+41.23}$	${0.234}_{-0.006}^{+0.002}$
Two spot	0.2705 ± 0.0004	416.80 ± 5.10	−3.11${}_{-0.05}^{+0.04}$	${3702.22}_{-36.41}^{+37.76}$	0.203 ± 0.002
Multiple spot—random	0.2721 ± 0.0004	${400.92}_{-6.68}^{+5.81}$	−2.97 ±0.05	${3701.30}_{-34.57}^{+34.30}$	0.178 ± 0.002
Multiple spot—active latitude	0.2716 ± 0.0004	${405.52}_{-5.58}^{+5.07}$	−3.01 ±0.05	${3700.95}_{-34.57}^{+36.00}$	0.185 ± 0.002

Download table as:  ASCIITypeset image

These additional multiple-spot experiments highlight that the large, high-contrast single spots considered as the main focus of this study likely do represent the worst-case scenario when considering the additional complication arising from the limb-darkening effect. Observationally, the regimes in which the limb darkening–spot interplay will really start to matter will therefore be those in which the system geometries most closely replicate those in Cases 7 and 8. Large high-latitude and polar spots have been proposed to exist on several stars (e.g., Järvinen et al. 2018; Almenara et al. 2022; Strassmeier et al. 2023) although any biases that these could introduce will be lessened if the star's rotation axis is not significantly inclined with respect to the line of sight as the spot will appear at the limb. The worst-case scenario would occur for a system where the host star both possesses a polar spot and is inclined relative to the transiting planet and observer. In such a scenario the polar spot could manifest as a large central spot. Examples such as the HAT-P-11 system, which consists of an active K4-dwarf hosting two highly misaligned planets (e.g., Sanchis-Ojeda & Winn 2011; Morris et al. 2017; Yee et al. 2018) and the almost pole-on, solar-type star τ Ceti, which is frequently included in target lists for exoplanet searches (Korolik et al. 2023) both indicate that observing a system with such geometry in the future cannot be ruled out.

Our recommendation for dealing with this worst-case, high-activity regime, where it will be necessary to include the limb darkening–spot interplay in order to fully negate the bias introduced by stellar contamination, is to improve our retrieval stellar models so that they are capable of also parameterizing the position of spots on the stellar disk. This is the focus of ongoing work, however, as with increasing the dimensionality of any model, we cannot ignore the risk of injecting intrinsic bias if the additional complexity is not physically motivated. For real observations, there will also be the added challenge of not having any a priori knowledge of spot positions as the stellar disks are not spatially resolved. This work has shown that at first order using the simpler, two-parameter spot prescription as is done with ASteRA is sufficient to remove the majority of the bias and enable a good understanding of the planet's atmospheric properties. A push toward a better understanding of the limb-darkening effect, whether in the presence of stellar activity or not, will also be highly beneficial. This will be especially important for later-type stars where offsets due to the choice of limb-darkening treatment appear to be inevitable (Patel & Espinoza 2022).

All of the retrievals presented in the previous sections have been conducted assuming a single T_spot. In order to explore the validity of this assumption, we conducted preliminary studies into modeling the contamination introduced by spots with radial temperature variations. This section examines the resulting impact on retrieval performance. We assume the same spot geometry as was used in the highest contamination, worst-case scenario (Case 7). The spot is now separated into a cooler central umbra characterized by a temperature T_umbra = 3750 K, and a surrounding penumbra characterized by a temperature T_pen = 4250 K. This separation into two regions of distinct temperatures is consistent with observed sunspots and with the 3D MHD models for later-type stars produced by Panja et al. (2020). We explore two different cases in which the total spot-filling factor (F_spot) is kept constant but the relative fractions of the umbra (F_umbra) and penumbra (F_pen) are varied. Table 7 shows the input parameters for these two cases. In doing this, we explore the retrieval response to contamination effects resulting from a spot that is characterized by two different temperatures/SEDs when the retrieval model is only capable of fitting this contamination using a single T_spot value.

As expected, this further discrepancy between the forward model and the retrieval model introduces an additional source of error to the retrieved planetary parameters. The retrieved parameters are given in Table 8 alongside the equivalent single T_spot case (Case 7) for comparison. The posterior distributions for the two-spot temperature variation retrievals are given in Appendix C. The retrieved spot parameters provide insight into how the retrieval has attempted to correct for the dual temperature spot. In both cases, the retrieval is best able to fit for the contamination with a moderate T_spot close to T_pen, and an overestimated F_spot. The overestimation of the spot-filling factor accommodates the increased contamination due to the higher-contrast umbra. The posteriors show an increasing correlation and degeneracy between T_spot and F_spot compared to the single-T_spot cases. The introduction of radial temperature variation slightly enhances the biases in the retrieved planetary parameters that were observed for the highest-activity cases. Intuitively, the greater effect is seen when considering a spot with a higher umbra–penumbra ratio, resulting in an overestimation of the retrieved planetary temperature of ∼47 K and an underestimation of the water mixing ratio of 0.39 magnitude. Although this error is non-negligible, importantly it is not limiting. Retrievals conducted in the presence of such contamination would still be useful and informative.

Table 8. Retrieved Planetary and Spot Parameters for the Two Radial Temperature Variation Spot Cases Considered in Section 4.6 Compared to the Highest-activity, Single-temperature Spot Case (Case 7) and the Ground Truth Values

	R (RJup)	T (K)	$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})$	Tspot (K)	Fspot
Ground truth	0.273	400	−3	⋯	⋯
Case 7	${0.2698}_{-0.0005}^{+0.0004}$	${429.02}_{-5.45}^{+5.98}$	−3.22 ±0.05	${3787.50}_{-59.60}^{+41.23}$	${0.234}_{-0.006}^{+0.002}$
Temperature variation Case 1	${0.2721}_{-0.0003}^{+0.0004}$	${446.52}_{-5.66}^{+5.57}$	−3.39 ±0.04	${4201.09}_{-34.46}^{+33.84}$	0.254 ± 0.003
Temperature variation Case 2	${0.2715}_{-0.0004}^{+0.0003}$	${440.65}_{-5.84}^{+5.86}$	−3.33${}_{-0.04}^{+0.05}$	${4299.21}_{-33.68}^{+33.98}$	0.278 ± 0.004

Download table as:  ASCIITypeset image