© ESO 2020

1 Introduction

After the ground-breaking work done during the last decades in the search for exoplanets using both the Doppler technique and photometric transits (Mayor & Queloz 1995; Charbonneau et al. 2000; Anglada-Escudé et al. 2016; Gillon et al. 2017), the community has expanded its focus to include the characterisation of known transiting planets through the study of their atmospheres. The central idea behind transmission spectroscopy is that the planetary transit depth has a chromatic dependence since its atmosphere selectively absorbs certain wavelengths. The measured transit depth as a function of wavelength constitutes the planet’s transmission spectrum and opens up the possibility of studying the chemical (composition, abundances) and physical (thermal structure, pressure profile) properties of the atmosphere. This has been revealed as a very successful technique (Seager & Sasselov 2000; Tinetti et al. 2007; Fortney et al. 2010; Burrows 2014; Sing et al. 2015; Tsiaras et al. 2019). The optimal candidates for transmission spectroscopy are low-density planets with atmospheres of low mean molecular weight (hydrogen-dominated) and high equilibrium temperature, thus favouring hot, giant gaseous planets (Stevenson 2016; Kreidberg 2018). Even for these candidates, the observed transit depth effects induced by the planet’s spectral features are the order of 10⁻³, which makes their detection very challenging, particularly in the case of active host stars, where other sources of variations are expected.

Future space-based telescopes such as the James Webb Space Telescope (JWST; Gardner et al. 2006) and the Ariel (Tinetti et al. 2018) missions are designed for this purpose. The latter will carry out a comprehensive, large-scale survey of the chemical compositions and thermal structures of the atmospheres of ~1000 known transiting exoplanets in the optical and near-infrared (NIR) wavelengths (0.5–8 μm) following itsexpected launch in 2028 (Tinetti et al. 2018; Encrenaz et al. 2018).

Magnetically-driven active regions on the surface of a star, such as dark spots and bright faculae, can produce significant alterations on the measured planetary transit depths (Lagrange et al. 2010; Meunier et al. 2010a; Barros et al. 2013; Oshagh et al. 2014). Such inhomogeneities on the stellar photosphere can induce chromatic effects (Sing et al. 2015), which, in the case of unocculted spots, can be very similar to the signature of atmospheric Rayleigh scattering (Rabus et al. 2009; McCullough et al. 2014). Rackham et al. (2018) studied their influence in M dwarfs and find it to be up to ten times larger than the effect expected from atmospheric compounds even at NIR wavelengths. Although Rackham et al. (2019) find the influence of surface phenomena in FGK stars to be generally lower and possibly only measureable for active stars, those studies challenge the reliability of retrieved transmission spectra, at least under certain circumstances. Clouds and hazes can also produce chromatic effects, introducing a Rayleigh-like slope and grey opacity, respectively (Pinhas & Madhusudhan 2017). Distinguishing the different contributions to wavelength-dependent transit depth variations is therefore of crucial importance, and this calls for sophisticated modelling of photospheric inhomogeneities.

Several approaches to model the effects of stellar activity on time-series spectrophotometry have been presented. Examples are the SOAP (Boisse et al. 2012) and SOAP 2.0 (Dumusque et al. 2014) codes, which are conceived to provide estimates of the effects of spots and plages on photometry and radial velocities (RVs) by surface integration of a pixel grid using solar observed cross-correlation functions (CCFs). Similar to this is the maximum entropy regularisation method used by Lanza et al. (2011) to model a unique and stable solution for the flux variations due to discrete elements of active regions. An alternative empirical approach is the FF’ method (Aigrain et al. 2012), which uses the flux time series to predict the RV jitter. Herrero et al. (2016) presented the StarSim code as a more sophisticated methodology to simulate light curves in various passbands and spectroscopic time series of active rotating stars taking into account the effects of stellar spots and plages. StarSim is based on surface integration techniques, taking into account the particular properties of each element of a fine stellar surface grid such as the effective temperature, the limb darkening and convection effects.

We describe the new version of StarSim in Sect. 2, detailing the designed method to retrieve the properties of a star and its activity by solving the inverse problem with real time-series data. The aim of this study is to obtain an activity model which comprises 1) a stellar surface map and, 2) a set of stellar parameters. This methodology is implemented in a test case using observational multiband photometry of the planet host WASP-52, as shown in Sect. 3. In Sect. 4, we use the retrieved activity model for WASP-52 to study the different sources of chromatic effects on the transit depths of its known planet WASP-52 b. In Sect. 5 we conclude and discuss the possibilities of this methodology in correcting transit spectroscopy data affected by stellar activity.

2 StarSim model

2.1 The forward model

The StarSim modelling code generates the effects of evolving dark spots and bright faculae of a rotating star on photometric, RV and activity indices time-series data of magnetically-active stars. We explain the most important points of the modelling approach below, but we refer the reader to the more detailed description in Herrero et al. (2016). The code divides the stellar photosphere into a grid of elements, each of which is assigned an effective temperature value depending on its nature and according to the distribution of stellar activity features. We use T_ph for an immaculate element of the photosphere, T_sp for those occupied by cool spots, and T_fc for those corresponding to bright faculae, and assume generally that T_sp < T_ph < T_fc. We attribute to each surface element the corresponding theoretical BT-Settl spectrum (Allard et al. 2013) generated with the Phoenix code. In this approach the only difference between surface features is its effective temperature (and not other spectral effects related, e.g. to the stellar chromosphere), as this is the parameter that dominates their chromatic signature on the simulated data. Activity features on the photosphere are described by a set of magnetically-active elements. Each active element is represented with a circular spot that can be surrounded by a concentric facular area. The model can take into account arbitrary shapes of the active regions by grouping several smaller active elements together. Throughout the text, we indistinctly use the word spot or the expression active element.

As a simplification aimed at a faster execution of the code, the properties of the stellar surface, such as the limb darkening, the radial velocity, and the projection angle towards the observer, are assumed to be constant inside each active element. Therefore, active elements are always assumed to be small as they are associated with the properties of a single pixel on the surface grid.In this version of StarSim we employ active elements having circular shape. Besides the different temperature, a set of six parameters are used to characterise an active region: the time of appearance of the active element, its lifetime, its location on the surface of the star (latitude and longitude), its angular radius, and the facula-to-spot area ratio Q. The last one is the only parameter assumed to be the same for all active elements. The model includes all the relevant geometric, physical and wavelength-dependent features related to the presence of active regions, such as limb darkening effects computed from Kurucz ATLAS9 models (Kurucz 2017), limb brightening in the case of faculae (Meunier et al. 2010a), convection effects inside and outside the active regions (see Herrero et al. 2016, for a detailed description of the convection model), and the projection of each surface element towards the observer.

Photometric and spectroscopic time-series are obtained by integrating over all the grid elements at each time step, taking into account the rotation of the star and the time-evolution of the active regions. When simulating photometrictime-series, the integrated spectrum of the projected stellar surface is computed at each time step, and then multiplied by the defined filter passbands in order to retrieve multiband fluxes and produce light curves. By construction, the resulting light curves are relative to the immaculate photosphere (i.e. immaculate flux = 1). To compare with observational measurements, a zero-point correction needs to be applied to each band. This is done by normalizing both the observed light curve and the model light curve to the weighted average of all measurements. In the case of the simulated light curve, only the epochs with observational measurements are used. For the calculation of RVs and other spectral indices related to CCFs, StarSim initially generates the CCFs produced from the spectrum of a single photosphere, spot and facular element, and then integrates the entire visible surface using CCFs instead of spectra. This is necessary to speed up the execution of the code and make it possible to handle with normal desktop computers. A slow rotator template or a user-defined mask of spectral lines can be used to calculate the CCFs. In a first step of each time-series simulation, the contribution of the entire immaculate photosphere is computed without any active elements. Then, for each time step of the simulated data series, the contribution of the active regions is added by considering their individual location on the observed stellar disc.

2.2 Analytical foundations of the inverse problem

Before moving on to discuss how StarSim handles the inversion of photometric light curves, we shall present the analytical foundations of the model and the relevant variables. We consider a stellar surface with active regions covering a total projected filling factor δ_sp, which is made of two components: one that produces flux variations over time due to rotational modulation, δ_M, and another one that stays constant as the star rotates (e.g. a polar cap or a uniformly spotted latitudinal band), δ₀. To demonstrate analytically that light curves can carry information of both spottedness levels we assume a simple model consisting of a star with a circular spot on its equator, with a surface δ_M, and a polar cap covering a surface δ₀, as illustrated in Fig. 1.

The amplitude of the photometric modulation can be estimated as the ratio of the flux of the star when the modulating spot is out of view, (1)

with respect to the phase when the spot is at the centre of the visible stellar disc, (2)

where F_ph and F_sp are the surface fluxes of the spotted and immaculate surface of the star, respectively, which depend on the temperature of each surface element (T_ph and T_sp) and the spectral band (). The relative amplitude of the photometric variability can be computed as (f₁ − f₂)∕f₁, which is equivalent to the measurement provided by observational light curves. Rearranging Eqs. (1) and (2), the photometric amplitude can be written as, (3)

where (4)

with ΔT_sp being the difference between the photospheric and spot temperature (ΔT_sp = T_ph − T_sp). We assume that the effective temperature of the star can be estimated independently. These equations illustrate that the photometric amplitude increases linearly with the modulating filling factor δ_M, and decreases with a combination of the non-modulating filling factor δ₀ and the brightness contrast between the photosphere and the spots Φ, which is, in turn, a function of the photometric band () and the temperature contrast (ΔT_sp). We remark here that these quantities satisfy some constraints, namely, 0 ≤ δ₀ ≤ 1, 0 ≤ δ_M ≤ 1, 0 ≤ δ₀ + δ_M ≤ 1 and Φ > 1, for cool spots.

Figure 2 illustrates Eq. (3) and its dependence on the three independent parameters. If only information from a single band is available (e.g. typical data from exoplanet surveys), the linear dependence of the light curve amplitude with δ_M permits the determination of the modulating filling factor provided the non-modulating filling factor is neglected and the spot temperature contrast is adopted as an external constraint. This has been a common practice in the literature. In case two or more bands are available, the possibility of determining another variable such as ΔT_sp or δ₀ arises. For typical spot temperature contrasts and visible bands, it can be shown that Φ is significantly greater than δ₀ and therefore the former parameter dominates. Thus, from two or more photometric bands (preferably covering a large interval in wavelength, i.e. large Φ variation), the simultaneous determination of δ_M and ΔT_sp becomes possible. A practical application can be seen in Mallonn et al. (2018). In the particular case of two bands, the solution becomes bi-valuate, with two possible ΔT_sp reproducing the amplitude difference, as shown in Fig. 2b. However, if a third band is available, such degeneracy can be broken, as illustrated by the gray lines in the figure. Equation (3) further shows that, from three bands or more, one can theoretically determine at the same time the three relevant variables. For that to be possible, the photometric information needs to be of sufficient precision to discriminate the changes induced by each variable. The curvature in Fig. 2c (i.e. the variation in the amplitude difference) is what makes it possible to determine the non-modulating filling factor from multi-colour photometry. Nevertheless, the scale of the variations makes reliable estimates of δ₀ very challenging for typical ground-based photometric precisions.

We remark that we have presented a simplified version of the problem, defined by only two epochs (maximum and minimum light). However, photometric time series also carry information on the relevant parameters because of the correlations present among the different measurement epochs, thus adding additional constraints to the spot properties. The formulation discussed here proves that it is possible to simultaneously determine the total spot filling factor and spot temperature contrast as long as good multiband photometric data are available, therefore providing theoretical foundation to the inverse problem.

Fig. 1 Example of a stellar photosphere including modulating and non-modulating spots. The model is composed of a polar cap, which does not contribute to light curve modulation, with a projected filling factor δ0, and an equatorial small spot with a maximum projected filling factor δM. F(⋅) is the brightness of each surface element.

2.3 The inverse problem

The most recent StarSim version can perform the inverse problem. The goal is to obtain the underlying properties, that is to say a stellar activity model as described by the parameters of the star and its surface map, that reproduce the observed time-series data X. This type of non-linear problem can be expressed as (5)

where X is the time-series data, F is the activity model, and θ is a set of stellar parameters. The surface map is the set of parameters that describe the surface distribution, sizes, and lifetimes of all the active elements considered, each of them defined as a small circular spot surrounded by a bright facular region. Finally, ɛ is an additional noise term, or jitter, that we adopt as uncorrelated and following a Gaussian distribution (white noise).

The inverse problem consists in finding the surface map that best reproduces the data X for a given θ, expressed as (6)

where means the optimal surface map linked to a specific set of stellar parameters. Subsequently, this map and its associated parameters provide an optimal fit to the photometric and spectroscopic time-series data when applying the forward model (Eq. (5)).

Fig. 2 Graphical representation of Eq. (3) for the parameters of interest assuming a stellar photospheric temperature of Tph = 5000 K and a black-body for the spectral energy distribution of each surface element. (a) Photometric amplitude as a function of the modulating filling factor for different spot temperature contrasts. A non-modulating filling factor δ0 = 0.05 is employed, but adopting other values has negligible effect at this scale. The gray area covers the spectral range from 400 nm to 1 μm. (b) Difference in activity-induced amplitude for photometric bands BV R with respect to the I band as a function of spot temperature contrast. A value of 0.05 is adopted as modulating filling factor (δM) and examples for non-modulating filling factors of δ0 = 0 and 0.1 are shown. The gray lines illustrate how the bi-valuate nature of the temperature contrast effect can be resolved if more than two bands are used. (c) Difference in activity-induced amplitude for photometric bands BV R with respect to the I band as a function of the non-modulating filling factor. A modulating filling factor δM = 0.05 and a spot temperature contrast ΔTsp = 600 K are adopted.

2.3.1 Objective function

The statistical function to optimise, or figure of merit, is a linear combination of the logarithmic likelihood function of all the time-series data defined as (7)

where is the log-likelihood of the fitted model according to the observational data X_j for the jth observable. The quantities a_j are the weights associated to each set of observables (here we assume a_j = 1). As shown inEq. (5), we consider the simplest case of non-correlated Gaussian uncertainties (white noise). The likelihood function is then written as (8)

where is the StarSim-generated model of the jth time-series observable with N measurements, and y_i are the observational data points. σ_i is the nominal error of the measurement i, and s_j is a quadrature-added jitter, that accounts for a possible incompleteness of the model, underestimated uncertainties or traces of correlated noise.

2.3.2 Optimizing the surface distribution of active regions

Surface maps describing the distribution, size and evolution of the active elements are obtained through maximisation of the figure of merit (Eq. (7)) given a fixed set of stellar parameters, θ, typically including the rotation period, P_rot, the spot temperature contrast, ΔT_sp, and the facula-to-spot area ratio, Q. Thus, the optimal surface map is given by (9)

where is the optimal surface map constrained to a specific set of parameters and is any surface map among all possible maps, . In order to optimise the expression in Eq. (9), we implemented a Monte Carlo Simulated Annealing optimisation algorithm (hereafter MCSA, Kirkpatrick et al. 1983).

The surface map optimisation starts with a random distribution of a fixed number of active elements, each one characterised by 5 adjustable parameters: time of appearance, lifetime, latitude, longitude, and angular radius. Active elements do not appear or disappear abruptly. They are assumed to grow or shrink linearly in radius at a rate of 0.5 deg day⁻¹ (see Herrero et al. 2016, for additional details). The value of Q is assumed to be the same for all spot elements. The total number of active regions is a parameter of StarSim and can be defined by the user. However, its is advisable to limit the maximum number of active elements thus forcing the model to retrieve simpler surface maps, avoiding the appearance of a large number of very small spots, especially when modelling noisy time-series data.

The optimisation process follows successive iterations where the algorithm randomly selects one of the adjustable parameters from a randomly-selected active element and modifies it slightly. Then, with this new map, the forward problem is re-computed and is evaluatedfor each time-series dataset. If such perturbed map improves the quality of the fit as given by Eq. (7), the change is accepted. Otherwise, it is accepted with a probability , where β is a parameter that grows in each step. With this strategy, the optimiser avoids getting trapped in local maxima. A more detailed explanation of the implementation of MCSA in StarSim can be found in Herrero et al. (2016).

Due to the large number of parameters describing the stellar surface map and the intrinsic randomness of Monte Carlo algorithms, several maps retrieved with the same θ may not be identical in spite of producing very similar time-series datasets with the forward model. This effect can be mitigated by performing a number of solutions with different initial spot maps and subsequently exploring the variance of the likelihood statistic. For each target and set of observations, these tests can help the user to define parameters such as the number of iterations of the MCSA algorithm per β step and the number of active regions on the stellar surface. The execution time (which increases linearly with both parameters) and some regularisation criteria (selecting the minimum number of spots to avoid unnecessarily complex surface structures) need to be factored in when deciding on the optimal procedure. Once these parameters are adopted and a large number of inversions with varying initial spot conditions are run, the final product of the inversion for a fixed set of parameters, θ, is a spot map calculated by averaging the total sample of optimal maps. The result is a smooth time-evolving surface map that contains valuable information about the surface distribution of the active regions, their typical lifetimes and possible differential rotation.

2.3.3 Optimizing stellar parameters

Equation (9) describes the optimisation of the surface maps when fixing a set of stellar parameters θ. However, tests show that has a dependence on θ. Due to the non-linear and multimodal nature of the problem, and θ are strongly coupled, and small variations on the parameters potentially imply large changes in . Therefore, fitting both sets of parameters simultaneously is computationally challenging. Therefore, we choose a two-step approach. The retrieval of optimal parameters in this context belong to the class of noisy optimisation problems (e.g. Grill et al. 2015), which are characterised by long evaluation times and noisy outputs for the objective function. In our particular case, we consider a large number of randomly-generated parameter sets drawn according to a prior distribution (generally uniform) and calculate the inverse problem to produce an optimal surface map starting from random initial conditions for each of the parameter sets. Using their values, a certain interval enclosing the best statistically-equivalent solutions can be defined. From the selected best solutions, the optimal parameter set and the corresponding uncertainties can be determined. The procedure can be expressed as (10)

where is the threshold defining equivalent solutions. Then, the median and 68% confidence interval of are used as the best estimates of stellar parameters and their uncertainties.

2.3.4 Degeneracies of the inverse problem

The methodology of StarSim to calculate the inverse problem implies accounting for a large number of parameters and potential correlations among them, which may produce degenerate solutions. This is especially important when only a particular type of data is considered (i.e. photometry or spectroscopy). A number of different active region configurations can yield very similar simulated datasets, with figures of merit that are not significantly different. This can be easily illustrated by the example of considering a system with a single spot on an equator-on star (i_star = 90°). In such case, the location in either hemisphere produces identical results for all datasets (hemispherical degeneracy). Also, spots with different sizes at different latitudes can produce similar effects, as they may yield the same projected area towards the observer. This results in size-latitude correlations. Such type of weak degeneracies can be solved with an accurate modelling of limb darkening. However, the uncertainty of photometric and spectroscopic measurements usually make limb darkening variations indistinguishable. This is why, especially for stars that are nearly equator-on, we cannot retrieve the latitude distribution of the active regions but only the stellar longitudes occupied by spots. Therefore, surface maps can only be shown as the filling factor of active regions projected on the longitude axis (see Sect. 3.3.3 for the case study presented here). Also, when polar and circumpolar active regions are present, their integrated effect can result in a constant offset of the photometric observables, which we have shown that can be difficult to identify. The inverse problem in case of stars that are not equator-on is particularly challenging when only photometry is available.

Another important potential degeneracy is caused by the correlation between the size of the active regions and their temperature contrast. However, there is potential relief to such degeneracy because the temperature contrast of spots and faculae introduces a chromatic effect, as we have shown above. This can be measurable when multiband time-series photometry is available. When cool spots are present on the stellar surface, the amplitude of the photometric variability is larger at bluer wavelengths, and lower in the red and infrared bands. Such chromatic signature can be studied by performing the inverse problem with StarSim on multiband photometric measurements, thus allowing to derive the temperature contrast of the active regions with respect to the photosphere and thus break the temperature-size correlation.

Finally, different combinations of the filling factor of spots and the amount of faculae (Q) can also result in similar solutions for the modelled observables. This is because bright regions can partially compensate the effects of dark spots. However, the presence of faculae can also be identified when analyzing multiband photometric time series, with proper modelling oflimb brightening (as opposed to limb darkening for dark spots) and the corresponding chromatic signature, which differs from the one produced by cool spots. Also, due to the blocking of convection by active regions and the subsequent decrease of the net convective blueshift in RVs, the spot-facula degeneracy can also be broken when contemporaneous photometry and high-precision spectroscopy are available.

2.3.5 Toy model inversion

As a further test of the validity of our approach, we conducted an inversion exercise using data generated with the forward functionality of StarSim. We considered astar with T_ph = 5000 K, i = 90 deg, P_rot = 15 days and with its surface covered by 5spots having a ΔT_sp = 700 K and no faculae (Q = 0). We considered spots of various sizes, placing them at different latitudes and longitudes in such a way that a photometric modulation is present and also ensuring that spotted regions are visible at all times (i.e. projected spot filling factor neverbeing zero), to produce a nonzero non-modulating filling factor. We assumed the spots to be of constant size and static in the frame of reference of the rotating star. We generated synthetic light curves in the UBV RIJ bands covering a timespan of 90 days (6 rotations) and with a uniform observational cadence of one measurement every 0.5 days. Noise was added to the photometric measurements following a heteroskedastic scheme, consisting of non-correlated Gaussian noise with (i.e. ≈1000 ppm).

The resulting 6 light curves were inverted following the scheme sketched in Fig. 3. The fitted parameters were P_rot, Q, and ΔT_sp. We ran several thousand inversions exploring the relevant parameter space, all starting from a random map. We considered three different assumptionsregarding the number of active regions in the model, namely 5, 10, and 20, to evaluate the impact on the quality of the solution andon the resulting filling factor. We furthermore considered that spots have a finite lifetime and that their growth or decrease rate is 0.5 deg × day⁻¹. Each spot was characterised by six parameters (appearance time, lifetime, longitude, latitude, radius) that were allowed to vary during the optimisation process. We established a likelihood criterion to select the best maps and solutions based on the inherent scatter of the solution when performing a large number of inversions from random maps and assuming the correct parameters. This yielded a number of 100−150 good solutions for all three cases. The statistical results of such solutions are given in Table 1, and a graphical illustration for the 10-spot case is shown in Fig. 4.

The results of the inversion tests allow assessing the retrieval performance of the StarSim model. As can be seen in Table 1, all fitted parameters are retrieved within the uncertainties regardless of the number of spots assumed. No significant bias is observed except for a slight tendency to underestimate the spot temperature contrast. However, additional tests showed that this may be caused by the relatively low number of solutions used and therefore it should not be reason for concern. It is interesting to note that the spot filling factor is also accurately obtained, including the non-modulating fraction. This indicates that the algorithm is able to reduce the size of the spots and place them appropriately to reproduce the correct spottedness of the star. We do not see a dependence on the number of active regions except for a slight tendency to overestimate the filling factor when 20 spots are used. The graphics in Fig. 4 show the high quality of the multiband fits and the low dispersion of the solutions. Also, the comparison of the input latitude-projected filling factor with the inversion results on the right panels of the figure indicates a very precise retrieval of both the longitudes and sizes of the active regions. We additionally ran tests using only the BV RI passbands, thus restricting the wavelength baseline. The solutions converged well only with slightly larger uncertainties, as expected. In all cases, the input parameters were well within the error bars.

Fig. 3 Flow chart of the parameter retrieval approach as presented in Sect. 2.3.

3 An activity model for WASP-52

As an example of the use of StarSim for the inverse problem, we use multiband photometric time-series data of the active planet host star WASP-52 to retrieve an optimal set of stellar parameters and reconstruct a time-variable probability map of the distribution of active regions.

3.1 The WASP-52 exoplanetary system

WASP-52 isa young and active K2 V star hosting an inflated Jupiter-sized exoplanet with an orbital period of 1.75 days (Hébrard et al. 2013). The inclination of the planetary orbit i is such that results in a transit of the planet. Hébrard et al. (2013), Louden et al. (2017), Mancini et al. (2017), May et al. (2018), Oshagh et al. (2018), and Öztürk & Erdem (2019) described and refined the different planetary parameters using light curves and the Rossiter-McLaughlin effect (see Table 2), whereas Kirk et al. (2016), Louden et al. (2017), Chen et al. (2017), and Alam et al. (2018) described its cloudy sodium-bearing atmosphere. In many of those studies, in-transit anomalies appearing on the transit light curves and the effects of active regions on the stellar photosphereare discussed. May et al. (2018) and Bruno et al. (2018) quantify the differences of the chromatic effect on transit depths from a spotted and unspotted photosphere. The different values found for the stellar rotation period P_rot and the spot temperaturedifferences ΔT_sp, which are fitted parameters of our inverse problem, are given in Table 3.

In Alam et al. (2018), a transmission spectroscopy study was done using three transits observed in the visible and NIR wavelengths with HST/STIS, combined with Spitzer/IRAC photometry. A stellar activity correction was applied by fitting the baseline flux from ground-based photometry using a quasi-periodic Gaussian process. The effective temperature of the star spots was assumed to be 4750 K, that is, 250 K cooler than the photosphere. A recent atmosphere study of WASP-52b by Bruno et al. (2020) combine the STIS and IRAC data from Alam et al. (2018) and HST/WFC3 from Bruno et al. (2018). The spot effects were corrected using the approach described in Rackham et al. (2017) for an heterogeneous photosphere. In both cases, the stellar photosphere was assumed to be dominated by cool spots. It is worth emphasising that Bruno et al. (2020) present a joint fit of the atmospheric model and a stellar contamination correction, parameterised by the fraction of stellar surface occupied by activity features and their temperature. Such approach may lead to a biased solution since genuine planetary features could be modelled as stellar contamination.

Fig. 4 Left panel: StarSim model fits to multiband synthetic light curves assuming the parameters in Table 1, for a subset of 50 days. Solid curves represent the mean of ~20 optimal solutions of the inverse problem. The gray line at the bottom is the projected spot filling factor of the maps, also showing the mean and 1σ shaded band. The error intervals of the photometry bands are very small and difficult to recognise. Right panels: longitudinal spot filling factor projected on the stellar equator as a function of time, for both the input (bottom) and the retrieved (top) spot maps. Only stellar longitudes visible at each time are shown (hence the band structure). The colour scale indicates the fractional latitude-projected spot coverage for each 1-degree longitude bin.

3.2 Photometric observations

Photometric light curves obtained from two different ground-based observatories, the 1.2-m STELLA telescope at Izaña Observatoryin Tenerife and the 0.8-m Joan Oró telescope (TJO) at the Montsec Astronomical Observatory in Catalonia, were used in this work. WASP-52 was observed as part of our monitoring survey VAriability MOnitoring of exoplanet host Stars (VAMOS, Mallonn et al. 2018). Measurements cover a time interval of more than 500 days in two different observing seasons (2016 and 2017), and were obtained contemporaneously with both telescopes using four different passbands (BV RI).

The STELLA telescope and its wide-field imager WiFSIP (Strassmeier et al. 2004, 2010) obtained nightly blocks of three individual exposures in Johnson B and three in Johnson V from May 2016 until January 2018, completed by three exposure in Cousins I since September 2017. The detector is a single 4k × 4k back-illuminated thinned CCD with 15 μm pixels, providing a field of view of 22 × 22 arcmin. Thetelescope was slightly defocused to minimise scintillation noise and improve the quality of the photometry (Mallonn et al. 2016). The data reduction employed the same ESO-MIDAS routines used for previous monitoring programmes of exoplanet host stars with STELLA/WiFSIP (Mallonn et al. 2015, 2018). Bias and flat-field correction was done by the STELLA pipeline. We used SExtractor (Bertin & Arnouts 1996) for aperture photometry and extracted the flux in elliptical apertures (SExtractor aperture option MAG_AUTO). Data points of low quality due to non-photometric conditions were discarded. Weighted nightly averages were considered for each filter, resulting finally in a total of 112 measurements in the B filter, 110 in V, and 16 in I. The use of nightly averages is justified because the noise floor for each night is typically limited by systematic errors in the measurement procedure and calibration of ground-based photometric observations. We avoid giving excessive weight to nights with larger number of measurements by computing the nightly averages and adding a jitter term to account for random errors not included in the model of the observations such as night-to-night calibration errors. This strategy is appropriate because we are interested in rotational modulation timescales, with the few-hour time domain being irrelevant.

The TJO telescope provided photometric data using its main imager MEIA2. The instrument has a single 2k × 2k back-illuminated thinned CCD yielding a field of view of 12.3 × 12.3 arcmin and a resolution of 0.36 arcsec pixel⁻¹. The images were calibrated with darks, bias and flat fields using the ICAT pipeline (Colomé & Ribas 2006) and differential photometry was extracted using AstroImageJ (Collins et al. 2017). A total of 66 and 77 epochs, that is, weighted nightly averages, were obtained with the Johnson-Cousins B and R filters, respectively.

3.3 Photometric inverse problem

3.3.1 Fixed StarSim parameters

The basic stellar parameters (relevant to select the appropriate atmosphere models) were adopted to be T_eff= 5000 K and logg = 4.5, which are close to the literature values (see Table 2). The stellar inclination was chosen to be 90°, which is a reasonable assumption given the measured values of the planet’s orbital inclination and spin-orbit angle. The lifetime of the active elements was assumed to be Gaussian-distributed around 250 ± 100 days. The choice of this value is done considering the full timespan of the series and the fact that there are two separate epochs covered by our datasets, each of ~200 days in duration. This allows for the presence of different active groups in each epoch, as well as some common regions in both, creating a flexible evolving map of the active regions. Our selection also matches the range of spot lifetimes from 70 to 350 days described for WASP-52 by Mancini et al. (2017). Several tests showed that those values are not critical in terms of fitting quality and do not produce any bias in other parameters. Furthermore, we did not consider fitting for differential rotation (which StarSim could) as this would add even more complexity to the parameter space. If significant differential rotation was present, this would be naturally accounted for by the fitting algorithm through the appearance and disappearance of active regions at slowly varying longitudes. The temperature contrast between faculae and the photosphere was fixed to ΔT_fc = 50 K. This is consistentwith the results from Solanki (1993) (see Herrero et al. 2016, for further discussion). As is shown in Sect. 3.3.2, our inversion process yields solutions with insignificant facular coverage independently of the temperature contrast ΔT_fc. Finally, we set a restriction to the active element colatitudes, namely that we only allowed their presence in one of the stellar hemispheres. This is possible in our case because of latitudinal symmetry. By doing so, we guarantee the absence of spot crossing events whenwe study planetary transits later in our analysis.

As discussed in Sect. 2.3.2, there are some other parameters that need to be set beforehand. One is the number of iterations of the MCSA algorithm. Ideally, this should be sufficiently large to ensure consistency in the resulting surface map but at the same time be computationally affordable. The number of iterations is also linked to the quantity of active elements considered on the stellar surface, since each element contributes 5 potential parameters. This number has to be sufficient large to generate the inhomogeneities of the stellar surface causing the photometric variability, yet not as large as to overfitinstrumental noise. It partially depends on the timespan of the dataset and the lifetime of spots, and the maximum radius of the active elements, which we set to 10 degrees to be consistent with the small spot approximation.

We ran a battery of tests considering different number of iterations and surface active elements. An initial exploration of the parameter space was used to obtain valid guesses of the optimal parameters P_rot, ΔT_sp, Q, and s (jitter). We systematically explored combinations of number of iterations, from 500 to 10 000, and number of spots from 60 to 150 and evaluated the resulting values. For each of the explored pairs, we ran 112 realisations starting with different random spot maps, except for the runs with 10 000 iterations, for which we considered 56 realisations.

The results of the tests are shown in Fig. 5. It is not surprising to see that the quality of the fits improves with the number of iterations, as it also does with the number of spots. The criterion to select an adequate number of iterations is mostly related to computational cost. For the present problem, numbers above 3000 iterations per MCSA step are prohibitive. Nevertheless, 3000 iterations already delivers very consistent solutions from random starting points. The final variance of the statistic from the sample converged solutions is about 5, which is sufficient to guarantee stable solutions. For lower numbers of iterations, this number increases to 7 (1500) and 17 (500). We note that for 10 000 iterations (3× longer computational time) the final variance is approximately 3. Regarding the number of active elements we see that less than 80 spots is clearly insufficient to fit our data (precision & time span); and we also see that the improvement in likelihood flattens out quite apparently beyond 120 spots. For simplicity arguments, and factoring in again computational costs, we find that using 100 spots delivers sufficiently reliable fits, and note as well that the average values for 150 and 100 spots overlap at the 1-sigma level. Thus, our adopted values regarding the number of MCSA iterations and spots was 3000 and 100, respectively, and this produces a optimal solutions with a variance of .

Fig. 5 Trial tests showing the statistic as a function of the number of spots and for several values of the number of iterations of the MCSA algorithm. The size of the error bars corresponds to the 1-sigma intervals around the optimal solutions found from 112 random initial spot maps, except for the case with 10 000 iterations, where we employed 56.

3.3.2 Fit results

The available datasets consist of multiband photometry and we are especially interested in exploring the parameters related to the chromatic properties of active regions, which play a major role on the characterisation of the effects of activity on transit spectroscopy measurements studied in Sect. 4. The parameters describing stellar properties must be fitted simultaneously to provide consistent solutions, as explained in Sect. 2.3.3. In our case, besides the 100-spot map, those parameters are the rotation period, P_rot, which drives the timescale of the variability, the facula-to-spot area ratio, Q, which determines the fractional coverage of bright active elements, the temperature contrast of the spots, ΔT_sp = T_ph − T_sp, and the additional noise term or jitter, s.

As discussed in Sect. 2.3.2 and shown in Fig. 3, the inversion procedure with StarSim consists of firstly selecting a parameter set drawn from a prior distribution covering the relevant parameter space, and subsequently optimizing random spot maps for each realisation. The resulting statistic J is then used to evaluate the relative quality of each set of parameters. We have shown before that the intrinsic variance of the statistic for our problem is , which indicates that differences of such order for different parameter sets are statistically indistinguishable. Therefore, the recipe that we adopt is taking the best (larger J) solution from the batch and then considering as statistically equivalent those that are within 3-sigma (). We are aware that this is not completely satisfactory because the interval may include some solutions with midly sub-optimal parameter sets, but we adopt this criterion to ensure a statistically significant collection of good solutions at the expense of some increase in the parameter uncertainties. Increasing the number of MCSA iterations or speeding up convergence could make this criterion more stringent. Admittedly, this is a limitation of the current algorithm imposed by the complexity of the optimisation method that should be improved in subsequent releases of StarSim.

After initial tests, for the case of WASP-52 we adopted uniform priors on the parameter sets within the ranges given in Table 3. The priors are quite narrow to avoid unnecessarily exploring irrelevant parameter space but, as seen below, the intervals are sampling regions beyond 3-sigma from the optimal parameters. A total of 21 296 realisations were performed and 108 satisfied the likelihood criterion from the highest likelihood one. In a second step, we further narrowed down the parameter priors to include only the ranges defined by the 108 good solutions. An additional 2576 realisations yielded another 257 good solutions, leading to a grand total of 365. These are the θ-configurations that were used to build the parameter distributions. The distributions are plotted in Fig. 6, and the relevant optimal parameters and their uncertainties are listed in Table 3, compared with other values from the literature.

Our model favours a heterogeneous surface dominated by dark spots with a temperature contrast of 575 ± 150 K with respect to the surrounding photosphere. This value is intermediate to those found in the literature, which can separated into two groups, one with low temperature contrasts of ~ 250 K and another one with contrasts ≳1000 K. The temperature contrast we find is generally lower than found by Andersen & Korhonen (2015) for K dwarfs and also lower than the prediction by the relationship of Berdyugina (2005), which would give ~ 1300 K for a K2 V star. This is, however, not surprising since the strong degeneracy between spot temperature contrast and filling factor can severely bias some determinations. We believe that our value, based on multi-colour photometry, provides a robust estimate. As already mentioned, the results from our analysis do not support a significant presence of bright faculae as expected for rapidly-rotating young K dwarfs (Radick et al. 1983; Lockwood et al. 2007). The 365 best solutions indicate a unilateral distribution consistent with Q = 0. We furthermore find a rotation period of days, which is in marginal agreement with the values found by Louden et al. (2017) and Bruno et al. (2020).

Fig. 6 Histograms and correlations of 365 statistically-equivalent solutions for the fitted StarSim parameters: the rotational period, Prot, the facula-to-spot area ratio, Q, the spot-photosphere temperature contrast, ΔTsp, and the photometric jitter, s. The dotted lines mark the median and 68% percentiles of the parameter distributions and the blue lines indicate the parameters of the solution with the best likelihood value found.

3.3.3 The evolving surface of WASP-52

In addition to the optimal parameters and their uncertainties, our inversion realisations also provide a picture of the stellar surface as a function of time. From the 365 accepted solutions we selected only those satisfying with respect to the best likelihood value. This yields 21 maps, which should be a good representation of the optimal model considering both the intrinsic randomness of the map inversion from the MSCA algorithm and the statistical variance of the stellar parameters. Figure 7 shows the longitudinal spot filling factor projected on the stellar equator, averaged for the 21 optimal maps, as a function of time. We perform a latitudinal projection because of the degeneracy present (see Sect. 2.3.4). The longitudes in the map are plotted in the reference frame of the rotational period found in our analysis (see Table 3). The map suggests that there is a clear dominant active region at a longitude of ~30 deg at the beginning of the first epoch, which changes into a more complex longitudinal pattern later in the season, and finally the spottedness level suffers a suddenly decrease. Active regions typically last for ~8 rotations or ~140 days. The second epoch appears to show a better-defined dominance of spot regions, with two alternating active longitudes at about approximately 130 deg and −50 deg, that is,in almost perfect opposition. The distribution does not indicate signs of differential rotation along the sampled period of time. The existence of measurable differential rotation would be seen in Fig. 7through the presence of two or more active regions with drifting longitudes. We do not see such effect in the time-series of WASP-52 thus probably indicating that either differential rotation is negligible or that dominant active regions appear at similar latitudes. Such dominant regions are used by the fitting algorithm to define the rotation period.

The best fits to the observational multiband photometry using the modelled parameters of the star and their associated spot maps are shown in Fig. 8, together with the multiband photometric measurements. The solid lines show the average models resulting from the optimal 21 surface maps, and the shaded bands show the 1σ range. The two epochs as shown in the top and bottom panel are separated by a ~140-day gap. The grey curves at the bottom of each panel show the projected spot filling factor, also showing the average model and the 1σ range. During our observational campaign, total filling factor values covered a range ~ 3–14% along an observational timespan of ~600 days in 2016–2017. This implies flux modulations of ~ 4–7% due to time-varying spot coverage.

Fig. 7 Longitudinal spot filling factor projected on the stellar equator of WASP-52 as a function of time, covering photometricepochs in 2016 and 2017 and corresponding to the average of 21 optimal solutions. Only stellar longitudes visible at each time are shown (hence the band structure). The colour scale indicates the fractional latitude-projected spot coverage for each 15-degree longitude bin.

4 Chromatic spot effects on simulated transits

4.1 Transit depth variability due to spots

For magnetically-active stars showing spots on their surface, the transit depth caused by a planet passing in front of the host star depends both on the planet effective radius and the possible inhomogeneities on the stellar photosphere. Using our nomenclature in Sect. 2.2, it is straightforward to show that a planet producing a depth when transiting an immaculate star, produces a transit depth when crossing a spotted star, where (11)

with δ_sp being the total spot filling factor. Such very simple model illustrates the chromaticity of the problem, since the function Φ is wavelength dependent. This has been identified as a source of systematic effects for transmission spectroscopy. For this reason, studying the imprint of active regions is of crucial importance for future instruments and projects aiming at the study of exoatmospheres such as, for example, the JWST and Ariel space missions.

We can use StarSim to estimate the impact of spots on the transit depth of WASP-52, by simulating transits of the planet. According to the ephemeris of WASP-52, a total of 112 and 104 transits occurred during the time span covered by the first and second seasons of our light curves, respectively. We used the StarSim model of stellar activity obtained in Sect. 3.3 for WASP-52. It is important to emphasise that we solely focus on the effects of unocculted star spots, whose properties cannot be constrained from the transit itself. Therefore, in our simulations we avoided planet-spot crossings by positioning active regions outside the swath covered by the planet during the transit. Occulted spots during transits may also have a significant impact on the planet properties. However, if photometry of sufficient precision and time cadence is available (such as expected for JWST and Ariel), spot-crossing events could be identified from the transit observations themselves and potentially modelled. Unocculted spots, on the contrary, leave no visible signature.

We define the transit depth as the difference between the out-of-transit and the in-transit (at transit mid-time) observed flux. The out-of-transit is the flux of the model without transit calculated at the mid-transit time. Hence, both stellar spots and limb darkening effects cause variations that depend on wavelength. We take this approach for computational simplicity, because we only need toevaluate the flux at the transit mid point and not the full transit event. The transit depth can then be written as (12)

where, SP(λ) and LD(λ) are the additive chromatic effects induced by spots and by stellar limb darkening, respectively, and is the non-chromatic (asymptotic) transit depth, which we assume to be constant, that is, we do not consider here the planet atmosphere.This can be written as (13)

where r and R_* are the planetary and stellar radii, respectively. In this paper, we assume for WASP-52 b (Öztürk & Erdem 2019).

The transit depth including chromaticity effects can then be expressed as (14)

where Δr accounts for the effect of spots and limb darkening on the planetary radius as a function of wavelength. Expanding Eq. (12), and neglecting second order terms, we derive (15)

Therefore, we can separate the chromatic contributions due to limb darkening (Δr_LD) and spots (Δr_SP) as, (16)

and (17)

Figure 9 shows the results of the simulations. The shaded gray region in panel a displays the range of transit depth values as a function of wavelength for the simulated transits comprised within the timespan of our datasets. Such transits differ in the amount and distribution of spots in the projected stellar surface. As a reference, the red line shows the transit depth in the case of an immaculate star. In this latter case, , the wavelength dependence is solely due to the limb darkening effect, which for WASP-52 produces transit depths from 0.0285 to 0.0256 in flux ratio units from the visual to the infrared bands, respectively. As reference, the dot-dashed blue line indicates the expected transit depth for a uniform planetary disc corresponding to . The shaded grey region in panel (b) shows the same results when removing the contribution of limb darkening, , therefore only including the chromatic effect due to the presence of spots. The results indicate that transit depth variations of WASP-52 b due to spots vary from ~2 × 10⁻⁴ to ~3 × 10⁻³ in the visual region, and from ~10⁻⁴ to ~7 × 10⁻⁴ in the NIR region. The range in the variations is caused by both different spot filling factors during each transit, and the uncertainty on the exactdistribution of spots. In other words, our results show that the transit depth of WASP-52 b can be overestimated by up to ~12 and ~3% in the visualand NIR bands, respectively, thus potentially affecting the retrieval of planetary atmosphere parameters, which could be one or two orders of magnitude smaller (Kreidberg 2018).

To understand the impact of different spot configurations we selected two transit events representing low and high-activity levels at well-sampled epochs of the light curves. These are labelled in Fig. 8 as TR-1 at BJD = 2 457 739.475 (16 Dec. 2016, low-activity case), when the projected filling factor of spots is estimated to be close to its minimum value around 5%, and as TR-2 at BJD = 2 457 954.698 (20 Jul. 2017, high-activity case), when spot projected filling factor reaches ~12%. Panel b of Fig. 9 shows the transit depth variation of TR-1 and TR-2 in blue and red, respectively, for the 21 optimal maps, with vertical error bars indicating the 1σ variance. These test cases clearly show that the strength of the chromatic effect is ~3 times larger when the projected spot filling factor is larger by about the same factor. This implies that the effect is also correspondingly larger for stars showing higher levels of stellar activity.

Table 4 provides the statistics of the solutions for the TR-1 and TR-2 transit cases for different wavelengths, encompassing the visual channels of the Ariel mission (~550 and ~705 nm), the infrared (~975 nm) and the low resolution spectrographs (NIR, AIRS). We computed the mean transit depth variation induced by spots, ⟨SP (λ)⟩, and its corresponding planetary radius variation ⟨Δr∕R_*⟩, along with the standard deviation σ induced by the uncertainty in the determination of the spot map. The statistics reveal that, even for the relatively low stellar activity level of WASP-52, with a spot filling factor of ~5%, the signature of stellar activity on the transit depth can be close to 10⁻³ of the flux in the visual bands and ~2 × 10⁻⁴ in the NIR bands. This translates into relative planetary radius changes from ~2 × 10⁻³ to ~5 × 10⁻⁴, depending on the wavelength considered, which are of the same order as the typical signature of exoplanet atmospheres on transmission spectra. Furthermore, at higher activity level such as that of TR-2, the effect of the starspots is about twice as large. This illustrates the complication of studying atmospheres of exoplanets orbiting relatively active host stars.

Fig. 8 StarSim model fits to multiband photometric ground-based observational data as described. Solid curves represent the mean of 21 optimal solutions of the inverse problem. The shaded bands indicate the 1σ ranges. The gray line at the bottom is the projected spot filling factor of the maps, also showing the mean and 1σ ranges. TR1and TR2 indicate transit events as discussed in Sect. 4.

Fig. 9 (a) Transit depth as a function of wavelength for the WASP-52 system simulated assuming an immaculate surface (red) and including spots in the photosphere (gray band), avoiding transit spot crossings. The gray area illustrates the total coverage from all simulated transits for the 21 optimal spot configurations, clearly showing the transit depth chromatic bias produced bynon-occulted spots. The dotted blue line is the non-chromatic transit depth assuming a uniform photosphere () according tothe parameters in Öztürk & Erdem (2019). (b) Chromatic signature on transit depths due to spots after subtracting limb-darkening effect (in red in top panel). TR-1 and TR-2 show two example transits corresponding to low- and high-activity – marked with vertical lines in Fig. 8 – respectively, and the gray band represents the area covered by all simulated transits. The black solid line with symbols shows the chromatic effect estimated using stellar spot parameters derived by Bruno et al. (2020).

4.2 Correcting transit depth for stellar activity

Various of methodologies to correct the impact of stellar activity on transmission spectroscopy of transiting planets have been suggested (see e.g. Knutson et al. 2012; Pont et al. 2013; Alam et al. 2018). Typically, they are based on measuring the photometricvariability of the host star, from which the spot filling factor and the temperature contrast or a combination of both are estimated. Furthermore, Rackham et al. (2018) claim that the methodologies applied to calculate these corrections could still be biased as they are usually computed assuming a single spot on the surface of the star. The problem is that photometric variations do not generally provide a realistic representation of the spot filling factorthroughout the rotation cycle. For example, a uniformly densely spotted star could produce a bias on the transit depth while not producing any observable photometric variability. Using semi-empirical relations between thefilling-factor of spots and the photometric variability of stars, Rackham et al. (2018) conclude that the chromatic effect on the transit depth due to spots could be several times larger than the signal of atmospheric features, thus greatly difficulting the general use of transmission spectroscopy for stars with some level of magnetic variability.

We demonstrate here that using multiband photometric light curves contemporaneous to transit observations, we can overcome this problem, at least partly. The advantage of using several photometric bands is that, as discussed in the previous section (see also Mallonn et al. 2018), it is possible to independently estimate the zero point with respect to the unspotted photosphere (z), and thus the absolute spot filling factor (including any persistent level of spottedness during the rotation cycle), and the temperature contrast. This provides extremely valuable information to correct transit depth variations due to spots. We can actually use our StarSim simulations of WASP-52 to investigate at which level of accuracy we can correct the effects of starspots on transit data.

The StarSim simulations of the WASP-52 planetary system discussed in the previous section reveal that, without additional information allowing us to infer the distribution and filling factor of spots during the transit, the planet-to-star radius cannot be determined with an accuracy better than a few percent. This is because in that case, we do not know at which rotation phase of the star the transit is observed, and therefore, we can only estimate a range of possible filling factors based on the amplitude of the photometric modulation. However, if we know the rotation phase at the time of transit, and we can estimate the spot properties and distribution from light curves, we can infer the correction that needs to be applied to the relative radius of the exoplanet (or transit depth) with an accuracy of order 10⁻⁴, as shown in Table 4.

Figure 10 illustrates the sequence of spot corrections on the transit depth as a function of time. We plot the planet radius variation due to unnoculted spots, Δr_SP∕R_*, following our StarSim model describing the photometric light curve in Fig. 8. Each panel corresponds to a different wavelength band of interest for the Ariel mission, as an example. To test all possible filling factors, we consider transits occurring at any time, but actual WASP-52 b transits are indicated as filled dots. The average value of all the simulations performed is displayed as a solid line. It is clear from this figure that the variation of the apparent planet radius depends on the projected spot filling factor and follows the stellar rotation period. The relative effect on the radius diminishes towards longer wavelengths due to the lower flux contrast of spots, from an average of ~4 × 10⁻³ to ~8 × 10⁻⁴ at 550 nm (VIS-Phot) and 6 μm (AIRS-Ch1), respectively. In all cases it is possible to estimate the unocculted spot correction to the planetary relative radius with a precision close to 10%, estimated as the standard deviation (σ) of the sample of 21 optimal fits to the light curves. As explained above, this residual uncertainty is due to the variance of the stellar parameters and the randomness of the determination of the active region map.

It is interesting to quantify the improvement on a realistic transit depth determination that we can reach using this approach. We measurethis improvement through what we call activity attenuation factor, which we compute as the ratio between the effect of spots on the transit depth and the uncertainty of the correction resulting from the StarSim model (⟨SP(λ)⟩/σ_SP). The distribution of attenuation factors for the different simulated transits is shown in Fig. 11, and this is a measure of the ability of the StarSim model to correct starspot effects. The plots indicate that, using simultaneous photometry, we can correct the effects of spots on the relative radius of WASP-52 b by factors between 2 and 15 depending on the spot filling factor and the accuracy and coverage of the contemporaneous photometric monitoring. Lower attenuation factors mainly correspond to transits occurring when the projected filling factor of spots is smaller and having poor photometric coverage. On the other hand, larger factors are mainly due to transits with greater spottedness at times of well-sampled photometriclight curves. On average, we can expect an attenuation factor around 10 on the relative radius measured for WASP-52 b.

The simulations of this planetary system show that, using contemporaneous photometric data, we are able to estimate transit depth corrections due to unnoculted spots with uncertainties of around few times 10⁻⁵ in the NIR, regardless of the fraction of photosphere covered by spots (see error bars in the first two columns of Table 4). And we recall that this corresponds to a star displaying a modulation of about 7% in flux in the visible band, and a spot filling factor of about 3–14%.

Fig. 10 Effect of unocculted dark spots on the planetary radius ΔrSP∕R* over time fordifferent wavelength bands of the Ariel mission. The solid lines correspond to the average value of 21 optimal fits, and the 1σ confidence level is shown as a light-shaded band. Solid dots indicate the timestamps of WASP-52 b transits according to the available ephemeris. Note the different vertical scale for each wavelength panel. Gray dots at the bottom of the panels represent the variance (σ, as a percentage value) of the results calculated from the different maps, as labelled in the right axis. The gap in the middle of each plot corresponds to the period between observational seasons without available photometric data (see Fig. 8).

Fig. 11 Distribution of activity attenuation factors from StarSim modelling considering all transit simulations. This factor is calculated as the ratio between the mean effect of spots on the transit depth and the uncertainty of the correction given by StarSim model (⟨SP(λ)⟩/σSP). Blue and red dots indicate the position in the histograms of the example transits TR-1 and TR-2, respectively.

4.3 Comparison with alternative activity correction methods for WASP-52 b

Several approaches have been presented to account for and correct out the chromatic effects of activity in transit observations. Instead of fitting photometric light curves, Rackham et al. (2019) estimate the filling factor of FGK spectral type stars from the peak-to-valley variability measured from Kepler photometry using the simple analytical form in Eq. (11). The temperature contrast of the active regions and the projected filling factor of spots are strongly degenerate when only a single passband is used to measure photometric variability. To overcome this problem, Rackham et al. (2019) assume a temperature for the spots consistent with empirical values reported by Berdyugina (2005). The authors also assume a uniform distribution of dark spots over the stellar photosphere. If we apply this approach to WASP-52, the temperature contrast of spots would be 1290 K. However, in spite of being significantly larger than the value we find in our fits, with this value it is not possible to reproduce the ~7% peak-to-peak photometric variability of our light curves (see Fig. 2 in Rackham et al. 2019) when assuming uniform spot coverage. Unphysically high spot coverage levels, a much larger spot temperature contrast or non-uniform distributions would be needed to explain the large observed amplitude. In addition to the different temperature contrast we obtain, the surface distribution of active regions derived with StarSim is distinctly non-uniform (see Fig. 7).

Bruno et al. (2020) also estimate the properties of starspots for WASP-52, following the procedure of Huitson et al. (2013) and the light curve normalisation factor suggested by Aigrain et al. (2012). Based on data in Alam et al. (2018), the authors fit both parameters from ASAS-SN and AIT photometry and transit spectroscopy from HST/STIS, WFC3 IR, and Spitzer/IRAC. They obtain temperatures <3000 K for the starspots and a ~5% filling factor. Although the filling factor is consistent with our results, the temperature of spots is much larger than the value we measure from multiband photometry. The chromatic signature produced by the parameters of Bruno et al. (2020) is plotted in panel b of Fig. 9 (black solid line) using Eq. (11). At wavelengths below ~1.5 μm), the resulting values are within the region allowed by our StarSim runs, and below 1 μm the predicted corrections are not far from the TR-2 (~12% spot filling factor) case. The chromatic dependence, however, is significantly different. It is therefore possible that the transmission spectrum of Bruno et al. (2020) is somewhat affected by the different chromatic slope, although the average correction value should not be strongly biased. Note, however, that if we were to use the spot parameters adopted by Bruno et al. (2020) (ΔT_sp = 2200 K, compared with 575 K in our work) at NIR wavelengths, very strong deviations would occur, with spot corrections overestimated by factors of 2–5.

The differences found highlight the importance of extracting the chromatic signature from multiband photometric datasets so that the degeneracy between spot temperature and filling factor can be broken, and to reproduce realistic spot distributions. The determination of physical properties of active regions is a key point for the study of their chromatic signature on transit observations. Besides, using independent data, instead of the same transmission spectroscopy used for exoplanet characterisation, prevents misinterpreting atmospheric features.

We caution that the procedure employed in the literature above corrects the chromatic effects of active regions in transits by estimating the spot filling factor and temperature. However, the actual spot distribution on the stellar surface is still relevant because of limb darkening effects. Figure 12 shows a comparison between two ad-hoc spot maps (MAP-1 and MAP-2), a fitted map of WASP-52 at epoch BJD = 2 457 645 and the theoretical depth variation with homogeneous distribution (Eq. (11)). All these models have been calculated with the same projected spot coverage of 5%. MAP-1 has a single spot at the centre of the disc while MAP-2 has spot group close to the limb. The differences are quite significant, especially in the visible. This simple exercise illustrates the need to consider the full stellar surface at the time of transit to properly correct for spot effects.

Fig. 12 Chromatic effects due to spot distributions with identical filling factors. MAP-1 and MAP-2 are two ad-hoc generated spot maps to show the differences in their depth profiles and the comparison with that provided with Eq. (11), shown in red. The green curve is a fitted map for WASP-52 at epoch BJD = 2 457 644.95 with the parametric distribution error bands.

5 Conclusions

In this study we present the StarSim code, which is able to model the surface of active stars through the inversion of photometric time-series. With the resulting model, we can investigate the influence of the chromatic effect produced by activity features for transmission spectroscopy. As a test case, we use simultaneous BV RI-filter photometry of the young active planet-host star WASP-52 covering ~600 days of 2016 and 2017 from the STELLA and TJO telescopes. The data are used to constrain an activity model with a heterogeneous surface composed of dark spots with a temperature difference of 575 ± 150 K with respect the surrounding photosphere. Our modelling rejects the hypothesis of the presence of bright facular regions, as expected for young FGK stars. We determine the rotation period of the star to be days. The fit also yields a probabilistic map of active regions, showing two or three prominent spot complexes at stable stellar longitudes.

It is important to emphasise that the ultimate goal of our effort is to correct the transit depth as a function of wavelength, and not necessarily to obtain precise reconstructions of the stellar surface features. That is, although there might be degeneracies in the number of spots and their exact distribution, it is the overall photometric effect at the time of transit what really matters. To study the influence of the derived activity model on the observables of transmission spectroscopy, we produced simulated transits of WASP-52 b by avoiding spot-crossing events. The chromatic effect of the spot map results in uncertainties of ~10% in transit depth (SP(λ)/) and of up to ~5% in the planetary radius (Δr_SP∕r), at visible wavelengths. After correcting for spot effect using the StarSim model, we are able to reduce residual depth uncertainties down to ~10⁻⁴ at 550 nm (Ariel/VIS-Phot) and ~ 3 × 10⁻⁵ at 6 μm (Ariel/AIRS-Ch1).

The remaining uncertainty of the correction factor, computed as the 1σ interval of the depth variation of the simulated transits for the fitted spot maps, encloses several effects, most importantly: 1) the photometric precision and phase coverage of the monitoring along time, 2) the incompleteness of the model implemented and, 3) the uncertainties in retrieving optimal surface maps and stellar parameters in the context of the photometric inverse problem.

The main advantage of the approach we present to compute the chromatic effect on the depth of exoplanet transits is that it is based on an independent and consistent deterministic method allowing to accurately determine the stellar parameters, the filling factor and the distribution of spots. Thus, the effect of spots and their different positions on the stellar disc are also taken into account. We have only considered here the effect of unocculted spots, which increase the transit depth. On the other hand, occulted spots produces the opposite effect. In our previous work on GJ 1214, we described the difficulty to achieve a unique solution for the planetary optical spectral slope when the transit observation is affected by a mixture of non-occulted and occulted spots, even if multi-colour monitoring of the star is available (Mallonn et al. 2018). However, very-high precision transit photometry could help to detect the effects of occulted spots and to derive their properties.

With the expected launch of JWST in 2021 and Ariel mission in 2028, the ability to observe transmission spectra of transiting planets will increase both in the number of stars and precision. In this context, our results show that contemporaneous ground-based photometric monitoring will be crucial to adequately model spot-driven stellar activity variations and correct out their chromatic effects on the transit depth measurements. In essence, our method exploits the chromaticity of spot-induced stellar time variability. Activity correction may be essential to attain the required accuracy in transmission spectroscopy so that meaningful constraints on the planetary atmospheres can be defined. We note that the correction is necessary irrespective of the precision of the space-based photometry since it deals with a systematic effect of astrophysical origin. We expect that higher activity attenuation factors could be achieved by optimizing the strategy of the photometric observations, that is, making sure that multi-colour measurements are obtained before and after the transits to better constrain the model. Photometricobservations from space would even boost the accuracy in light curve modelling, and hence, the correction capabilities.

Acknowledgements

We are grateful to the referee for their valuable and insightful comments and suggestions that helped to improve the paper significantly. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme.

References

All Tables

All Figures

	Fig. 1 Example of a stellar photosphere including modulating and non-modulating spots. The model is composed of a polar cap, which does not contribute to light curve modulation, with a projected filling factor δ0, and an equatorial small spot with a maximum projected filling factor δM. F(⋅) is the brightness of each surface element.
In the text

Fig. 2 Graphical representation of Eq. (3) for the parameters of interest assuming a stellar photospheric temperature of Tph = 5000 K and a black-body for the spectral energy distribution of each surface element. (a) Photometric amplitude as a function of the modulating filling factor for different spot temperature contrasts. A non-modulating filling factor δ0 = 0.05 is employed, but adopting other values has negligible effect at this scale. The gray area covers the spectral range from 400 nm to 1 μm. (b) Difference in activity-induced amplitude for photometric bands BV R with respect to the I band as a function of spot temperature contrast. A value of 0.05 is adopted as modulating filling factor (δM) and examples for non-modulating filling factors of δ0 = 0 and 0.1 are shown. The gray lines illustrate how the bi-valuate nature of the temperature contrast effect can be resolved if more than two bands are used. (c) Difference in activity-induced amplitude for photometric bands BV R with respect to the I band as a function of the non-modulating filling factor. A modulating filling factor δM = 0.05 and a spot temperature contrast ΔTsp = 600 K are adopted.

In the text

	Fig. 3 Flow chart of the parameter retrieval approach as presented in Sect. 2.3.
In the text

Fig. 4 Left panel: StarSim model fits to multiband synthetic light curves assuming the parameters in Table 1, for a subset of 50 days. Solid curves represent the mean of ~20 optimal solutions of the inverse problem. The gray line at the bottom is the projected spot filling factor of the maps, also showing the mean and 1σ shaded band. The error intervals of the photometry bands are very small and difficult to recognise. Right panels: longitudinal spot filling factor projected on the stellar equator as a function of time, for both the input (bottom) and the retrieved (top) spot maps. Only stellar longitudes visible at each time are shown (hence the band structure). The colour scale indicates the fractional latitude-projected spot coverage for each 1-degree longitude bin.

In the text

	Fig. 5 Trial tests showing the statistic as a function of the number of spots and for several values of the number of iterations of the MCSA algorithm. The size of the error bars corresponds to the 1-sigma intervals around the optimal solutions found from 112 random initial spot maps, except for the case with 10 000 iterations, where we employed 56.
In the text

Fig. 6 Histograms and correlations of 365 statistically-equivalent solutions for the fitted StarSim parameters: the rotational period, Prot, the facula-to-spot area ratio, Q, the spot-photosphere temperature contrast, ΔTsp, and the photometric jitter, s. The dotted lines mark the median and 68% percentiles of the parameter distributions and the blue lines indicate the parameters of the solution with the best likelihood value found.

In the text

Fig. 7 Longitudinal spot filling factor projected on the stellar equator of WASP-52 as a function of time, covering photometricepochs in 2016 and 2017 and corresponding to the average of 21 optimal solutions. Only stellar longitudes visible at each time are shown (hence the band structure). The colour scale indicates the fractional latitude-projected spot coverage for each 15-degree longitude bin.

In the text

	Fig. 8 StarSim model fits to multiband photometric ground-based observational data as described. Solid curves represent the mean of 21 optimal solutions of the inverse problem. The shaded bands indicate the 1σ ranges. The gray line at the bottom is the projected spot filling factor of the maps, also showing the mean and 1σ ranges. TR1and TR2 indicate transit events as discussed in Sect. 4.
In the text

Fig. 9 (a) Transit depth as a function of wavelength for the WASP-52 system simulated assuming an immaculate surface (red) and including spots in the photosphere (gray band), avoiding transit spot crossings. The gray area illustrates the total coverage from all simulated transits for the 21 optimal spot configurations, clearly showing the transit depth chromatic bias produced bynon-occulted spots. The dotted blue line is the non-chromatic transit depth assuming a uniform photosphere () according tothe parameters in Öztürk & Erdem (2019). (b) Chromatic signature on transit depths due to spots after subtracting limb-darkening effect (in red in top panel). TR-1 and TR-2 show two example transits corresponding to low- and high-activity – marked with vertical lines in Fig. 8 – respectively, and the gray band represents the area covered by all simulated transits. The black solid line with symbols shows the chromatic effect estimated using stellar spot parameters derived by Bruno et al. (2020).

In the text

Fig. 10 Effect of unocculted dark spots on the planetary radius ΔrSP∕R* over time fordifferent wavelength bands of the Ariel mission. The solid lines correspond to the average value of 21 optimal fits, and the 1σ confidence level is shown as a light-shaded band. Solid dots indicate the timestamps of WASP-52 b transits according to the available ephemeris. Note the different vertical scale for each wavelength panel. Gray dots at the bottom of the panels represent the variance (σ, as a percentage value) of the results calculated from the different maps, as labelled in the right axis. The gap in the middle of each plot corresponds to the period between observational seasons without available photometric data (see Fig. 8).

In the text

	Fig. 11 Distribution of activity attenuation factors from StarSim modelling considering all transit simulations. This factor is calculated as the ratio between the mean effect of spots on the transit depth and the uncertainty of the correction given by StarSim model (⟨SP(λ)⟩/σSP). Blue and red dots indicate the position in the histograms of the example transits TR-1 and TR-2, respectively.
In the text

	Fig. 12 Chromatic effects due to spot distributions with identical filling factors. MAP-1 and MAP-2 are two ad-hoc generated spot maps to show the differences in their depth profiles and the comparison with that provided with Eq. (11), shown in red. The green curve is a fitted map for WASP-52 at epoch BJD = 2 457 644.95 with the parametric distribution error bands.
In the text