© The Authors 2025

1 Introduction

Since the first exoplanet discoveries (Wolszczan & Frail 1992; Mayor & Queloz 1995), thousands of exoplanets have been found, and the focus of the scientific community has moved from discovery to characterisation. Common methods such as transits (e.g. Deeg & Alonso 2018) and radial velocities (e.g. Wright 2018) to derive planets’ properties require the analysis of the stellar signal to correctly detect, and thereafter to characterise, exoplanets. Stellar activity is the main factor limiting our ability to characterise exoplanets or even to detect them (e.g. Cameron 2020).

The activity of a star is defined by the evolution of its magnetic field. An enhancement of the intensity or a change in the local magnetic-field configuration on the stellar surface could be the origin of different processes that alter the stellar signal in different spectral bands and/or could directly impact the properties of a nearby exoplanets throughout the so-called star-planet interaction (SPI) processes (Vidotto 2019). One of the manifestations of stellar activity is the appearance of photospheric inhomogeneities, primarily spots and faculae, which enhance the variability of the stellar signal being, respectively, colder and hotter than the stellar photosphere (e.g. Pagano 2013).

In the case of transit method, one of the issues is the contamination of the stellar activity in the transit-depth determination (e.g. Ballerini et al. 2012; Rackham et al. 2023). When a transiting planet covers a spot or a series of spots, the difference in flux between photosphere and spot artificially alters the transit depth and shape causing an incorrect estimation of the planet radius and/or transit timing (e.g. Oshagh et al. 2013). This also applies when a spot’s effect is omitted from the out-of-transit correction of in-transit studies (Czesla et al. 2009). Analogously, in transmission spectroscopy studies, the stellar contamination dramatically affects our ability to detect chemical species in the upper atmosphere of the transiting planets (e.g. Edwards et al. 2020; Thompson et al. 2024) with spots that could mimic planetary spectral features (Micela 2015; Cracchiolo et al. 2021a,b; Rackham et al. 2023).

The description of stellar activity is crucial to determining exoplanet parameters or even detecting them. Several authors have proposed different approaches to describe stellar activity. On the one hand, analytical and/or numerical models describe stellar variability in terms of spot and faculae filling-factor variations without the need to know the position of the inhomogeneities on the stellar disc (e.g. Kipping 2012; Lanza et al. 2003). They are extremely fast in describing the photometric variability, but, due to the lack of the geometrical information, they are not able to describe spectroscopic effects. On the other hand, several models do include stellar surface reconstruction, and they are able to give a deeper insight into the stellar variability at the cost of more time-consuming retrieval processes (e.g. Boisse et al. 2012; Dumusque et al. 2014; Herrero et al. 2016; Ikuta et al. 2020; Chakraborty et al. 2024). Among those, they distinguish different geometrical descriptions or the inclusion of particular physical effects; for example, the blueshift inhibition in active regions (Dumusque et al. 2014; Herrero et al. 2016) due to the strong local magnetic field (e.g. Dravins et al. 1981) and spot temporal evolution (Herrero et al. 2016; Ikuta et al. 2020).

Here, we propose PAStar as an alternative model to describe photospheric activity in the form of spots and faculae, which represents a compromise between physical description and speed of evaluation. Due to its flexibility, it could easily be improved. It is different from present models, particularly with regard to the description of the inhomogeneities and on the geometrical construction. As for the former, the model could simultaneously describe spots surrounded by faculae, isolated faculae, and isolated spots, without any relation between their radii. This flexibility allows the description of complex stellar surfaces, thus expanding our ability to describe and correct for the effect of the stellar activity on observed data. On the latter, the inhomogeneities and – more generally – the stellar surface are built on the projected stellar disc, therefore avoiding approximations on the projection of the surface-area elements and, consequently, improving the precision of the synthetic products. Moreover, the modularity of the model allows us to modify and, thus, study the effect of many aspects of the model, e.g. the limb darkening, in detail.

The model could be applied to a wide range of stellar observations, from photometric to spectroscopic time series. The occulted flux by a transiting planet is saved to allow further and specific studies of the transmission spectrum.

We describe the general properties and assumption of the model in Sect. 2. We present several synthetic products in Sect. 3 to show the capabilities of the model. In Sect. 4, we present a comparison of photometric light curves from our model with respect to the SOAP code (Dumusque et al. 2014) to check the consistency of the results. In Sect. 5, we show the validation of the model against photometric solar data, while in Sect. 6 we present a validation of the model in the case of a more active facula-dominated synthetic star generated by the SOAP code. Finally, in Sect. 7 we present the discussion and conclusions.

2 Model

In our model, the stellar atmosphere is a combination of three components: the quiet photosphere, spots, and faculae. Each component is characterised by considering a stellar atmosphere at a given temperature. Our framework is not bound to any specific collection of stellar photospheric models. In this work, we used Phoenix models, which are known for their ability to reproduce the atmosphere of low-mass stars (e.g. Maldonado 2012, PhD thesis). Nevertheless, we note that other families of models, or even the assumption of a black-body spectra, can be considered.

The number of the inhomogeneities is arbitrary, as is their distribution across the stellar surface; i.e. we do not impose any preferential location. The model takes into account the effects of star inclination, Doppler shifts due to stellar rotation, and limb darkening. Also in this case, the model is not linked to a specific limb-darkening model, but it requires the limb darkening coefficient(s) as a function of the wavelength of the input stellar spectrum as input. Two formalisms have been implemented: a linear one, and the four-coefficients formula from Claret (2000).

The rotational period of the star and its radius and temperature are imposed. The model does not describe the intrinsic evolution of spots and/or faculae, only the variation in the position due to the stellar rotation.

2.1 Geometry

The model builds the stellar photosphere in spherical coordinates, with unitary radius and equally spaced in colatitude and longitude. From the spherical grid (R=1, Θ, Φ), we compute the related Cartesian coordinate grid (X, Y, Z) as follows: $$X=cos(i_{\star }-\pi /2)\cdot sin(\mathrm{\Theta })\cdot cos(\mathrm{\Phi })-sin(i_{\star }-\pi /2)\cdot cos(\mathrm{\Theta }),$$(1) $$Y=sin(\mathrm{\Theta })\cdot sin(\mathrm{\Phi }),$$(2) $$Z=sin(i_{\star }-\pi /2)\cdot sin(\mathrm{\Theta })\cdot cos(\mathrm{\Phi })+cos(i_{\star }-\pi /2)\cdot cos(\mathrm{\Theta }),$$(3)

where i_⋆ is the stellar inclination with respect to the line of sight, Θ is the colatitude, and Φ is the longitude. Star inclination ranges from 0 to 90°, which correspond, respectively, to a pole-on and equatorial view.

The inhomogeneities are assumed to be circular with the faculae not necessarily correlated to the spot, which is the case for isolated faculae. Any distortion effect on the visible disc arises from their projection on the disc; therefore, the inhomogeneities become spherical caps. This allows us to better describe their effect when they approach the limb during the stellar rotation in time series. faculae are built as circles surrounding spots, which result in isolated faculae when the spot radius is zero.

With these assumptions, for a given spot-facula configuration the stellar photosphere configuration is built once for all the times of a series. A window rotating oppositely to the star selects the visible surface portion (i.e. X>0) whose longitudinal disc centre can be evaluated as follows: $${\varphi }_c(t)=-2\pi \cdot (t-t_0)/P_{\star },$$(4)

where ϕ_c(t) is the longitude of the centre of the visible disc in the stellar surface, t is the time, t₀ is a reference time, and P_⋆ is the star’s rotational period. Therefore, the visible portion of the surface is selected considering the star inclination.

The flux from each surface element is evaluated as follows in the case of linear limb-darkening: $$\begin{array}{rl}F(\lambda ,t,j,k)=& Mask(t,j,k)\wedge \frac{f(\lambda )\cdot d{A}_{j,k}}{2\pi \cdot (1/2-c_1(\lambda )/6)}.\\ & [1-c_1(\lambda )\cdot (1-\mu (j,k))],\end{array}$$(5)

or in the four coefficient formalism $$\begin{array}{rl}F(\lambda ,t,j,k)=& Mask(t,j,k)\wedge \frac{f(\lambda )\cdot d{A}_{j,k}}{2\pi \cdot (1/2-\sum _{i=1}^4{a}_i\cdot c_i(\lambda ))}.\\ & [1-\sum _{i=1}^4{c}_i(\lambda )\cdot (1-\mu (j,k{)}^{i/2})],\end{array}$$(6)

where j and k are, respectively, the indexes of colatitude and longitude, t is the time, Mask(t, j, k) is a matrix of integers to select the proper flux (f(λ)) accordingly based on the element area occupation. μ(j, k) is the cosine of the heliocentric angle computed as $\sqrt{1-Y^2(j,k)-Z^2(j,k)}$. Coefficients a_i are normalisation constants of the four coefficient limb-darkening formulae whose values are a₁ = 1/10, a₂ = 1/6, a₃ = 3/14, and a₄ = 1/4 and obtained by integrating the limb-darkening profile over the stellar disc to give the total disc-integrated flux (e.g. Cracchiolo et al. 2021a). Coefficients c_i(λ) are the limb-darkening coefficients, and dA_j,k is the projected surface-element area, which is computed as follows: $$\begin{array}{rl}d{A}_{j,k}& =1/2\cdot |dl1\cdot (dl2-dl3)-dl4\cdot (dl5-dl6)|\\ dl1& =Z_{j,k+1}-Z_{j,k}\\ dl2& =Y_{j+1,k+1}-Y_{j,k+1}\\ dl3& =Y_{j,k}-Y_{j+1,k}\\ dl4& =Y_{j,k+1}-Y_{j,k}\\ dl5& =Z_{j+1,k+1}-Z_{j,k+1}\\ dl6& =Z_{j,k}-Z_{j+1,k}\end{array}$$(7)

where Z_j,k and Y_j,k are Cartesian coordinates of the trapezoid vertices on the projected Cartesian grid.

The flux from the single surface element is then added to a specific spectral bin, depending on the Doppler shift due to the stellar rotation. Therefore, to evaluate the wavelength shift the projected disc is divided into strips whose width (ΔY) is determined by the binning of the input spectrum as in the following: $$\mathrm{\Delta }Y(\lambda )=\mathrm{\Delta }\lambda /\delta \lambda$$(8)

with δλ given by $$\delta \lambda =sin(i_{\star })\cdot 2\pi \cdot Y/(c\cdot P_{\star })\cdot \lambda ,$$(9)

where Δλ is the binning of the spectrum, c is the speed of light, and P_⋆ is the rotational period of the star. Finally, the flux from each surface element is added to the corresponding spectral bin. This model can be used to generate forward models or to retrieve stellar inhomogeneity properties from observations.

3 Synthetic observations

In the following, we present several possible synthetic models generated as described in Sect. 2: a photometric light curve (Sect. 3.1), a planetary transit (Sect. 3.2) reproducing HD 209458 b planet characteristics (del Burgo & Allende Prieto 2016), and a high-resolution spectrum (Sect. 3.3) with a resolving power of 115 000. For this exercise, we fixed the stellar characteristics, which are presented in Table 1, and we show the differences by considering an immaculate (quiet) photosphere, with respect two different activity scenarios. The first scenario considers that stellar activity is dominated by spots, whereas the second also takes into account the effect of the faculae that encircle the spots. In Fig. 1, a stellar configuration is presented, in the case of the most complex scenario. The initial observed stellar surface of the time series has been marked by two dashed black lines in the surface configuration map whose view centre is selected by a dashed red line. As explained in Sect. 2, this is a moving window that selects the observed surface portion at a specific time: in this case t=0. The geometrical grid consists of 500 points in latitude and 1000 points in longitude, and it results in a pixel size of 0.36°, which represents a good compromise between computational cost and accuracy of the solution in typical photometric and spectroscopic studies. In the high-resolution spectral evaluation, we selected the same configuration as for the initial time of the photometric light curve. In the transiting planet case, the first mid-transit time in the series corresponds to t=1d.

3.1 Photometry

Here, we present a synthetic photometric time series generated from the configuration shown in Fig. 1. For its evaluation, the black-body flux σT⁴ has been considered for both photosphere, spots and faculae, with σ being the Stefan Boltzmann constant and T the temperature, considering the temperatures defined in Table 1. In Fig. 2 (left), we present six days of the synthetic time series, corresponding to 1.5 periods of stellar rotations. The initial star configuration matches the one shown in Fig. 1 (right), and the star has been rotated clockwise. Together with the white light curves, in Fig. 2 (right) we also present the projected filling factor as a function of time (or rotation), defined as the total covered area by inhomogeneities on the projected visible disc divided by π, i.e. the total projected area, being R_⋆ = 1. In the case of the complex scenario, the filling factor is due to both spots and faculae.

The initial star configuration marks a condition of very high activity in which inhomogeneities cover more than half of the visible disc. As the rotation proceeds, large spots approach the limb, and the flux rises substantially due to the effect of the limb darkening, reaching its maximum after half of the rotation; there, only few small spots are visible on the disc with a projected filling factor of ~0.1. At later times, a conglomerate of inhomogeneities approaches the visible disc, drastically increasing the filling factor (up to 0.6). Alternatively, the flux falls towards its minimum at t=3.7 d. At t=4 d, the star completes a rotational period (as does the configuration), and, therefore, the time series repeat themselves.

The presence of the faculae affects the flux. Although faculae have only been built with an increment radius that is one over three of the radius of the spots, their higher temperature with respect to the spots and photosphere compensates for their smaller filling factor. This effect leads to a maximum increase in the total flux of ~10%, with respect of the spot case, when their projected filling factor is 0.2.

3.2 Photometric planetary transit

We present the synthetic transit of a planet reproducing the characteristic of HD 209458 b but facing our synthetic active star. We considered the black-body fluxes for the stellar and inhomogeneity fluxes, as was done for the photometric time series in Sect. 3.1. To evaluate the synthetic transit, PAStar requests the coordinates of the planetary orbit in the projected disc; i.e. Y and Z coordinates and the ratio of the planet to the stellar radius (R_p/R_⋆). All of this information is taken directly from the python package PyLightcurve v4.0.1 (Tsiaras et al. 2016), which allows us to evaluate the orbits (and many other properties) of all the planets in the Exoplanet Characterisation Catalogue within the ExoClock project (Kokori et al. 2022) directly; this project also contains HD 209458 b. Although our synthetic star radius (R=1.3 R_⊙) differs slightly with respect to the one of HD 209458 (R=1.2 R_⊙; see e.g. del Burgo & Allende Prieto 2016, and references therein), for the sake of simplicity we did not change our star parameters; this results in an overestimation of the planetary radius (i.e. ~10%) when the R_p/R_⋆ is fixed. Although this discrepancy would give incorrect results in a specific study, here we aim to show the general effect of spots and faculae crossing events; therefore, the exact value of the planetary radius is not relevant. In this configuration, during the transit the star rotates by ~10°. For this reason, we considered the evolution of the inhomogeneities caused by rotation in the evaluation.

In Fig. 3, we present a single transit of HD 209458 b together with two snapshots of the stellar configuration at two times during the transit. As in the previous sections, we compared the immaculate case with the two activity cases (i.e. spots and spots+faculae). During the transit, the planet crosses two distinct spots (and faculae), allowing us to disentangle their effects in the transit depth. The first spot is located close to the centre, while the second crosses the star limb.

When the planet crosses a spot, we observe a rise of the flux during the transit that is more pronounced when faculae are not considered. This is because spots have a lower flux with respect to the photosphere; therefore, their occultation rises the relative flux with respect to the average out-of-transit flux. Faculae alternatively have a higher flux with respect to the photosphere, and therefore their occultation leads to a fall of the flux during the transit. These effects in combination with high activity cases or high spot coverage along the transit chord could lead to an incorrect estimate of the transit depth and, consequently, of the planet radius. If the spot occulted is close to the limb, the high spot contrast (or low flux) leads to a distortion of the transit shape for the same effect described above, which artificially lowers the transit duration. These are quite typical observed effects (e.g. Bruno et al. 2018).

The presence of spots also affects the transit depth when they are not crossed. They lower the flux during the transit substantially with respect to the value when an immaculate photosphere is considered (see the green line in Fig. 3 between the crossed spots). This effect can be mitigated when an accurate out-of-transit correction is implemented.

Fig. 1 Configuration of inhomogeneities on stellar surface. Left: flux on latitude–longitude (θ − ϕ) domain colour-coded and ranging from 0 (in dark violet, corresponding to the presence spots) and 1 (in orange, corresponding to the photosphere), to 1.1 (in bright yellow, corresponding to the presence of faculae). A dashed red line marks the centre of the visible disc, while the two adjacent dashed black lines mark the boundary of the observed disc. Right: projection of left mask in the visible disc in which the colour of each surface element accounts for the (left) flux and the pixel area, but not for the limb darkening. White lines mark a low resolution spherical grid of 3.6°, which is ten times lower than the actual grid spacing.

Fig. 2 Photometric time series computed by the model. Left panel: photometric time series obtained with parameters listed in Table 1 and considering black-body emission, in the case considering both spots and faculae (orange) or only spots (green) and in the case of an immaculate photosphere (dashed blue). Right panel: projected filling factors as function of time for the two activity scenarios considered, and colour-coded as in the left panel.

Fig. 3 Photometric planetary transit computed by the model. Left: transit of HD 209458 b in case of (dashed blue) quiet photosphere, (green) spot-only activity, and (orange) spot+faculae scenario. Fluxes are normalised by their out-of-transit fluxes. Middle and right: snapshots of stellar configuration at two times marked by the blue arrows in the left panel. The transit cord is marked by a dashed black line, while the planet position is marked by a transparent and white circle, with the planet extent correctly scaled to the stellar radius. Stellar flux is coloured-coded as in Fig. 1.

3.3 High-resolution spectra

To calculate high-resolution spectra, the code makes use of predefined spectral energy distributions (SEDs), to be combined depending on the configuration of surface inhomogeneities. Here, we used Phoenix spectra from the Husser et al. (2013) database¹, with the stellar characteristics listed in Table 1.

This database covers the wavelength range from 500 Å to 55 000 Å with resolutions of 0.1 Å in the ultraviolet (UV) bands (500–3000 Å), R500 000 in the optical and near-infrared (NIR) bands (3000–25 000 Å), and R100 000 in the infrared band (25 000–55 000 Å). The parameter space covers from 2300 K to 8000 K for the temperature, 0.0–6.0 for the logarithm of gravity, −4.0 to 1.0 for the Fe-to-H ratio, and −0.2 to 1.2 for the α to H ratio.

The high resolution of the Phoenix spectra in the optical and infrared bands results in a very computationally expensive evaluation of synthetic products. For this reason, we lowered the resolution of the spectra to match the one of real case studies. We chose to reproduce data from the high-resolution spectrograph HARPS-N (Cosentino et al. 2012), which has a resolution of R115 000 in the optical band between 3800 and 6900 Å. Before performing this operation, spectra from the database were modified in order to describe observed spectra from Earth rather than a vacuum; therefore, wavelengths were shifted following Ciddor (1996) by applying $${\lambda }_{air}={\lambda }_{phoenix}/f,$$(10)

where λ_air is the wavelength as if the star has been observed from the ground, λ_phoenix is the vacuum wavelength of the synthetic Phoenix spectra and f is the correction factor (Ciddor 1996) given by $$\begin{array}{rl}f=& 1.0+0.05792105/(238.0185-({10}^4/{\lambda }_{phoenix}{)}^2)\\ & +0.00167917/(57.362-({10}^4/{\lambda }_{phoenix}{)}^2).\end{array}$$(11)

After the correction on the wavelengths, we performed the spectrum degradation by applying a Gaussian filter on the spectra. We used the procedure gaussian_filter in the public available python package scipy to perform this task. Since the procedure works with pixels in images, we had to convert the sigma, needed by the filter, from angstrom to pixel. We selected a very strong line, i.e. H_α absorption line (6562.8 Å), to calculate the full width at half maximum (FWHM) and the resulting sigma by the following: $$\begin{array}{rl}FWHM& =\sqrt{{\left({\lambda }_0/R_{harpsn}\right)}^2-{\left({\lambda }_0/R_{phoenix}\right)}^2}\\ \sigma & =FWHM/(2\sqrt{2}\ln 2)/dA,\end{array}$$(12)

where λ₀ is the reference wavelength of the H_α line, R_harpsn and R_phoenix are the resolving power of, respectively, the HARPSN instrument and the Phoenix spectra, and dA is the bin width of the Phoenix spectra. We note that this procedure makes sense only in the case of non-saturated absorption lines; therefore, in cases where an H_α absorption line is saturated, another choice should be made.

After the previous steps, the stellar spectrum was calculated in the case of the configuration shown in Fig. 1 and in the two scenarios presented in the previous section, which are the complex activity (spots and faculae) case and the only-spots case. In Fig. 4, we present a selected and narrow spectral window of the synthetic spectrum, ranging from 5461 Å to 5465 Å. This spectral window is populated by FeI lines (5461.550 Å, 5462.959 Å, 5463.277 Å, 5464.280 Å, Nave et al. 1994), by a NiI (5462.493 Å, Litzèn et al. 1993) line, and by a blend of CrI (5463.928 Å, Kiess 1953) and NiI (5463.920 Å, Litzèn et al. 1993).

The primary effect of the inhomogeneities is to shift the level of flux in the selected range, depending on their projected filling factor. Secondarily, as can be noticed from the top right panel, the line centroid has been blueshifted by the presence of a conglomerate of spots on the right portion of the visibile disc (see Fig. 1); these are coming towards the observer. When a time series is considered, this effect results in a periodic shift of the line centroid correlated to the star rotation. This is the very typical effect in which radial velocities of a star are modulated by the stellar activity, in this case mostly dominated by spots.

Other important effects include, from one side, the fact that the broadening due to stellar rotation blends all absorption lines that lie within the maximum Doppler shift (~0.3 Å in this range or 16.5 km/s). On the other hand, the line broadening coherently builds the shift of the line centroid with respect the non-rotational case (lower panels), resulting in a much more evident modulation of the centroid shift (or radial velocity modulation) in the case of high rotational speed, or lower rotational period, when a time series is being considered.

Fig. 4 High-resolution spectra computed by the model. Top rows: broadened spectra by star rotation in the case of immaculate photosphere (blue), the complex activity scenario (orange), and the spots-only scenario (green). Bottom row: same as top rows, but without considering the rotational broadening in the computation. For comparison purposes, absolute fluxes (left panels) are presented together with the normalised ones (right). The normalisation trend has been evaluated with a Gaussian filter of 1200 Å as sigma.

4 Comparison with literature: SOAP code

In order to check the consistency of the synthetic products of our forward model, we compared PAStar outputs with an alternative model present in literature: SOAP code (v2.0, Dumusque et al. 2014). To derive this comparison, a synthetic photometric light curve has been generated with SOAP code, which evaluates photometric and radial velocity variations induced by active regions, in the form of spots and faculae. SOAP models faculae as spatially uncorrelated bright features with respect to spots with a centre-to-limb temperature variation that matches that of the solar case (e.g. Meunier et al. 2010). For this test, our synthetic star reproduces a solar-type star whose photometric variability is dominated by four faculae and whose characteristics are reported in Table 2.

We modified our code to employ the same centre-to-limb temperature variation of the faculae temperature of the SOAP model, as well as for the flux contrast evaluation, which considers a black-body radiation at 5293.4115 Å for both spots and faculae. These changes expand the flexibility of the code to adapt to different model requirements depending on the scientific cases to be studied.

SOAP employs a quadratic formula for the limb-darkening description whereas, in PAStar we employ a more complex fourcoefficient formula (Claret 2000). To derive a coherent comparison, we matched the limb-darkening description by simplifying it in both models with regard to the linear case.

Therefore, we compared the two codes on three photometric light curves, which differ on the stellar inclination but share the same spatial configuration for the four faculae considered. The result of this comparison is presented in Fig. 5. We also tested the difference in the modelled light curves due to the grid resolution in one of the presented cases; i.e. for i_⋆ = 45°. In this case, we present the comparison using the same default number (N) of grid points that SOAP uses to build the visible stellar surface – i.e. N = 300² – and the number we employed as a default for colatitude points, which is N = 500². We remark that in our geometrical construction, when the star inclination differs from 90°, with all longitudes being visible close to the pole, the resulting number of grid points in the visible disc is much higher than in the SOAP code due to the projection of the stellar surface on the Cartesian grid.

Differences between the two models fall below the ~0.004% threshold with our model, which, systematically, gives higher fluxes. We notice a trend in the discrepancies with the amplitude of residuals, which correlate with the star inclination, resulting in the polar-on solution (i_star=0°) to have the greatest deviation. This effect could arise from the different geometrical construction and/or resolution chosen in the two models; however, due to the very low amplitude of the residuals, we assert a good match between the two models.

The difference between the two models increases at the lowest resolution. The SOAP code benefits greatly from the increase in the grid points giving smoother profiles, but this is not requested by PAStar, which gives nearly identical results for low and high resolution. This is primarily due to the geometrical construction of our code that relies on the Cartesian projected grid to integrate the stellar disc, and it results in higher precision with respect to the SOAP code at similar resolution.

Fig. 5 Top panel: comparison of SOAP forward model with respect ours, with, respectively, coloured lines marking the SOAP solutions and circles marking ours; parameters listed in Table 2 are used. The equatorial-on solution is presented in blue, an intermediate inclination of 45° is in orange, while the polar-on solution is presented in green. Fluxes have been normalised to the photospheric and unperturbed flux. An insert in the panel shows a comparison of the two models using different numbers of grid points to build the stellar surface. The insert follows the same notation of the main panel, but two profiles are added to mark the solution obtained using 5002 grid points in the case of the SOAP code (dashed line) and ours (x marker). Bottom panel: residuals of two models evaluated as (FSOAP − FPAStar)/FSOAP · 100, considering the lowest grid resolution.

5 Validation of the model on solar data

The aim of this section is to test the ability of the model to describe the stellar activity in real scientific cases in which we know (almost) nothing about the star. However, to validate the model, we aim to reproduce both the star variability and the actual surface inhomogeneities’ configuration (locations and sizes). Although in the stellar case there are different attempts to retrieve the spot-facula configuration (e.g. Donati & Brown 1997), knowing the spot-facula configuration in an unequivocal way is the chimera of stellar science due to the lack of sufficient spatial resolution in current instrumentation. The only star in which this exercise can be done unambiguously is the Sun.

5.1 The Sun as a star

The Sun has been observed with a multitude of instruments from both space (e.g. Ogawara et al. 1991; Domingo et al. 1995; Kosugi et al. 2007; Pesnell et al. 2012) and the ground (e.g. Scharmer & Lofdahl 1991; Cosentino et al. 2012), and in many spectral bands, from radio to X-rays. Among all the possible data available for the Sun observed as a star, we selected an optical photometric time series from the Variability of solar IRradiance and Gravity Oscillations (VIRGO) Sun photometer (SPM) (Fröhlich et al. 1995, 1997; Jiménez et al. 2002) on board the Solar and Heliospheric Observatory (SOHO) mission (Domingo et al. 1995). Despite the simplicity of the information retrieved from a photometric time series with respect to a high-resolution spectrum, it still represents a viable test to efficiently probe many aspects of the model, the geometric construction of spots-faculae, the flux evaluation, and the limb darkening. Other effects such as the Doppler broadening due to stellar rotation, can only be probed with a spectral synthesis. However, considering the characteristics of the Sun, such an aspect can be neglected due to its slow rotation, which gives vsini =~ 1.6 km/s (Pavlenko et al. 2012).

VIRGO-SPM measures the spectral irradiance in different spectral ranges at the temporal cadence of 60 seconds. We analysed data for the measurement of the spectral irradiance at 5000 Å with a bandwidth of 50 Å and at the highest level of calibration (flag L2 mission long; publicly available at the SOHO Science Archive²).

Although the Sun offers unique opportunities to test and validate stellar models, it is the worst case when it comes to activity, because the Sun is not an active star. Here, ‘active’ refers to the filling factor of surface inhomogeneities such as spots and faculae. In active stars, the spot’s filling factor could reach, or even surpass, the 50% of the visible disc (e.g. O’Neal et al. 2004), while in the Sun few are reached only during solar cycle maxima. Therefore, this section does not represent an exhaustive validation of the model, it is only its first step, and it is important to determine the minimum activity that our procedure can retrieve. An application of the model to active stars is presented in Sect. 6. To derive this validation, data from the 24th solar cycle maximum were selected, in the window that starts from 2013-12-30T04:13:08.920 (international atomic time; TAI), and it extends for 41.2 days thereafter. Data in the selected temporal window are presented in Fig. 6. Due to the high cadence of the VIRGO-SPM instrument (i.e. 60 s) raw data have been binned down to 12 h. In each bin, the mean of the raw data is taken as bin value, whereas the standard deviation is taken as bin error, corresponding to an average error of 70 ppm.

Fig. 6 Solar flux in VIRGO-SPM green channel as a function of time. Raw data are shown in blue, while binned data are shown by the dashed orange line, together with the error bars. Reference time (Tref) is 1996-01-23T00:00:04.46 TAI.

5.2 Assumptions

To proceed with the validation, we made the hypothesis that the Sun is dominated by dark spots in the selected temporal window; i.e. most of the variability is related by the variation in the location of sunspots with zero temperature. Therefore, we also neglect the effect of faculae. In the case of the VIRGO green channel, this choice is supported by measurements of the intensity ratio between umbra and the photosphere by Albregtsen et al. (1984), which show a weak variation of the ratio through the disc at 5790 Å and values below 0.1, decreasing for shorter wavelengths. Moreover, for the sake of simplicity, we neglect the inclination of the Sun with respect to the line of sight, and we set i_⋆ to 90°. Although these assumptions could be awkward for a well-known object, which the Sun is, the reason for the validation is to treat the Sun as a star in a typical and almost blind search for activity characteristics. Furthermore, the rotational period is imposed at 26.5 d, which was retrieved from data as a time between two periodic minima.

A limb-darkening four-coefficient formalism was set; the values of the coefficients were calculated with ExoTETHyS v2.0.10 python package (Morello et al. 2020). The relevant parameters used for the calculation are presented in Table 3 together with the resulting coefficients. For this task, we do not set the exact values of the Sun’s parameters, but the closest values to a point grid of the Phoenix_2012_13 (Claret et al. 2013) database to avoid interpolation on the results.

5.3 Retrieval framework

The model has been coupled to a retrieval framework in order to obtain parameters best-fit values and their errors. As for retrieval framework, we use Multinest v3.10 (Feroz et al. 2009) to derive Bayesian inference through the python interface offered by the PyMultiNest v2. 12 (Buchner et al. 2014) package and with the following likelihood: $$\ln p(y_n,t_n,\sigma ,v)=-1/2{\displaystyle \sum _n}[{(y_n-F(t_n,v))}^2/{\sigma }^2+\ln 2\pi {\sigma }^2],$$(13)

where y_n and t_n are, respectively, data fluxes and times, $\sigma =\sqrt{{\sigma }_n^2+{\sigma }_j^2}$ with σ_n as data errors and σ_j as the white noise jitter term. These were introduced to account for an unknown source of errors. Observed data are modelled as follows: $$\frac{S(t)-S_0}{S_0}[\mathrm{ppm}]={10}^6\cdot (f(t_n,v)\cdot f_{\text{scale }}-1),$$(14)

where S(t) is the solar flux at time t; S₀ is a flux value with respect to the data that are shifted and scaled by the L2 mission long correction³; ppm stands for part-per-million; f(t_n, v) is the integral of Eq. (6) over the disc, i.e. j and k indexes, and supposing a unitary total photospheric flux; and f_scale accounts for data corrections and was chosen considering that the maximum model value, i.e. when spots are not present, should match the maximum possible value in the observed data. For this reason, we set f_scale to 1 + 10^–6 · y_max, where y_max is the data maximum value of the observed tine series, and the factor 10^–6 arises from the data scaling which are expressed in parts-per-million. Considering all the assumptions we have made, the parameter vector (v) is made by the white noise jitter term (σ_j), latitudes (θ), longitudes (ϕ), and radii (R) of the m spot being considered. Prior values for the free parameters are presented in Table 4. Due to the defined assumptions on i_⋆, the prior for the latitude has been restricted to the (0, 90)° range. This choice allows the presence of spots only in the northern stellar hemisphere, and it avoids the degeneracy in the solution of two spots that are symmetric with respect to the stellar equator.

Our strategy consists of running the retrieval repeatedly by increasing the number of spots (m) and by comparing the Bayesian evidence between two subsequent iterations. We take the iteration that results in a difference in Bayesian evidence with respect to the previous step of at least five as the best-fit solution; this is considered as statistically strong evidence following Jeffreys & Lindsay (1963). Subsequent fits will not be considered statistically relevant. To ensure a good sampling of the parameter space, we selected 4000 live points in the MultiNest algorithm. The sampling efficiency was set to 0.5 to achieve a good compromise between model selection and parameter exploration. Moreover, a multi-modal search was set to ‘true’ to allow efficient sampling in the case of the presence of multi-modal distributions, while all other MultiNest parameters are left at their default values. To remove initial and low probability values from posteriors, a cleaning procedure was applied to the results by selecting the samples with a weight over the 0.75 quantile of the weight distribution.

5.4 Results

As explained in previous subsections, the retrieval was performed starting by considering only one spot and then repeated by increasing their number. Moreover, since we want to not only explain the photospheric variability, but also to validate PAStar, we compared the surface configuration and the light curve from the model with the actual Sun spot configuration at similar times. To do this, we selected the available images of the Sun – on the date of the selected VIRGO observation – from the Helioseismic and Magnetic Imager (HMI) (Schou et al. 2012) on board the Solar Dynamic Observatory (SDO) mission (Pesnell et al. 2012) publicly available at the SDO-NASA database⁴. Within this database, the HMI continuum images were downloaded, which scans the FeI absorption line at 6173 Å with a FHWM of 75 mÅ. Although the images are in a slightly different band with respect to the temporal series to be fitted, we do not expect relevant differences in spot size and position, and we only made use of them for comparison purposes.

The Bayesian log evidences (log Z) for the different models we tested are presented in Table 5. As a result, we selected the configuration with two spots to better describe the data since they represent the best evidence among the explored cases. The last tested scenario, i.e. the three-spot model, has a log Z lower than the two-spot model; therefore, we stopped increasing the complexity of the model here (i.e. the number of spots). In Fig. 7, we present the best-fit model considering the maximum a posteriori probability (MAP) values for parameters together with the model spot configuration and the Sun spot configuration at similar times of HMI-continuum images. The corresponding posterior distribution is presented in Fig. 8.

The solution found by the fit consists of two well-defined spots that alternate in the visible disc. Spot longitude is retrieved with a high level of precision even for the smaller spot. However, latitudes are less constrained, and the posterior indicates a band where spots can be found rather than a precise and well-located position, as for the longitudes. This is not surprising because, while the longitude defines how the spot configuration modulates the photometric variation due to the rotation, the latitude controls the amount of the variation imposed by the spot presence due to the limb-darkening effect. This effect should correlate with the size of the spot; therefore, we would expect a correlation between the latitude and radius parameters, and this is exactly what the posterior shows. Another parameter that could correlate with radius and latitudes would be the spot temperature; however, this is set to zero here.

We report a difference in the shape of the synthetic profile with respect to the observed one in the two main minima. This effect could arise from the simplified assumption of circular spots, which in this case cannot reproduce the highly irregular shape of the observed spots, and it could be particularly crucial for the Sun in which spots are small and irregular (at best, they are of the size shown during the selected observation) and less evident for active stars, in which their size increases. This discrepancy is overcome by high values of the white-noise jitter, which has a MAP comparable with the standard deviation of residuals (σ_res = ~240 ppm). Another possibility to explain this discrepancy is the neglect of the spot flux and the fact that the limb darkening could contribute to characterising the variation of the photometric profile. However, finding the temperature of sunspots is beyond the scope of this work, in which the Sun is treated as a star in the validation process.

5.5 Applicability of the model

In the previous subsection, the Sun is considered as a star from a methodology point of view. We derived the surface inhomogeneity reconstruction as we would apply the method to a far away star of which we do not know the spot temperature or star inclination. We note that in solar data, systematic errors dominate measurements due to a very high signal-to-noise ratio (S/N), evaluated as < y/σ >, where y is the stellar signal (flux) and σ is its error. Here, we tested the model with all the assumptions above, but in a much lower S/N scenario, which could mimic that of stellar observations. The question we addressed therefore concerns the ability of our method to retrieve the spot configuration when the S/N decreases to levels of the typical photometric observations from space (e.g. Transiting Exoplanet Survey Satellite; Ricker et al. 2015). To artificially decrease the S/N, we multiplied data errors by a factor of 30 and repeated the iterative procedure explained in Sect. 5.3. In Fig. 9, we present the result of this process. In this exercise, the solution considering one spot is the most statistically significant, only the bigger spot is retrieved with a broader constraint in the latitude parameter with respect the previous case. This can be seen from the posterior distribution presented in Fig. 10. We also tested other S/N scenarios; however, for factors over 30, the MAP parameters retrieved do not recover the actual position of the spot on the solar disc correctly. Therefore, this configuration is about the minimum activity level that PAStar can analyse due the combination of a small sunspot filling factor and the imposed S/N. For the sake of completeness, in Fig. 11 the retrieved solar configuration increasing the error by a factor of 100 is presented. The difference between the actual latitude of the spot with respect to the retrieved one increases with the increase of data errors, and the posterior distributions (not shown here) appear broader than the previous case, particularly in the longitude parameter.

Fig. 7 Comparison between best-fit (MAP) model of two-spot configuration versus Sun data. Left panels: synthetic photometric light curve is presented as blue line; binned solar data are shown via orange line with errors. A black cross marks the time at which the model (middle) and actual Sun (right) spot configurations have been shown. Times are shifted with respect to the first of the series. The spot position is reflected with respect to the equator for visualisation purposes, and stellar flux is coloured-coded as in Fig. 1.

Fig. 8 Posterior distribution of two-spot model. MAP solution is marked as a green line, while the median is marked as a red line. Dashed black lines in diagonal histograms mark, respectively, the 16th and 84th quantiles.

Fig. 9 Same as top panel of Fig. 7, but for the low-S/N case when one spot is considered.

Fig. 10 Same as Fig. 8, but for low-S/N case.

Fig. 11 Same as top panel of Fig. 9, but considering a lower S/N; i.e. 100 times the data errors.

6 Application to active stars: A facula-dominated scenario

In this section, we test the model and the retrieval framework against an activity scenario dominated by faculae. The synthetic stellar configuration was retrieved following the methodology described in Sect. 5, which considers an iterative procedure to determine the number of appropriate inhomogeneities on the stellar surface and, thereafter, the surface configuration. In this approach, the same retrieval process is derived for an increasing number of inhomogeneities, which halts when Bayesian evidence of the best-fit solution does not increase more than five with respect to the solution of the previous step. As a photometric light curve generated by SOAP, the equatorial-view case (i_⋆=90°) of the synthetic light curves presented in Sect. 4 is selected, which represents the case of an active star. This configuration employs four faculae distributed across the northern stellar hemisphere and whose characteristics are listed in Table 2. We derived this validation in the following two different scenarios:

1. We assume to know that the star is facula-dominated, and therefore we suppress the presence of spots on the surface (Scenario 1).

2. We do not have enough information on the star to formulate a hypothesis, and consequently we search for a general spot and facula configuration (Scenario 2).

In both scenarios, we retrieved the surface configuration by arbitrarily varying data errors (σ) between 10% and 40% of the photometric standard deviation in order to study limits and degenerations of the retrieval process in two S/N cases. Bayesian evidence of these retrievals are presented in Table 6. The best-fit configuration in the two scenarios are presented, respectively, in Figs. 12 and 13, while posterior distributions in the case of the best-fit solution for Scenario 1 are presented in Fig. 14.

As a result of the Bayesian evidence comparisons, the retrieval process favours the solution of Scenario 1, which matches the number of input faculae. Moreover, all the solutions of Scenario 1 are statistically more significant than their counter-parts in Scenario 2 as they give systematically higher Bayesian evidences. In the case of lower S/Ns, in both scenarios the best-fit solutions are characterised by only one inhomogeneity, with no matches with respect to the input solution. The obtained best-fit solutions in the two scenarios (Figs. 12 and 13) are similar in terms of surface configuration, but differ in the presence of very small spots in Scenario 2, which never exceed 0.04 R_⋆ and act as a source of noise by increasing the data jitter. Consequently, they result in a less constrained posterior distribution (not shown here). In general, the obtained surface configuration resembles the input one, but some discrepancy is present, as is evident from lower panels of Figs. 12 and 13. The smallest faculae (ΔR_faculae = 0.1 R_⋆) of the input configuration are missing from the obtained solution in favour of a non-existent inhomogeneity at a high latitude (θ ~ 60°) that has a smaller radius. The effect of the missing inhomogeneity was overcome by the model by placing the input equatorial faculae at a higher latitude to average the position of two close inhomogeneities. This behaviour is the same as the one shown in Sect. 5, where spots average the position of the actual spots on the Sun photosphere. Here, this effect applies to faculae that are separated by less than ~30° (centre-to-centre). When facula separation is higher, the correct position is retrieved, within a conservative 5 σ level, as obtained from the posterior distributions.

We extended this application to the case of i_⋆=60°, but only for Scenario 1 and for the highest S/N; i.e. considering data errors equal to 10% of the standard deviation of data. The result of the retrieval and the posterior distribution are presented, respectively, in Figs. 15 and 16. Differently from the previous case, although the four-facula best-fit solution has a greater Bayesian evidence (log Z = 806.5) than the three-facula case (log Z = 801.7), following our selection criterion we chose the latter as the most significant one. As a consequence, the smallest faculae are missing from the reconstructed surface with the others that are correctly retrieved within a conservative 10 σ level.

7 Discussion and conclusions

In this paper, we present a model to describe the stellar photosphere on particularly active stars. They are characterised by surface inhomogeneities, spots, and faculae, which define the observed variability. We modelled the stellar surface as a superposition of components, photosphere, spots, and faculae. They are spherical caps whose shape is distorted by their projection on the stellar surface, and they are described throughout independent stellar atmosphere.

The model describes the fundamental effects to be taken into account in order to describe the observed variability. These effects are the stellar inclination with respect to the line of sight, the limb darkening (independent for each component), and the broadening of spectral line shape due to the stellar rotation. It is differentiated from other models present in literature, primarily in the geometrical construction and the flexibility in describing different kinds of inhomogeneities on the same stellar surface, without any assumption on their radii.

The model produces several products as photometric and spectroscopic light curves, also in the presence of a transiting exoplanet, and for which the occulted flux by the planet is evaluated to allow further studies regarding the exoplanet atmosphere. This variety allows the model to be applied to several scientific cases, thus expanding its applicability.

To check the consistency of our synthetic products, a comparison with an analogous model present in literature was made. For this purpose, a synthetic facula-dominated Sun was generated with the SOAP code (Dumusque et al. 2014) for three stellar inclinations. We report a good match between solutions generated by the two codes with residuals within 0.004% with respect the SOAP solution, which can relate to the different geometrical construction. We show that our model gives smooth and reliable light curves at low resolution, while the SOAP code needs a much higher number of grid points.

The model has been tested against solar data in order to validate the geometrical construction and the evaluation of the photometric light curve. Although the Sun is not a particularly active star, it is the only object in which this validation could be made uniquely. For this reason, it offers a great challenge for the validation because of the dimension of the spots, which is generally small compared to other active stars.

By combining the model with a retrieval framework, we were able to constrain the latitude, longitude, and radius of the spots on the solar surface. We find that a model with at least two spots is able to explain most of the observed variability for the specific temporal window chosen, obtaining a good match between the modelled versus the actual spot location and size. Due to the simplicity of the spot geometry, modelled spots represent an ‘effective’ spot of the conglomerates observed. This effect leads to differences in the modelled light curve with respect to the observed one, and it is enhanced by the small filling factor of the observed spot. By testing the model in a real S/N scenario, we demonstrate the ability of the model to retrieve spot information when spots are as small as 0.05 R_spot/R_star. However, we stress that targets of applicability of the model are rather active stars in which spots’ sizes and their filling factors are much greater than in the solar case, and for those stars the effect of non-circular spots would be less evident.

The validation process was also extended to active facula-dominated stars by using synthetics photometric light curves generated by the SOAP code. The facula configuration has been retrieved in two scenarios in terms of modelled inhomogeneities, i.e. facula-dominated and spot+faculae scenarios, and in two S/N cases, i.e. 10% and 40% of the photometric variability, in the equatorial view case. As a result, the retrieved configuration matches the input one, provided that the S/N is sufficiently high. However, although the input configuration is globally retrieved in the case of the facula-dominated scenario, discrepancies in facula properties between the input and retrieved configuration could be obtained when the input faculae are closer than 30°. This behaviour is also observed in the solar case, but on a smaller scale, where the retrieved spot is an average of the observed spot configuration. We extended this analysis to the case of i_⋆ = 60° in the facula-dominated scenario and for the highest S/N. In this case, we were also able to retrieve the faculae out of the four of the input configuration, with only the smallest missing.

The model does not include differential rotation, which may be relevant depending on the stellar type being studied. In solar-type stars, equatorial rotation could differ significantly from the polar one. In the Sun, this effect results in a difference of ~10 days (e.g. Thompson et al. 2003) between rotational periods. As in the case of the Sun validation, when the effect of the differential rotation is neglected, but present, a different rotational period should be employed to correctly describe the photometric variability, which reflects the rotational period of the active latitude band in which inhomogeneities lie. In this validation process, it was evaluated from data as the difference between minima. However, our primary targets of interest are active stars in which the fast rotation limits the effect of differential rotation (Reinhold, Timo et al. 2013), which although present could be subsequently lowered by limiting the analysis to short temporal windows shorter than the equatorial rotational period.

Another physical assumption that could affect our results is the imposition of the spot and/or facula temperature to values that could not be correct and, consequently, leading to an imprecise determination of latitude and radius of inhomogeneities due to the strong degenerations expected. This degeneracy could be break by a simultaneous multi-band analysis that takes advantage of the different flux-contrast ratio of inhomogeneities, with respect to the photosphere, in different spectral bands (Ballerini et al. 2012).

In conclusion, the model we developed offers a flexible and valuable tool to describe the activity of stars when it is dominated by spots and faculae. It has a wide range of applicability due to the wide variety of products. The high flexibility in the inhomogeneity description allows our model to describe a wider range of activity cases with respect to the models present in literature. The description of stellar activity is a fundamental step in several astrophysical contexts and is covered by the method we presented. In terms of future developments, the flexibility of the model allows us to add effects to the physics description of the stellar surface (differential rotation of the star, spot and facula temperature...) with particularly ease, enhancing the ability of the model to describe more scientific cases.

Fig. 12 Comparison between best-fit (MAP) model of the four-facula configuration versus input SOAP configuration. Left panels: best-fit photometric light curve presented as blue line; SOAP synthetic light curve shown as orange line with its errors, obtained by adding in quadrature data errors and jitter. A black cross marks the time at which the model (middle) and input SOAP (right) facula configurations are shown. Stellar flux is colour-coded as in Fig. 1, with the exception of the facula flux that has been modified to 103 in order to enhance the visibility.

Fig. 13 Same as Fig. 12, but in the case of Scenario 2.

Fig. 14 Same as Fig. 8, but in the case of the facula-dominated Sun of the validation process presented in Fig. 12.

Fig. 15 Same as Fig. 12, but considering i⋆ = 60°. A white cross marks the pole of the star.

Fig. 16 Same as Fig. 14, but considering i⋆ = 60°. Although multiple modes were present in the posterior distributions, only the one containing the MAP is shown and used to evaluate uncertainties.

Acknowledgements

The authors thank the referee for the useful comments and suggestions that improved deeply our manuscript. The authors acknowledge the support of the ASI-INAF agreement 2021-5-HH.0. J.M. acknowledges support from the Italian Ministero dell’Università e della Ricerca and from the European Union – Next Generation EU through project PRIN MUR 2022PM4JLH “Know your little neighbours: characterizing low-mass stars and planets in the Solar neighbourhood”. A.P. acknowledge support from the INAF Minigrant of the RSN-2 no. 16 “SpAcES: Spotting the Activity of Exoplanet hosting Stars” according to the INAF Fundamental Astrophysics funding scheme. A.P. and G.M. acknowledge support from the European Union – Next Generation EU through the grant n. 2022J7ZFRA – Exo-planetary Cloudy Atmospheres and Stellar High energy (Exo-CASH) funded by MUR – PRIN 2022.

References

All Tables

All Figures

Fig. 1 Configuration of inhomogeneities on stellar surface. Left: flux on latitude–longitude (θ − ϕ) domain colour-coded and ranging from 0 (in dark violet, corresponding to the presence spots) and 1 (in orange, corresponding to the photosphere), to 1.1 (in bright yellow, corresponding to the presence of faculae). A dashed red line marks the centre of the visible disc, while the two adjacent dashed black lines mark the boundary of the observed disc. Right: projection of left mask in the visible disc in which the colour of each surface element accounts for the (left) flux and the pixel area, but not for the limb darkening. White lines mark a low resolution spherical grid of 3.6°, which is ten times lower than the actual grid spacing.

In the text

Fig. 2 Photometric time series computed by the model. Left panel: photometric time series obtained with parameters listed in Table 1 and considering black-body emission, in the case considering both spots and faculae (orange) or only spots (green) and in the case of an immaculate photosphere (dashed blue). Right panel: projected filling factors as function of time for the two activity scenarios considered, and colour-coded as in the left panel.

In the text

Fig. 3 Photometric planetary transit computed by the model. Left: transit of HD 209458 b in case of (dashed blue) quiet photosphere, (green) spot-only activity, and (orange) spot+faculae scenario. Fluxes are normalised by their out-of-transit fluxes. Middle and right: snapshots of stellar configuration at two times marked by the blue arrows in the left panel. The transit cord is marked by a dashed black line, while the planet position is marked by a transparent and white circle, with the planet extent correctly scaled to the stellar radius. Stellar flux is coloured-coded as in Fig. 1.

In the text

Fig. 4 High-resolution spectra computed by the model. Top rows: broadened spectra by star rotation in the case of immaculate photosphere (blue), the complex activity scenario (orange), and the spots-only scenario (green). Bottom row: same as top rows, but without considering the rotational broadening in the computation. For comparison purposes, absolute fluxes (left panels) are presented together with the normalised ones (right). The normalisation trend has been evaluated with a Gaussian filter of 1200 Å as sigma.

In the text

Fig. 5 Top panel: comparison of SOAP forward model with respect ours, with, respectively, coloured lines marking the SOAP solutions and circles marking ours; parameters listed in Table 2 are used. The equatorial-on solution is presented in blue, an intermediate inclination of 45° is in orange, while the polar-on solution is presented in green. Fluxes have been normalised to the photospheric and unperturbed flux. An insert in the panel shows a comparison of the two models using different numbers of grid points to build the stellar surface. The insert follows the same notation of the main panel, but two profiles are added to mark the solution obtained using 5002 grid points in the case of the SOAP code (dashed line) and ours (x marker). Bottom panel: residuals of two models evaluated as (FSOAP − FPAStar)/FSOAP · 100, considering the lowest grid resolution.

In the text

	Fig. 6 Solar flux in VIRGO-SPM green channel as a function of time. Raw data are shown in blue, while binned data are shown by the dashed orange line, together with the error bars. Reference time (Tref) is 1996-01-23T00:00:04.46 TAI.
In the text

Fig. 7 Comparison between best-fit (MAP) model of two-spot configuration versus Sun data. Left panels: synthetic photometric light curve is presented as blue line; binned solar data are shown via orange line with errors. A black cross marks the time at which the model (middle) and actual Sun (right) spot configurations have been shown. Times are shifted with respect to the first of the series. The spot position is reflected with respect to the equator for visualisation purposes, and stellar flux is coloured-coded as in Fig. 1.

In the text

	Fig. 8 Posterior distribution of two-spot model. MAP solution is marked as a green line, while the median is marked as a red line. Dashed black lines in diagonal histograms mark, respectively, the 16th and 84th quantiles.
In the text

	Fig. 9 Same as top panel of Fig. 7, but for the low-S/N case when one spot is considered.
In the text

	Fig. 10 Same as Fig. 8, but for low-S/N case.
In the text

	Fig. 11 Same as top panel of Fig. 9, but considering a lower S/N; i.e. 100 times the data errors.
In the text

Fig. 12 Comparison between best-fit (MAP) model of the four-facula configuration versus input SOAP configuration. Left panels: best-fit photometric light curve presented as blue line; SOAP synthetic light curve shown as orange line with its errors, obtained by adding in quadrature data errors and jitter. A black cross marks the time at which the model (middle) and input SOAP (right) facula configurations are shown. Stellar flux is colour-coded as in Fig. 1, with the exception of the facula flux that has been modified to 103 in order to enhance the visibility.

In the text

	Fig. 13 Same as Fig. 12, but in the case of Scenario 2.
In the text

	Fig. 14 Same as Fig. 8, but in the case of the facula-dominated Sun of the validation process presented in Fig. 12.
In the text

	Fig. 15 Same as Fig. 12, but considering i⋆ = 60°. A white cross marks the pole of the star.
In the text

	Fig. 16 Same as Fig. 14, but considering i⋆ = 60°. Although multiple modes were present in the posterior distributions, only the one containing the MAP is shown and used to evaluate uncertainties.
In the text