Issue 
A&A
Volume 641, September 2020



Article Number  A82  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202037586  
Published online  11 September 2020 
Correcting for chromatic stellar activity effects in transits with multiband photometric monitoring: application to WASP52
^{1}
Institut de Ciències de l’Espai (ICE, CSIC),
Campus UAB, Can Magrans s/n,
08193
Bellaterra,
Spain
email: rosich@ice.cat
^{2}
Institut d’Estudis Espacials de Catalunya (IEEC),
08034
Barcelona, Spain
^{3}
LeibnizInstitut für Astrophysik Potsdam,
An der Sternwarte 16,
14482
Potsdam, Germany
Received:
27
January
2020
Accepted:
29
June
2020
Context. The properties, distribution, and evolution of inhomogeneities on the surface of active stars, such as dark spots and bright faculae, significantly influence the determination of the parameters of an orbiting exoplanet. The chromatic effect they have on transmission spectroscopy, for example, could affect the analysis of data from future space missions such as James Webb Space Telescope and Ariel.
Aims. To quantify and mitigate the effects of those surface phenomena, we developed a modelling approach to derive the surface distribution and properties of active regions by modelling simultaneous multiwavelength timeseries observables.
Methods. We present an upgraded version of the StarSim code, now featuring the capability to solve the inverse problem and derive the properties of the stars and their active regions by modelling timeseries data. As a test case, we analyse ~600 days of BV RI multiband photometry from the 0.8m Joan Oró (TJO) and 1.2m STELLA telescopes of the K2 V exoplanet host star WASP52. From the results, we further simulated the chromatic contribution of surface phenomena on the observables of its transiting planet.
Results. Using StarSim we are able to determine the relevant activity parameters of WASP52 and reconstruct the timeevolving longitudinal map of active regions. The star shows a heterogeneous surface composed of dark spots with a mean temperature of 575 ± 150 K lower than the photospheric value, with filling factors ranging from 3 to 14%. We used the results to study the chromatic effects on the depths of exoplanet transits obtained at different epochs and corresponding to different stellar spot distributions. In the case of WASP52, which has peaktopeak photometric variations of ~7% in the visible, the residual effects of dark spots on the measured transit depth of its giant planet, after applying the calculated corrections, are about 10^{−4} at 550 nm and 3 × 10^{−5} at 6 μm.
Conclusions. We demonstrate that by using contemporaneous groundbased multiband photometry of an active star, it is possible to reconstruct the parameters and distribution of active regions over time, thus making it feasible to quantify the chromatic effects on the planetary radii measured with transit spectroscopy and mitigate them by about an order of magnitude.
Key words: techniques: photometric / methods: data analysis / methods: statistical / stars: individual: WASP52
© ESO 2020
1 Introduction
After the groundbreaking 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; AngladaEscudé 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 lowdensity planets with atmospheres of low mean molecular weight (hydrogendominated) 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^{−3}, which makes their detection very challenging, particularly in the case of active host stars, where other sources of variations are expected.
Future spacebased 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, largescale survey of the chemical compositions and thermal structures of the atmospheres of ~1000 known transiting exoplanets in the optical and nearinfrared (NIR) wavelengths (0.5–8 μm) following itsexpected launch in 2028 (Tinetti et al. 2018; Encrenaz et al. 2018).
Magneticallydriven 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 Rayleighlike slope and grey opacity, respectively (Pinhas & Madhusudhan 2017). Distinguishing the different contributions to wavelengthdependent 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 timeseries 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 crosscorrelation 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 timeseries 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 WASP52, as shown in Sect. 3. In Sect. 4, we use the retrieved activity model for WASP52 to study the different sources of chromatic effects on the transit depths of its known planet WASP52 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 timeseries data of magneticallyactive 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 BTSettl 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 magneticallyactive 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 faculatospot 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 wavelengthdependent 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 timeseries are obtained by integrating over all the grid elements at each time step, taking into account the rotation of the star and the timeevolution of the active regions. When simulating photometrictimeseries, 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 zeropoint 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 userdefined mask of spectral lines can be used to calculate the CCFs. In a first step of each timeseries 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), δ_{0}. 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 δ_{0}, 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_{1} − f_{2})∕f_{1}, which is equivalent to the measurement provided by observational light curves. Rearranging Eqs. (1) and (2), the photometric amplitude can be written as, (3)
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 nonmodulating filling factor δ_{0} 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 ≤ δ_{0} ≤ 1, 0 ≤ δ_{M} ≤ 1, 0 ≤ δ_{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 nonmodulating 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 δ_{0} arises. For typical spot temperature contrasts and visible bands, it can be shown that Φ is significantly greater than δ_{0} 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 bivaluate, 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 nonmodulating filling factor from multicolour photometry. Nevertheless, the scale of the variations makes reliable estimates of δ_{0} very challenging for typical groundbased 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 nonmodulating 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 timeseries data X. This type of nonlinear problem can be expressed as (5)
where X is the timeseries 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 timeseries 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 T_{ph} = 5000 K and a blackbody 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 nonmodulating 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 activityinduced 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 nonmodulating filling factors of δ_{0} = 0 and 0.1 are shown. The gray lines illustrate how the bivaluate nature of the temperature contrast effect can be resolved if more than two bands are used. (c) Difference in activityinduced amplitude for photometric bands BV R with respect to the I band as a function of the nonmodulating filling factor. A modulating filling factor δ_{M} = 0.05 and a spot temperature contrast ΔT_{sp} = 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 timeseries data defined as (7)
where is the loglikelihood 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 noncorrelated Gaussian uncertainties (white noise). The likelihood function is then written as (8)
where is the StarSimgenerated model of the jth timeseries 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 quadratureadded 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 faculatospot 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^{−1} (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 timeseries data.
The optimisation process follows successive iterations where the algorithm randomly selects one of the adjustable parameters from a randomlyselected active element and modifies it slightly. Then, with this new map, the forward problem is recomputed and is evaluatedfor each timeseries 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 timeseries 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 timeevolving 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 nonlinear 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 twostep 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 randomlygenerated 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 statisticallyequivalent 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 equatoron 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 sizelatitude 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 equatoron, 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 equatoron 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 timeseries 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 temperaturesize 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 spotfacula degeneracy can also be broken when contemporaneous photometry and highprecision 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 nonmodulating 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 noncorrelated 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^{−1}. 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 10spot 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 nonmodulating 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 latitudeprojected 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.
Results of inversion tests using synthetic data with StarSim.
3 An activity model for WASP52
As an example of the use of StarSim for the inverse problem, we use multiband photometric timeseries data of the active planet host star WASP52 to retrieve an optimal set of stellar parameters and reconstruct a timevariable probability map of the distribution of active regions.
3.1 The WASP52 exoplanetary system
WASP52 isa young and active K2 V star hosting an inflated Jupitersized 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 RossiterMcLaughlin 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 sodiumbearing atmosphere. In many of those studies, intransit 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 groundbased photometry using a quasiperiodic 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 WASP52b 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 latitudeprojected spot coverage for each 1degree longitude bin. 
Important parameters of WASP 52 (top) and WASP 52b (bottom).
Priors and results from the StarSim inversion of WASP52 photometry.
3.2 Photometric observations
Photometric light curves obtained from two different groundbased observatories, the 1.2m STELLA telescope at Izaña Observatoryin Tenerife and the 0.8m Joan Oró telescope (TJO) at the Montsec Astronomical Observatory in Catalonia, were used in this work. WASP52 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 widefield 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 backilluminated 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 ESOMIDAS routines used for previous monitoring programmes of exoplanet host stars with STELLA/WiFSIP (Mallonn et al. 2015, 2018). Bias and flatfield 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 nonphotometric 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 groundbased 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 nighttonight calibration errors. This strategy is appropriate because we are interested in rotational modulation timescales, with the fewhour time domain being irrelevant.
The TJO telescope provided photometric data using its main imager MEIA2. The instrument has a single 2k × 2k backilluminated thinned CCD yielding a field of view of 12.3 × 12.3 arcmin and a resolution of 0.36 arcsec pixel^{−1}. 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 JohnsonCousins 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 spinorbit angle. The lifetime of the active elements was assumed to be Gaussiandistributed 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 WASP52 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 1sigma 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 1sigma 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 100spot map, those parameters are the rotation period, P_{rot}, which drives the timescale of the variability, the faculatospot 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 3sigma (). We are aware that this is not completely satisfactory because the interval may include some solutions with midly suboptimal 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 WASP52 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 3sigma 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 multicolour 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 rapidlyrotating 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 statisticallyequivalent solutions for the fitted StarSim parameters: the rotational period, P_{rot}, the faculatospot area ratio, Q, the spotphotosphere temperature contrast, ΔT_{sp}, 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 WASP52
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 betterdefined 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 timeseries of WASP52 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 ~140day 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 timevarying spot coverage.
Fig. 7 Longitudinal spot filling factor projected on the stellar equator of WASP52 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 latitudeprojected spot coverage for each 15degree longitude bin. 
4 Chromatic spot effects on simulated transits
4.1 Transit depth variability due to spots
For magneticallyactive 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 WASP52, by simulating transits of the planet. According to the ephemeris of WASP52, 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 WASP52. 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 planetspot 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), spotcrossing 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 outoftransit and the intransit (at transit midtime) observed flux. The outoftransit is the flux of the model without transit calculated at the midtransit 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 nonchromatic (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 WASP52 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)
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 WASP52 produces transit depths from 0.0285 to 0.0256 in flux ratio units from the visual to the infrared bands, respectively. As reference, the dotdashed 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 WASP52 b due to spots vary from ~2 × 10^{−4} to ~3 × 10^{−3} in the visual region, and from ~10^{−4} to ~7 × 10^{−4} 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 WASP52 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 highactivity levels at wellsampled epochs of the light curves. These are labelled in Fig. 8 as TR1 at BJD = 2 457 739.475 (16 Dec. 2016, lowactivity case), when the projected filling factor of spots is estimated to be close to its minimum value around 5%, and as TR2 at BJD = 2 457 954.698 (20 Jul. 2017, highactivity case), when spot projected filling factor reaches ~12%. Panel b of Fig. 9 shows the transit depth variation of TR1 and TR2 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 TR1 and TR2 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 WASP52, with a spot filling factor of ~5%, the signature of stellar activity on the transit depth can be close to 10^{−3} of the flux in the visual bands and ~2 × 10^{−4} in the NIR bands. This translates into relative planetary radius changes from ~2 × 10^{−3} to ~5 × 10^{−4}, 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 TR2, 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 groundbased 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 WASP52 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 bynonocculted spots. The dotted blue line is the nonchromatic 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 limbdarkening effect (in red in top panel). TR1 and TR2 show two example transits corresponding to low and highactivity – 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). 
Mean values of transit depth variations due to spots, , on WASP52 simulated transits.
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 semiempirical relations between thefillingfactor 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 WASP52 to investigate at which level of accuracy we can correct the effects of starspots on transit data.
The StarSim simulations of the WASP52 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 planettostar 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^{−4}, 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 WASP52 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^{−3} to ~8 × 10^{−4} at 550 nm (VISPhot) and 6 μm (AIRSCh1), 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 WASP52 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 wellsampled photometriclight curves. On average, we can expect an attenuation factor around 10 on the relative radius measured for WASP52 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^{−5} 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 Δr_{SP}∕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 lightshaded band. Solid dots indicate the timestamps of WASP52 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 TR1 and TR2, respectively. 
4.3 Comparison with alternative activity correction methods for WASP52 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 peaktovalley 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 WASP52, 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% peaktopeak 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 nonuniform 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 nonuniform (see Fig. 7).
Bruno et al. (2020) also estimate the properties of starspots for WASP52, 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 ASASSN 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 TR2 (~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 adhoc spot maps (MAP1 and MAP2), a fitted map of WASP52 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%. MAP1 has a single spot at the centre of the disc while MAP2 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. MAP1 and MAP2 are two adhoc 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 WASP52 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 timeseries. 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 RIfilter photometry of the young active planethost star WASP52 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 WASP52 b by avoiding spotcrossing 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^{−4} at 550 nm (Ariel/VISPhot) and ~ 3 × 10^{−5} at 6 μm (Ariel/AIRSCh1).
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 nonocculted and occulted spots, even if multicolour monitoring of the star is available (Mallonn et al. 2018). However, veryhigh 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 groundbased photometric monitoring will be crucial to adequately model spotdriven stellar activity variations and correct out their chromatic effects on the transit depth measurements. In essence, our method exploits the chromaticity of spotinduced 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 spacebased 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 multicolour 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 PGC2018098153BC33, as well as the support of the Generalitat de Catalunya/CERCA programme.
References
 Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147 [NASA ADS] [CrossRef] [Google Scholar]
 Alam, M. K., Nikolov, N., LópezMorales, M., et al. 2018, AJ, 156, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Allard, F., Homeier, D., & Freytag, B. 2013, Mem. Soc. Astron. It., 84, 1053 [NASA ADS] [Google Scholar]
 Andersen, J. M., & Korhonen, H. 2015, MNRAS, 448, 3053 [NASA ADS] [CrossRef] [Google Scholar]
 AngladaEscudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Barros, S. C. C., Boué, G., Gibson, N. P., et al. 2013, MNRAS, 430, 3032 [NASA ADS] [CrossRef] [Google Scholar]
 Berdyugina, S. V. 2005, Liv. Rev. Sol. Phys., 2, 8 [Google Scholar]
 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boisse, I., Bonfils, X., & Santos, N. C. 2012, A&A, 545, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bruno, G., Lewis, N. K., Stevenson, K. B., et al. 2018, AJ, 156, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Bruno, G., Lewis, N. K., Alam, M. K., et al. 2020, MNRAS, 491, 5361 [NASA ADS] [Google Scholar]
 Burrows, A. S. 2014, Nature, 513, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45 [Google Scholar]
 Chen, G., Pallé, E., Nortmann, L., et al. 2017, A&A, 600, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Colomé, J., & Ribas, I. 2006, IAU Special Session, 6, 11 [Google Scholar]
 Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Encrenaz, T., Tinetti, G., & Coustenis, A. 2018, Exp. Astron., 46, 31 [CrossRef] [Google Scholar]
 Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Grill, J.B., Valko, M., & Munos, R. 2015, in Neural Information Processing Systems, Montréal, Canada [Google Scholar]
 Hébrard, G., Collier Cameron, A., Brown, D. J. A., et al. 2013, A&A, 549, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Herrero, E., Ribas, I., Jordi, C., et al. 2016, A&A, 586, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252 [NASA ADS] [CrossRef] [Google Scholar]
 Kirk, J., Wheatley, P. J., Louden, T., et al. 2016, MNRAS, 463, 2922 [NASA ADS] [CrossRef] [Google Scholar]
 Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. 1983, Science, 220, 671 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Kreidberg, L. 2018, Exoplanet Atmosphere Measurements from Transmission Spectroscopy and Other Planet Star Combined Light Observations (Berlin: Springer International Publishing AG), 100 [Google Scholar]
 Kurucz, R. L. 2017, Astrophysics Source Code Library [record ascl:1710.017] [Google Scholar]
 Lagrange, A.M., Desort, M., & Meunier, N. 2010, A&A, 512, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Boisse, I., Bouchy, F., Bonomo, A. S., & Moutou, C. 2011, A&A, 533, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lockwood, G. W., Skiff, B. A., Henry, G. W., et al. 2007, ApJS, 171, 260 [NASA ADS] [CrossRef] [Google Scholar]
 Louden, T., Wheatley, P. J., Irwin, P. G. J., Kirk, J., & Skillen, I. 2017, MNRAS, 470, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Mallonn, M., von Essen, C., Weingrill, J., et al. 2015, A&A, 580, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mallonn, M., Bernt, I., Herrero, E., et al. 2016, MNRAS, 463, 604 [NASA ADS] [CrossRef] [Google Scholar]
 Mallonn, M., Herrero, E., Juvan, I. G., et al. 2018, A&A, 614, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mancini, L., Southworth, J., Raia, G., et al. 2017, MNRAS, 465, 843 [NASA ADS] [CrossRef] [Google Scholar]
 May, E. M., Zhao, M., Haidar, M., Rauscher, E., & Monnier, J. D. 2018, AJ, 156, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
 McCullough, P. R., Crouzet, N., Deming, D., & Madhusudhan, N. 2014, ApJ, 791, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Meunier, N., Desort, M., & Lagrange, A.M. 2010a, A&A, 512, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oshagh, M., Santos, N. C., Ehrenreich, D., et al. 2014, A&A, 568, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oshagh, M., Triaud, A. H. M. J., Burdanov, A., et al. 2018, A&A, 619, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Öztürk, O., & Erdem, A. 2019, MNRAS, 486, 2290 [CrossRef] [Google Scholar]
 Pinhas, A., & Madhusudhan, N. 2017, MNRAS, 471, 4355 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Sing,D. K., Gibson, N. P., et al. 2013, MNRAS, 432, 2917 [NASA ADS] [CrossRef] [Google Scholar]
 Rabus, M., Alonso, R., Belmonte, J. A., et al. 2009, A&A, 494, 391 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rackham, B., Espinoza, N., Apai, D., et al. 2017, ApJ, 834, 151 [Google Scholar]
 Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ, 853, 122 [Google Scholar]
 Rackham, B. V., Apai, D., & Giampapa, M. S. 2019, AJ, 157, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Radick, R. R., Mihalas, D., Lockwood, G. W., et al. 1983, PASP, 95, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916 [Google Scholar]
 Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2015, Nature, 529, 59 [Google Scholar]
 Solanki, S. K. 1993, Space Sci. Rev., 63, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, K. B. 2016, ApJ, 817, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Strassmeier, K. G., Granzer, T., Weber, M., et al. 2004, Astron. Nachr., 325, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Strassmeier, K. G., Granzer, T., Weber, M., et al. 2010, Adv. Astron., 2010, 970306 [Google Scholar]
 Tinetti, G., VidalMadjar, A., Liang, M.C., et al. 2007, Nature, 448, 169 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Exp. Astron., 46, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J., & Yurchenko, S. N. 2019, Nat. Astron., 3, 1086 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Mean values of transit depth variations due to spots, , on WASP52 simulated transits.
All Figures
Fig. 1 Example of a stellar photosphere including modulating and nonmodulating 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 T_{ph} = 5000 K and a blackbody 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 nonmodulating 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 activityinduced 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 nonmodulating filling factors of δ_{0} = 0 and 0.1 are shown. The gray lines illustrate how the bivaluate nature of the temperature contrast effect can be resolved if more than two bands are used. (c) Difference in activityinduced amplitude for photometric bands BV R with respect to the I band as a function of the nonmodulating filling factor. A modulating filling factor δ_{M} = 0.05 and a spot temperature contrast ΔT_{sp} = 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 latitudeprojected spot coverage for each 1degree 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 1sigma 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 statisticallyequivalent solutions for the fitted StarSim parameters: the rotational period, P_{rot}, the faculatospot area ratio, Q, the spotphotosphere temperature contrast, ΔT_{sp}, 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 WASP52 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 latitudeprojected spot coverage for each 15degree longitude bin. 

In the text 
Fig. 8 StarSim model fits to multiband photometric groundbased 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 WASP52 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 bynonocculted spots. The dotted blue line is the nonchromatic 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 limbdarkening effect (in red in top panel). TR1 and TR2 show two example transits corresponding to low and highactivity – 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 Δr_{SP}∕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 lightshaded band. Solid dots indicate the timestamps of WASP52 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 TR1 and TR2, respectively. 

In the text 
Fig. 12 Chromatic effects due to spot distributions with identical filling factors. MAP1 and MAP2 are two adhoc 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 WASP52 at epoch BJD = 2 457 644.95 with the parametric distribution error bands. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.