Forward modelling of Kepler-band variability due to faculae and spots

Variability observed in photometric lightcurves of late-type stars (on timescales longer than a day) is a dominant noise source in exoplanet surveys and results predominantly from surface manifestations of stellar magnetic activity, namely faculae and spots. The implementation of faculae in lightcurve models is an open problem, with scaling typically based on spectra equivalent to hot stellar atmospheres or assuming a solar-derived facular contrast. We modelled rotational (single period) lightcurves of active G2, K0, M0 and M2 stars, with Sun-like surface distributions and realistic limb-dependent contrasts for faculae and spots. The sensitivity of lightcurve variability to changes in model parameters such as stellar inclination, feature area coverage, spot temperature, facular region magnetic flux density and active band latitudes is explored. For our lightcurve modelling approach we used actress, a geometrically accurate model for stellar variability. actress generates 2-sphere maps representing stellar surfaces and populates them with user-prescribed spot and facular region distributions. From this, lightcurves can be calculated at any inclination. Quiet star limb darkening and limb-dependent facular contrasts were derived from MURaM 3D magnetoconvection simulations using ATLAS9. 1D stellar atmosphere models were used for the spot contrasts. We applied actress in Monte-Carlo simulations, calculating lightcurve variability amplitudes in the Kepler band. We found that, for a given spectral type and stellar inclination, spot temperature and spot area coverage have the largest effect on variability of all simulation parameters. For a spot coverage of 1%, the typical variability of a solar-type star is around 2 parts-per-thousand. The presence of faculae clearly affects the mean brightness and lightcurve shape, but has relatively little influence on the variability.


INTRODUCTION
Over its 9.5 year operation, the NASA Kepler satellite (Borucki et al. 2010) collected high precision photometry on ∼ 5 × 10 5 stars, resulting in the confirmed detection of > 2500 exoplanets. Activityinduced stellar variability is a dominant source of the intrinsic noise in both transit , Oshagh et al. 2014and Kirk et al. 2016) and radial velocity surveys (Saar & Donahue 1997, Desort et al. 2007, Lagrange et al. 2010, Meunier et al. 2010a, Haywood et al. 2014and Oshagh et al. 2017) that can hinder the accurate detection and characterisation of planets with active stellar hosts. As a result, Kepler has detected fewer low-mass planets than expected . In the process however, the high-precision photometry collected has provided information on stellar photometric variability (Basri et al. 2013), which will continue to be an important consideration when analysing data from Kepler successors such as TESS, CHEOPS, PLATO and ARIEL.
Lightcurves of late-type stars show a quasi-periodic behaviour on timescales relevant for planetary transits (hours to days). This variability is mainly caused by the presence of photospheric surface fea-E-mail: l.johnson17@imperial.ac.uk tures, manifestations of surface magnetic activity, which evolve in time and appear and disappear on a star's visible hemisphere as it rotates. The most prominent of these features are dark spots and bright faculae (Berdyugina 2004, Domingo et al. 2009, Solanki & Unruh 2013and Radick et al. 2018.
The link between the presence of surface features and stellar variability can be uncovered by observing the Sun as a star (Willson &Hudson 1991 andFröhlich 2013). By comparing high-resolution solar disc images and simultaneous solar photometry, spot and facular signatures can be identified in solar lightcurves. Large dips are observed in solar lightcurves at times corresponding to spot crossings, and double-peaked signatures are observed during the passage of faculae across the solar disc (Fligge et al. 1998, Fligge et al. 2000and Fröhlich 2002. The double-peaked shape results from solar faculae having large intensity contrasts near the limb of the disc, but being difficult to distinguish from the quiet photosphere close to disc centre in the visible broadband spectrum.
For more distant stars, photometric lightcurve analysis is often used to infer the presence and surface area coverage of surface features. However, there are limitations to this approach: photometry alone cannot provide reliable information on feature latitudes and is prone to underestimating coverage areas as it is insensitive to rotationally invariant features such as latitudinal bands and polar fea-tures (see Solanki & Unruh 2004, and references therein). This can be overcome by taking measurements in more than one band, allowing colour excesses to be detected and pointing to the presence of rotationally invariant features. There also exists a degeneracy between spot areas and contrasts Silva-Valio et al. 2010). This is because a larger lightcurve dip can result from either a spot with larger area or lower brightness/temperature (although the width can help constrain the feature size). Basri & Nguyen (2018) noted that almost all stellar lightcurves can be fitted using a '2-spot' model (see also Torres et al. 1972, van Leeuwen et al. 1987and Luger et al. 2021. The fact that we never see 'simple' feature distributions like this on the Sun raises the question of whether accurate spatial information on stellar surface feature distributions can be inferred through lightcurve analysis. In order to obtain a metric of a star's activity level from photometric observations, Basri et al. (2011) introduce the range variability, a simple, but useful quantity computed as R var = P(95) − P(5), where P(5) and P(95) are the 5 th and 95 th percentiles of the lightcurve amplitude respectively (calculated over a chosen time interval). In Basri et al. (2013), R var was calculated for the Sun during Solar Cycle 23 over 30 day intervals (just over one solar rotation period, P ≈ 24.5 days at the equator). R var correlates very well with other activity indices on the Sun (Salabert et al. 2017) and reflects the solar cycle. Stars with larger amounts of surface features are generally expected to have higher R var values. While R var is sensitive to instrumental noise when at comparable levels to that of the stellar activity or higher, it is much more robust than the untrimmed lightcurve amplitude (see Appendix A).
Several recent studies have focused on forward modelling lightcurves of active late-type stars, e.g. Rackham et al. (2018Rackham et al. ( , 2019 and Sarkar et al. (2020). These studies use synthetic spectra derived from 1D model atmospheres to approximate intensities for the quiet star and surface features at chosen effective temperatures. While this approach is sufficient to model spots, it neglects the strong dependence of facular intensities on their position on the stellar disc, treating them as 'hot spots'. At most wavelengths, solar faculae appear brighter when viewed at a greater angle near the disc limb because the magnetic flux tubes that manifest as faculae at the photosphere radiate mostly through their hot walls (Spruit 1976). In Herrero et al. (2016) and Cauley et al. (2017), the empirically-constructed facular limb brightening law of Meunier et al. (2010b) is used to reproduce this effect. This was extrapolated to a solar model (Borgniet et al. 2015) and later to stars of spectral type F6-K4 ) based on preliminary runs of the simulations used in this paper (Norris 2018).
Since there are no direct observations of faculae on stars other than the Sun, the only possible solution is to rely on simulations that are capable of taking the interaction between the matter and magnetic field into account, i.e. magnetohydrodynamics (MHD) simulations (e.g. Salhab et al. 2018). 3D magnetoconvection simulations were used to calculate both quiet star (Danilovic et al. 2010, Hirzberger et al. 2010and Beeck et al. 2013) and facular (Afram et al. 2011 andBeeck et al. 2015) intensity contrasts for a range of stellar spectral types. In Beeck et al. (2013Beeck et al. ( , 2015, the 3D radiation-MHD code MURaM (Vögler et al. 2005) was used to simulate 3D regions of a solar twin atmosphere for a range of surface field strengths. Norris et al. (2017) used the ATLAS9 spectral synthesis code (Castelli & Kurucz 2003) emergent intensities from these regions. Consequently, mean intensity contrasts were derived for the quiet photosphere and facular regions with a range of mean magnetic field strengths in the wavelength range 149.5 nm ≤ λ ≤ 160000.0 nm, finding good agreement with solar observations at dif-ferent heliocentric angles. Facular contrast values based on MURaM calculations were used in Schrijver (2020) to model photospheric facular appearances on Sun-like stars in the ultraviolet (UV), visible and infrared (IR).
We developed the stellar variability code actress to use quiet star intensities and facular contrasts derived from MURaM simulations  in forward modelling the lightcurves of active G2, K0, M0 and M2 stars. We first determine the lightcurve signatures resulting from individual surface features for different spectral types and viewing angles in the Kepler band. Next, we use a Monte-Carlo approach, populating actress model stars with Sunlike feature distributions and calculating rotational lightcurves and their corresponding R var to quantify the possible effect of varying our physical parameter inputs on variability levels.
The paper is structured as follows: in Sect. 2, our lightcurve modelling approach is detailed. The simulation inputs and setup are described in Sect. 3. Results of the investigations performed are presented in Sect. 4 and discussed in Sect. 5. The conclusions are given in Sect. 6.

MODELLING STELLAR LIGHTCURVES
actress calculates model lightcurves for rotating, magnetically active stars. The active stellar surface is modelled using a spherical HealPix map and populated with circular spots and facular regions with defined radii and positions (see Sect. 3). HealPix allows 2D orthographic hemispherical projections to be taken from any angle, representing the stellar disc as viewed by an observer. Limbdependent intensities are assigned to each pixel for the quiet photosphere and surface features (see Sect. 2.1). Lightcurves are calculated by rotating the star at a chosen inclination i (the angle between the line of sight and the rotation axis, i = 90 • corresponds to an equator-on view of the star) and accounting for the contributions of all pixels on the visible disc throughout a full rotational phase.
An example of the lightcurve modelling approach is shown in Fig. 1 for a Sun-like star. In the top panel, intensity maps of the stellar disc are shown, with stellar inclinations i = 90 • and 30 • . Dark spots are visible regardless of disc position (or viewing angle) but bright facular regions are more prominent towards the disc limbs. The resulting lightcurves shown in the middle panel reflect the surface distribution of features throughout a rotation. As the features are mostly positioned close to the equator in latitude, fewer surface features are visible at i = 30 • than at i = 90 • . The i = 30 • lightcurve is both smoother and has a smaller amplitude as a result (see Fig. 7 for equivalent lightcurves calculated with spots only). The HealPix map of the model photosphere in the bottom panel shows that the largest lightcurve dips correspond to the longitude regions most densely populated with spots. Numerical noise is present in the simulated lightcurves, resulting from the pixellation of features on both the surface map and disc projection. No additional noise is present.

Intensity calculations
To model active stellar discs, emergent intensities are needed as a function of disc position, or limb distance µ = cos θ where θ is the viewing angle. Limb-dependent intensities for the quiet photosphere and surface features are obtained by calculating local thermodynamic equilibrium (LTE) specific intensity spectra with ATLAS9 (Castelli &Kurucz 2003 andKurucz 2017) at 9 viewing angles, corresponding to limb distances 0.2 ≤ µ ≤ 1.0 Norris 2018). Figure 2 (inset) illustrates how µ corresponds to disc HealPix map representing the active stellar surface, cosine-scaled in latitude and flattened in longitude to resemble a solar synoptic map. The quiet photosphere is displayed in orange, facular regions are bright yellow and spot regions are dark blue. The crosses represent the centres of the stellar discs in the top panel. We note that the brightness of the features on the HealPix map does not correspond to the actual brightness (which is a function of observing angle).
position. The intensity spectra for the quiet photosphere and facular regions are calculated along sight lines through MURaM simulation cubes at different viewing angles, taking 3D geometry into account. Intensity spectra for the spots are calculated using 1D model atmospheres instead. This approach has been well-tested on the Sun and shown to reproduce total solar irradiance with very high accuracy .
In a given wavelength band (the Kepler band was used in this study), we fit a 3-parameter non-linear limb darkening relation, to the 9 angles from computations, following Sing et al. (2009) (see also Sing 2010). I(1) is the intensity at disc centre and a, b and c are limb darkening coefficients. This allows us to obtain the intensity at

MURaM-derived intensities for quiet stars and faculae
Solar faculae form when small magnetic flux tubes intersecting the stellar surface are swept by convective motion into intergranular lanes and the field is there intensified by the process of convective collapse (Parker 1978, Spruit 1979and Grossmann-Doerth et al. 1998; see Solanki 1993 for a review). Plasma is evacuated from within the flux tubes, decreasing the opacity. The decrease of the convective heating is overcompensated by the radiative heating from hot walls, resulting in regions brighter than the surrounding photosphere (Spruit 1976). At wavelengths dominated by continuum radiation in the visible, strong faculae appear brightest when viewed side on (near the disc limb) as the hot walls are most visible (Yeo et al. 2013). An approach that accurately accounts for 3D geometry is required to reproduce this behaviour (Keller et al. 2004 andCarlsson et al. 2004). Simulations are the only way to obtain reliable facular contrasts for stars of spectral types other than G2. A stellar disc without magnetic field is devoid of bright and dark surface features with lifetimes of days or longer 1 and will appear limb darkened. For the quiet photosphere, the MURaM code (Vögler et al. 2005 andBeeck et al. 2013) was used to run magnetoconvection simulations without magnetic fields (referred to as hydrodynamic, or field-free). Once relaxed, several simulation snapshots were taken at constant time intervals. To obtain quiet star intensities for all spectral types in our modelling approach, spectra were calculated along rays through the simulation cube snapshots at 9 viewing angles and averaged to obtain mean spectra at all angles, following Norris et al. (2017). We then convolve the spectra with the transmission curve of a given passband and integrate to obtain intensities (in photon units). Equation 1 was fitted with the Kepler-band intensities. The resulting quiet star intensity relations are shown in Fig. 2 Figure 2. Plot of MURaM-derived average intensity fits for the quiet photosphere (solid lines), weak-field facular regions (dashed lines) and strong-field facular regions (dotted lines) against limb distance in the Kepler band for G2, K0, M0 and M2 stars. Black crosses represent the 9 computed quiet star intensity values for each spectral type. Inset: actress-generated G2 stellar disc with lines of constant µ. The central point represents µ = 1.0, and the dashed lines (moving radially outwards from the centre) represent µ = 0.9, 0.8, 0.7, 0.5 and 0.2.
(solid lines), for all 4 spectral types considered in our models. We provide Kepler-band quiet star fit coefficients in Table. 1, labelled 'hydro' (see Appendix B for quadratic fit coefficients and equivalent TESS-band coefficients).
For the facular regions, vertical fields of mean strength B z = 100 G and 500 G were injected into the field-free MURaM simulations 2 . Once again, the simulations were allowed to relax before a series of snapshots were taken and mean spectra calculated. The facular region intensities derived from the 100 G and 500 G runs are used to model facular regions of low and high activity respectively, and are hereafter referred to as weak-field and strong-field facular regions. We note that B z is the mean over the simulated cube and does not correspond to the magnetic field strengths measured for faculae in surface magnetograms.
The MURaM-derived facular region intensities are shown in Fig. 2 alongside the quiet star intensities. We provide Kepler-band fit coefficients for the weak and strong-field facular regions in Table. 1, labelled 100 G and 500 G respectively. As the functional dependency between the quiet photosphere and facular region intensities (I phot , I fac ) is difficult to distinguish in Fig. 2 (especially for later spectral types), we show facular intensity contrasts, I fac /I phot − 1 in Fig. 3 (top and middle panels). The weak-field facular regions (top panel) have positive contrasts for all spectral types and µ, arising from facular brightening. They are brightest at the limbs of the stellar disc Table 2. Stellar parameters used for all available spectral types in our simulations. T phot are the equivalent effective temperatures of the quiet photosphere calculated from the MURaM-derived spectra. ∆T fac are the effective temperature differences from the photosphere calculated from MURaM 100 G and 500 G runs (for facular regions of weak and strong-field respectively). The spot temperature differences ∆T spot (based on Panja et al. 2020) and surface gravities log g are input parameters of the 1D radiative equilibrium model atmospheres used to calculate ATLAS9 spot spectra. (low µ) due to the limb-angle dependence of faculae . Strong-field facular regions (middle panel) have strong, positive contrasts at the limbs of the stellar disc for all spectral types. On Sun-like (G2 and K0) stars, facular region contrasts decrease to close to zero (cf. Yeo et al. 2013) at disc centre (µ = 1.0), whereas the M dwarf (M0 and M2) facular contrasts become negative, exhibiting spot-like behaviour near disc centre. This results from the presence of dark, pore-like structure in the simulated regions overwhelming the brightness contributions of the small-scale magnetic features forming the faculae. In general, peak contrast levels for the strong-field regions are much higher than for the weak-field regions (roughly 5 times higher for the G2 star). Both weak and strong-field facular contrasts tend to increase with stellar effective temperature (higher for earlier spectral types) at most limb distances. Solar facular contrasts (derived from 1D empirical models, Unruh et al. 1999) used in the Spectral And Total Irradiance REconstruction (SATIRE; see Fligge et al. 2000 andKrivova et al. 2003) model are also plotted for comparison with the G2, strong-field MURaM results. The contrasts are of roughly the same magnitude for both models, though these 1D contrasts used in SATIRE are lower near the limb and slightly higher near disc centre. Table 2 lists effective temperatures for the quiet star MURaM simulations T phot , as well as the effective temperature differences for the faculae ∆T fac (= T fac − T phot ) and spots (see Sect. 2.3). The effective temperatures for the MURaM simulations (quiet stars and faculae) represent averages over several snapshots. They were derived from disc-integrated intensity spectra where Equation 1 was used to extrapolate for µ < 0.2. Temperatures calculated for the several simulation snapshots in each run differ due to oscillations but agree well overall, with standard deviations of the order of 6 K for the G2 and K0 runs, and 2 K for the M0 and M2. T phot and ∆T fac are listed in Table 2. ∆T fac is negative for the strong-field facular regions on M2 stars.

Spot contrasts from 1D model atmospheres
Sunspots are formed when a high concentration of magnetic field suppresses convective motion at the photosphere over a large horizontal area. If this area is sufficiently large, the reduction in convective energy transport is greater than the radiative flow through the side walls, resulting in a locally cool, dark region. On the Sun, spots are made up of two components: darker umbra at the centre surrounded by less-dark penumbra. Solar umbra and penumbra have effective temperatures 1000 − 1900 K and 250 − 400 K cooler than the quiet photosphere respectively (Solanki 2003). Spot tem- Figure 3. Intensity contrasts between photospheric features and the quiet photosphere against limb distance in the Kepler band for G2, K0, M0 and M2 stars. Top and middle panels: MURaM-derived average facular region contrasts for runs with weak-field and strong-field facular regions, respectively. SATIRE model facular contrasts are also plotted (Unruh et al. 1999). Bottom panel: spot contrasts with temperature, T spot (see Table 2). peratures recorded for stars other than the Sun are subject to large uncertainties due to the limitations of (and discrepancies between) the indirect techniques required to observe them (see Appendix C for examples, andBerdyugina 2005 andStrassmeier 2009 for reviews). Spots have been simulated using the MURaM code for the Sun (Rempel et al. 2009a,b) and other main sequence stars (Panja et al. 2020), providing umbral and penumbral effective temperatures for G2, K0 and M0 stars alongside other fundamental parameters.
In our model, we do not distinguish between umbra and penumbra and instead assign an average effective temperature T spot across whole spot areas, as in Chapman (1987), Unruh et al. (1999) and Herrero et al. (2016). Solar irradiance variability modelling (Wenzler et al. 2006) indicates that this approach does not significantly affect the resulting variability. We use umbral and penumbral temperatures based on the results of Panja et al. (2020) for G2, K0 and M0 stars and extrapolate temperatures for the M2 by fitting quadratic relations to the existing data. A 4:1 penumbral-to-umbral ratio is assumed to calculate the average effective spot temperature. For the spectral types considered in our models, the effective temperature differences between spots and the photosphere ∆T spot (= T phot − T spot ) are listed in Table 2.
We follow Unruh et al. (2008) and use a 1D radiative equilibrium stellar atmosphere model with a cooler effective temperature than the quiet photosphere to represent spots. Mean MURaM quiet star spectra and 1D model spectra with equivalent effective temperature agree very well, although small differences are present for all spectral types. To minimise the effect of model differences on the resulting spot intensities, ATLAS9 spectra were calculated from 1D model stellar atmospheres at effective temperatures T phot and T spot for the quiet star and spots, respectively. The ratio of intensities was multiplied with the quiet star intensities from MURaM to provide spot intensities for our models. We show Kepler band spot intensity contrasts, I(T spot )/I(T phot ) − 1 in Fig. 3 (bottom panel) for all spectral types. In the Kepler band, spots always have negative contrasts several times larger than that of facular regions, regardless of limb distance. Spot contrasts become more negative with increasing stellar effective temperature (a direct consequence of the spot temperature differences chosen).

MODEL SETUP
We performed a number of investigations using the lightcurve modelling approach described in Sect. 2: obtaining typical variability levels assuming a Sun-like feature distribution and evaluating the impact of spot temperature, active region latitude and facular-to-spot coverage fraction on variability. Here, we describe how Sun-like feature distributions are generated and the assumptions made.

Latitude and longitude distributions
We aim to model brightness variability on rotational timescales. Focusing on the Sun, spots and facular regions are observed to occur in latitudinal bands in both the northern and southern hemispheres. The central latitudes of these 'active bands' are known to change throughout the solar activity cycle, beginning at higher latitudes at the start of each 11-year cycle and moving towards the equator as the cycle progresses (Hale et al. 1919). To model this, we use typical solar cycle averaged latitudinal distributions. These distributions are based on Borgniet et al. (2015), in which average spot and facular latitude distributions are presented for Solar Cycle 23 (May 1996 -October 2007), derived from USAF/NOAA sunspot group data and MDI/SOHO magnetograms (Scherrer et al. 1995) respectively. Transit mapping observations of the K4V star HAT-P-11 show a very similar spot distribution to the Sun (Morris et al. 2017). This result indicates that the solar cycle averaged distributions in our simulations can also apply to some K-type stars.
We fitted the solar spot and facular region latitude data with equatorially symmetric Gaussian distributions (for spots, φ = ±16 • , Fig. 4 (top panel). Feature longitudes are randomised. As the number of spots (and the activity level) increases, this latitudinal distribution results in an active band between approximately ±30 • . If spots (or facular regions) completely filled this band, surface area coverage would be 50% and disc coverage for a non-inclined star would be 61%. In the following we define the area coverages of spots A spot and facular regions A fac as percentages of the total stellar surface.

Size distributions
In Bogdan et al. (1988), it was found that the size distribution of sunspots obeys a lognormal profile (cf. Baumann & Solanki 2005). For the size distribution of spots in our simulations, we use a lognormal distribution fitted by Solanki & Unruh (2004) to sunspot umbral areas recorded at Mt Wilson between 1921 and 1982 in units of parts-per-million (ppm) of a solar hemisphere's surface area (10 −6 A 1 2 ). The size distribution was multiplied by a factor of 5 to represent total sunspot area (assuming a 4:1 penumbral-to-umbral area ratio) and the minimum spot size is set as 7.5 × 10 −6 A 1 2 (Solanki & Unruh 2004).
While individual faculae are much smaller than spots, solar faculae are grouped in active regions. Our models use average intensities from simulated 'facular active regions' which are not completely filled with individual faculae. We obtain our facular region size distribution by multiplying the spot size distribution by an arbitrary factor of 2.5. This rudimentarily mimicks how active regions observed on the solar surface emerge in groups, as well as increasing the computational efficiency of the simulation runs (by requiring fewer facular regions to reach the desired coverage level). The feature size distributions are shown in Fig. 4 (bottom panel).

Facular-to-spot ratio and spatial association
On the Sun, the area coverage of faculae, A fac,obs is much greater than that of spots, A spot . The observed facular-to-spot coverage fraction, Q obs = A fac,obs /A spot varies throughout the solar activity cycle: Q obs ≈ 10 at activity maximum, and Q obs > 40 at activity minimum (Chapman et al. 2001 andShapiro et al. 2014). The value of Q obs is heavily dependent on how solar faculae are identified.
As stated in Sect. 3.2, we model facular active regions that are not completely filled with individual faculae. Thus, the facular-tospot ratio Q in our models is not equivalent to Q obs . Magnetic pixels in solar images from the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics Observatory (Schou et al. 2012) are estimated to be 40% filled with faculae (Yeo et al. 2014, Chatzistergos 2017and Chatzistergos et al. 2019. We assume that MURaM boxes are equivalent to HMI magnetic pixels in facular filling, therefore using Q = 10/0.4 = 25 for solar activity maximum.
Following the spatial association of spots and faculae observed in solar images, we place half of the facular regions (by area) around spots. Every spot has a grouping of facular regions around it, with a coverage area roughly proportional to the spot area. The remaining facular regions are placed independently. A similar approach was implemented in Rackham et al. (2018Rackham et al. ( , 2019.

Monte-Carlo procedure
In our Monte-Carlo simulations, we sequentially added spots and facular regions to actress model stars, pseudo-randomly sampling feature latitudes, longitudes and sizes for 0 ≤ A spot ≤ 20%. Once added, the spots and facular regions do not evolve in time, remaining present throughout the time series of a simulated lightcurve (see Nèmec et al. 2020a,b for solar lightcurve modelling accounting for the temporal evolution of magnetic features). Adding more features allows us to probe higher activity levels. An initial facular-to-spot coverage fraction Q (which holds until 100% coverage is reached) was used to update the facular region coverage A fac after each spot is added and lightcurves were calculated at a series of regular intervals (finer gridding at low coverage levels 3 ) up to the maximum. The lightcurves were normalised by their means and the range variability R var was calculated for each. We typically generated N = 320 realisations for each parameter choice (all combinations of stellar spectral type and facular region field strength).
To allow the simulations to reach the desired coverage levels, the latitude constraint on facular regions is removed (enabling high latitude placement) once the total area coverage A tot = (A spot + A fac ) > 50% (as stated in Sect. 3.1, the low latitudes will be saturated with features at this coverage level). This is not entirely realistic as in practice, the whole active band latitude is expected to shift as activity level increases. High latitude facular region placement allows us to explore the fully-covered case and compare with other forward modelling studies (e.g. Rackham et al. 2018Rackham et al. , 2019. Once A tot = 100% is reached, Q will decrease from its initial value as spots continue to be added. In addition, we begin gradually increasing the size distribution of facular regions once A tot > 95%. This is because as the simulated stellar surface approaches complete coverage, the change in coverage with each new feature becomes increasingly small due to overlap with preexisting features. Increasing the feature sizes in this way ensures that total surface coverage is achieved within the alloted computation time. In this implementation, we do not account for superposition of features. The placement of a facular region will overwrite photosphere but not spots, and spot placement overwrites both photosphere and facular regions. See Appendix D for further details on how the feature distributions change from low to high activity throughout a typical simulation run.

Lightcurve signatures of single features
To illustrate the lightcurve variability contributions of spots and facular regions, we set up a toy model consisting of a single equatorial feature (radius r = 2 • ) on an otherwise quiet stellar surface. This single feature was modelled as a spot and a facular region of weak or strong field, for all available spectral types in the Kepler band. The resulting lightcurves (or feature signatures) are shown in Fig. 5 for i = 90 • , 60 • and 30 • and can be seen to reflect the intensity contrasts in Fig. 3. The weak-field facular region signatures (Fig. 5, top row) are observed to be bright for all spectral types. The signatures have flattopped profiles for the G2 and M0 stars, double-peaked profiles for the M2 star at i = 90 • and 60 • , and single-peaked profiles in all other cases. Decreasing the inclination generally reduces the signature sizes for the weak-field regions. This results from the equatorial feature being foreshortened towards the disc limbs (decreasing µ) and therefore having a smaller projected area on the disc, reducing its brightness contribution. This effect dominates over the increased brightness contribution due to the facular regions having higher contrasts towards the limbs.
For the strong-field facular regions (Fig. 5, middle row), doublepeaked signature profiles with central minima are observed for all spectral types when the star is viewed equator-on. For later spectral types, the peaks are less pronounced and spaced further apart and the central minimum has an increased depth. The Sun-like (G and K) facular signatures are brighter than the quiet photosphere everywhere except at disc centre. The M dwarf signatures exhibit brightening at the disc limbs, but are darker overall. As i decreases, the central minima become shallower for all types while the peak heights are preserved. This is due to the regions being positioned further from disc centre, where contrasts are darkest. Lightcurve signatures are also calculated using SATIRE model intensities for comparison with our results. The SATIRE signatures have less pronounced peaks and central minima than the G2 signatures.
For the spot regions, where we used the adopted relationship between the spot and photospheric temperature (Sect. 2.3), lightcurve dips are larger for hotter stars (Fig. 5, bottom row). As for the weakfield facular regions, signature sizes decrease as i decreases.

Lightcurve variability simulations
Monte-Carlo simulations were performed to investigate how lightcurve variability levels change with respect to stellar spectral type, feature area coverage, facular region field strength and inclination. We use a facular-to-spot area ratio of Q = 25 here, corresponding to Q obs = 10, representative of solar activity maximum (see Sect. 3.3).
In Fig. 6, R var against spot coverage is shown for all available spectral types and facular region field strengths at stellar inclination i = 90 • . The top, middle and bottom rows show the 'combined', 'spot only' and 'facular only' models respectively, for spot filling factors up to 5% of the total surface 4 Grey regions on the plot represent total area coverages A tot ≤ 50% in which facular regions are latitudinally constrained, and A tot ≤ 100% in which initial Q holds. Considering the 'combined' model (top row), mean range variability as characterised by R var increases steeply at first as A spot increases, before tapering off. Variability levels have a clear spectral type dependence, with higher variability exhibited for earlier types. This is a direct consequence of the spot and facular contrasts, relative effective temperature differences ∆T feature /T phot (see Table 2) and bolometric intensity contrasts become smaller for cooler stars (Panja et al. 2020). We observed a large spread in variability around the mean in all cases (more so at higher coverage levels) and show the interquartile ranges to represent this 5 . The spread is larger for earlier types (with higher variability). Increasing the facular field strength from weak (left column) to strong (right column) results in higher variability. The M2 star with strong-field facular regions exhibits a local maximum in variability around spot coverage 1%.
The inclusion of facular regions in our simulations has a small but noticeable effect on range variability in the Kepler band, shown here by comparing the 'combined' model results to the 'spot only' results ( Fig. 6, middle row). The inclusion of weak-field regions results in lower mean variability for all spectral types at coverages A spot 2%. At higher coverage levels, variability is slightly higher than the 'spot only' results for the Sun-like (G2 and K0) stars and roughly the same for the M-dwarfs (M0 and M2). For the strong-field case, the 'combined' model results are higher than the 'spot only' results for all spectral types at A spot 3%, and lower than the 'spot only' results at higher coverage levels. The local maximum in variability seen in the M2, strong-field case for the 'combined' model at lower coverage is absent in the 'spot only' results and is therefore a consequence of facular presence (M2 facular regions are dark, thus reinforcing rather than cancelling spot variability contributions).
The 'facular only' results shown in the bottom row of Fig. 6 reveal that for all spectral types and facular field strengths, facular-induced variability increases with feature area coverage up to A spot ≈ 1% (or A fac ≈ 25%), after which it begins to decrease again. This decrease results from the active latitude band (between approximately ±30 • ) being over half-filled, therefore additional facular region placement will make the band more uniform. The secondary peaks result from relaxing the latitudinal constraint on facular region placement, for A fac ≈ 50%. At A spot ≈ 3.8%, the star is completely covered by facular regions and variability therefore decreases to zero. The facular variability contribution is greater with strong-field facular regions, though still lower than when spots are included in the model in most cases (the G2 'facular only' variability with strong-field regions is comparable to the 'spot only' variability at A spot ≈ 1%).
While facular presence does not have a large impact on variability amplitudes (and therefore R var ) in the Kepler band, we found that inclusion of facular regions in our models dramatically influences the shape of lightcurves. Facular regions also affect the total brightness of the star (see Sect. 4.6). In Fig. 7, we show example lightcurves calculated with the 'combined' model with strong-field facular regions, alongside lightcurves with spots only. The 'combined' model lightcurves exhibit greater complexity, with more peaks and troughs than the 'spot only' lightcurves. This is more so the case for the lightcurves calculated at stellar inclination i = 90 • than at i = 30 • . Strong-field facular lightcurve signatures (see Fig. 5) have a highfrequency component at i = 90 • that is absent at i = 30 • . Increased lightcurve complexity resulting from facular presence has the potential to complicate the determination of stellar rotation periods from photometry (Shapiro et al. 2017.

Comparing simulations with observations
The spread in simulated variability around the mean is large at all coverage levels, as seen in Fig. 6. We examine this further by performing a rudimentary comparison between simulated variability 5 Interquartile ranges were plotted instead of 1-σ intervals because the variability distributions are asymmetric and non-Gaussian. distributions and observations 6 . Fig. 8 shows histograms of R var for the G2, weak-field facular region runs at i = 90 • , for a range of spot coverages. The observed 30-day Kepler-like variability over Solar Cycle 23 (the variability of a solar lightcurve created from composite VIRGO data, with a bolometric character comparable to Kepler data) and of the Kepler G-dwarf sample (Basri et al. 2013, Fig. 4) are plotted for comparison. On the right, the distribution of spot coverages for the simulation runs is shown. Each dashed line represents the centre of a spot coverage bin within which there are 320 simulation realisations. We note that the uniform distributions of spot coverages used have no physical basis (uniform coverage distributions are not seen throughout solar activity cycles, nor expected on other stars) and are used for illustrative purposes only.
In Fig. 8 (top panel), we show the simulated variability for spot coverages 0.05 − 0.25%, compared with the observed solar variability over Cycle 23. The high coverage cutoff in simulated spot coverage was chosen to reflect the observed maximal hemispheric area coverage from USAF/NOAA sunspot data for Cycle 23. The distributions appear slightly negatively skewed in log-space and therefore exhibit weak positive skewness in real space. This results from a higher probability of low-variability configurations in which features are near-uniformly distributed in longitude at all coverage levels. We also show the simulated variability for spot coverages 1.0 − 20.0% in Fig. 8 (bottom panel), alongside the observed variability of the Kepler G-dwarf sample. A low-coverage cutoff higher than for the Cycle 23 comparison was chosen because the Sun is thought to be relatively inactive compared to other solar-type stars (with known rotation periods close to solar, Reinhold et al. 2020). The highcoverage cutoff is also higher (the maximum spot coverage in our simulations).

Effect of T spot on lightcurve variability
The Kepler-band range variability is very sensitive to changes in T spot (at a fixed T phot ). To investigate the effect of different spot effective temperatures, we calculated mean R var of a G2 star for three different values. We chose 4041 K due to its use in Rackham et al. (2019) (representative of sunspot umbral temperatures) and 5100 K to represent average sunspot temperatures in agreement with SATIRE model fitting to SORCE/TIM data (Ball et al. 2011(Ball et al. , 2012. 4785 K is the G2 spot temperature based on an unweighted linear fit to all stellar observations in Appendix C. Increasing T spot (i.e. decreasing ∆T spot ) decreases the variability regardless of spectral type, facular field strength and stellar inclination angle i (although only i = 90 • results are shown here, we found this behaviour to be unchanged at different i). Changes in T spot result in an overall scaling of variability, with only small variations due to the variability contribution of the facular regions. This is illustrated by the thin grey lines in Fig. 9 which represent the T spot = 4785 K results scaled by factors of 1.43 and 0.69. These factors correspond to the Kepler-band flux differences from the ATLAS9 spectra at T phot and T spot . In the T spot = 5100 K model, scaling breaks down slightly at lower coverage levels (A spot < 8.5%) as facular contributions become more important when ∆T spot is smaller. The difference between the T spot = 4041 K results and the scaled variability (greater Figure 6. Range variability R var against A spot for the 'combined', 'spot only' and 'facular only' models (top to bottom) for all spectral types at i = 90 • with weak-field (left) and strong-field facular regions (right). The two 'spot only' panels are identical. Filled lines represent the means and shaded regions represent the interquartile ranges. Grey regions represent coverages A tot ≤ 50% in which facular regions are latitudinally constrained, and A tot ≤ 100% in which initial facular-to-spot ratio (Q = 25) holds. at higher coverage levels) is a consequence of the maximum of the Planck function moving sufficiently to the red at low temperatures such that the linearity assumed for the scaling is no longer given.

Effect of stellar inclination and active latitudes on lightcurve variability
In the initial simulation runs, lightcurves were simulated for stars viewed equator-on (i = 90 • ), with surface spots and facular regions distributed in active bands with central latitudes (φ = ±16 • ) based on the solar cycle average. Here, we investigate the effect of changing the stellar inclination angle and shifting the active bands on a G2 star further from the equator, with central latitudes (φ = ±30 • ) based on the maximum mean emergence latitude of the solar cycle. This typically occurs at the beginning of a cycle. The resulting plots are shown in Fig. 10.
Comparing the mean R var calculated for the 'solar average' case at i = 90 • , 60 • and 30 • reveals that variability decreases as inclination decreases. This results from the lightcurve signatures of spots and facular regions being smaller for lower i (see Sect. 4.1). This behaviour is observed for all spectral types and facular region field strengths. Changing the active band latitudes also has a noticeable effect on the variability: increasing the central latitudes from ±16 • to ±30 • results in variability levels roughly 15% lower on average at i = 90 • , 8% lower at i = 60 • , and 10% higher at i = 30 • . This is because variability is highest when the active latitude band in which spots occur coincides with disc centre at a given stellar inclination. While spots are dominant in the Kepler band, facular regions (which occur in broader latitude ranges than spots) also contribute. When the spot band coincides with disc centre, facular regions will be clustered nearer the 'top' and 'bottom' of the disc where contrasts are  high, resulting in more high-variability configurations than otherwise.

Effect of Q on lightcurve variability
Here, we investigate the effect of varying the initial facular-to-spot area coverage fraction Q on the variability. As stated in Sect. 3.3, the assumed relation between Q and the observed solar facular-to-spot coverage fraction Q obs is Q ≈ 2.5 × Q obs . The facular contribution to variability is largest for stars with relatively low activity levels and spot coverage fractions (Lockwood et al. 2007, Hall et al. 2009and Radick et al. 2018. We probe the low feature coverage range A spot ≤ 2.5% with a much finer grid for a G2 star, viewed equatoron at i = 90 • . The results for Q = 10, 25, 40, 50 and 60 are plotted in Fig. 11 for a G2 star with strong-field facular regions as both R var (top panel) and as mean brightness normalised to the quiet star (bottom panel).
Varying Q has a relatively small effect on mean variability levels in the Kepler band compared to other simulation parameters. For the small facular-to-spot ratios (Q = 10), mean variability is slightly lower than the Q = 0 (or 'spot only') results. The high facular-tospot ratio results (Q = 25, 40, 50 and 60) increases variability for A spot 1.5% and decreases variability at higher coverage levels. Typically, increasing Q results in larger deviations from the 'spot only' results. When the latitudinal constraint on facular regions is relaxed (illustrated by the change from solid to dashed lines), the variability slightly increases relative to what would otherwise be expected. This results from the previously unpopulated higher latitudes having inhomogeneities present. The local maximum in variability just before Q begins to decrease (illustrated by the change from dashed to dotted lines, easiest seen in the Q = 60 results at A spot ≈ 1.6%) is an artificial effect due to the increasingly large facular regions that are placed to ensure 100% coverage is reached, resulting in some longitudes being much brighter than others before the photosphere is completely covered. The effects of varying Q are similar for weak-field regions, though less pronounced.
As spots and facular regions are added to the simulated stellar surface, the mean brightness Φ of the modelled stars changes. In Fig. 11, we show the change in mean brightness ∆ Φ relative to the brightness of a featureless star for weak (top right) and strongfield facular regions (bottom right). In the presence of spots only (Q = 0), the mean brightness decreases linearly with A spot due to the increased presence of dark spots. The gradient (d Φ /dA spot ) becomes positive once both the facular-to-spot area ratio and facular field strength are sufficiently high for facular brightening to overcome spot darkening, resulting in an (initially linear) increase in mean brightness with spot coverage. A slight decrease in gradient is seen when the latitudinal constraint on facular regions is relaxed (dashed lines) (see Sect. 3.4); this results from foreshortening of facular regions near the disc limbs overcoming the limb-brightening. The gradients sharply decrease when Q begins to decrease as 100% coverage is reached (dotted lines) and dark spots overwrite bright facular regions.

TESS-band variability
The Transiting Exoplanet Survey Satellite (TESS) launched in 2018, aiming to discover sub-Neptunes via transit photometry (Ricker et al. 2015). The TESS detector bandpass spans roughly 600 − 1000 nm, slightly redder and narrower than the Kepler bandpass (roughly 400 − 900 nm). In Fig. 12, we show the TESS-band mean variability results calculated for spot coverage levels 0−20%, alongside the equivalent Kepler-band results for all spectral types with weak-field facular regions (as in Sect. 4.2). Compared to the Keplerband results, TESS-band variability is lower for all spectral types and coverage levels. This is a direct consequence of the spot and facular contrasts being smaller in redder wavelength intervals (Nèmec et al. 2020b).

DISCUSSION
Using actress, we calculated model Kepler-band lightcurves for active, late-type stars of spectral type G2, K0, M0 and M2, featuring Sun-like surface distributions of spots and facular regions. We assume solid body rotation and model single period lightcurves only. The profiles of actress lightcurves result from the superposition Figure 11. Plots of R var (top panel) and ∆ Φ (bottom panel) against A spot for a G2 star with facular-to-spot coverage fractions Q = 0, 10, 25, 40 and 60, at i = 90 • with strong-field facular regions in the Kepler band. ∆ Φ is the relative change in mean brightness with respect to a featureless star. Solid lines represent total coverage levels at which facular regions are latitudinally constrained and initial Q holds. Dashed and dotted lines represent total coverages ≥ 50% at which the latitudinal constraint is relaxed, and dotted lines represent 100% coverage at which Q decreases due to spots overwriting facular regions. The solid, dashed and dotted lines are equivalent to the dark grey, light grey and unshaded regions in Fig. 6 respectively. of the signatures of all features on the model stellar photosphere. The intensity contrasts between the quiet star and surface features directly influence the size and shape of individual feature lightcurve signatures. Quiet star intensities and facular region contrasts are derived from MURaM 3D simulations of stellar surface magnetoconvection. Spot contrasts are calculated using ATLAS9 synthetic spectra with effective temperatures based on data from MURaM spot simulations (Panja et al. 2020). We calculate the proxy for variability, R var as the outlier-trimmed amplitude of a lightcurve. As in Jackson & Jeffries (2012, 2013 and Rackham et al. (2018Rackham et al. ( , 2019, we found that spot coverage level, spot temperature and feature distribution are the main drivers of Kepler-band variability at a given stellar inclination. Variability levels are also affected by spectral type, facular region field strength and facular coverage level.
Considering a star with only one type of surface feature (faculae or spots), as feature area coverage increases we found that mean R var initially increases proportional to √ A feature before tapering off at higher coverage levels, eventually 'turning over' and decreasing to zero. This is due to the fact that a half-filled active band will become more uniform as features continue to be added, reducing the variability. In Rackham et al. (2018Rackham et al. ( , 2019, a comparable modelling approach is used: randomly distributing large spots with radius 2 • between latitudes ±30 • . The results are quantitatively similar to ours, finding an initial increase in variability with area coverage up to a turnover point, then decreasing again, reaching zero variability at A spot ≈ 33 % when the active band is full (and therefore rotation- ally invariant). Using the same simulation parameters, we found the active band to be fully populated at A spot ≈ 50% instead. The key difference in our approaches is in modelling the geometry of the system. actress stellar discs are created using a single projected hemisphere of a HealPix 2-sphere, whereas Rackham et al. (2018Rackham et al. ( , 2019 apply a double cosine weighting kernel to half of a 180 × 360 pixel rectangular grid. This accounts for linear limb darkening and feature foreshortening in a rudimentary sense, but does not account for the surface area of a stellar hemisphere.
For individual lightcurves, the change in R var when adding photospheric features is highly sensitive to the longitudinal placement. For example, a typical Kepler-band lightcurve will have its largest peak corresponding to the longitude region containing the fewest spots. Placing a spot in this region will reduce the size of the lightcurve peak, therefore reducing R var . Placing a spot in the longitude region containing the highest density of spots will increase the size of the corresponding lightcurve trough, increasing R var . Depending on their placement and contrast, facular regions can enhance or diminish the variability. Mostly-bright facular regions (weak-field on all spectral types, strong-field on G and K stars) will increase the size of a lightcurve peak, or diminish the size of a lightcurve trough. The 50% of facular regions placed near spots (see Sect. 3.3) will therefore generally lower the range variability. Mostly-dark strongfield features on M stars have a similar effect to spots, diminishing lightcurve peaks and accentuating troughs. This effect is seen in the M2 strong-field results in Fig. 6, where mostly-dark facular regions amplify the variability at low coverage levels. We also find that including facular regions in our models dramatically changes the shape of lightcurves. Increased lightcurve complexity resulting from facular presence therefore has the potential to complicate the determination of stellar rotation periods from photometry (Shapiro et al. 2017. While the sensitivity of individual lightcurve variability to the longitudinal distribution of features contributes significantly to the large spread in R var around the mean, the mean itself is mainly dependent on the location of the active latitudes (see Fig. 10). Using a Sunlike latitude distribution (active bands centred on φ = ±16 • ), R var is highest when the star is viewed equator-on (i = 90 • ). This effect was also seen by  who analysed synthetic activity time series for F6-K4 stars. Shifting the active band latitudes to φ = ±30 • causes R var to decrease when i = 90 • , and increase when i = 60 • . This is caused by the active latitudes being positioned closer to disc centre when the star is inclined. While high latitude features like this are not seen on the Sun, Doppler images of faster rotating stars (Strassmeier 2009) reveal high latitude (Wolter et al. 2005(Wolter et al. , 2008 and even polar (Barnes et al. 2001 andBarnes 2005) feature emergence. This can arise from the polewards deflection of the buoyant magnetic flux tubes (that form spots and facular regions at the photosphere) as they rise through the outer convective envelope of the star (Schüssler et al. 1996). Coupling our model with simulations of the rise of magnetic flux through the stellar convection zone (e.g. Işık et al. 2018) would allow us to implement realistic feature distributions for stars more active than the Sun, as well as providing spot and facular region filling factors.
Varying the facular-to-spot coverage fraction in the range 0 ≤ Q ≤ 60 was found to have a relatively small, but significant effect on mean R var , compared to the other parameters probed in our simulations (see Fig. 11, top panel). We found that, for Q = 25, 40, 50 and 60 and in the coverage range probed, the inclusion of facular regions in our models results in increased R var at coverages below A spot ≈ 1.5%, and lower variability at higher coverages. This can be attributed to large, bright regions increasing lightcurve amplitude up until the photosphere is over half-covered. For the Q = 10 runs, the presence of facular regions reduces variability levels at all coverage levels. This results from the spatial association between spots and facular regions in our simulations ensuring a concentration of bright regions around the dark spots, reducing lightcurve amplitudes. The increased coverage of unassociated facular regions at high Q will result in more high-variability configurations. However, as R var is calculated from mean-normalised lightcurves, facular brightening leads to a slight reduction in lightcurve amplitudes overall, lessening the aforementioned effect.
The high sensitivity of the variability proxy R var to spot temperature in our simulations (see Fig. 9) highlights the importance of using an accurate T spot determination for each spectral type. Underestimating spot temperatures (and therefore overestimating their variability contribution) might wrongly suggest that faculae can be neglected. The T spot to T phot relationship implemented in this paper is based on data from MURaM simulations (Panja et al. 2020). In Appendix C, we show results based on an unweighted linear fit to stellar observations, the majority of which were determined using TiO-band modelling and lightcurve analysis methods. Temperature determinations with these methods are highly uncertain and often inconsistent with each other. This inconsistency results from the different methods often reflecting different physical temperatures; TiOband intensity is weakly temperature dependent at T ≥ 3800 K (Neff et al. 1995 andO'Neal et al. 1996) and is therefore sensitive to spots on late-K and M dwarfs, but only umbral temperatures on hotter stars. Lightcurves are sensitive to umbra and penumbra regardless of stellar effective temperature. In addition, there exists a degeneracy between contrast and filling factor in lightcurves that cannot be broken when single bands (e.g. Kepler) are used. TiO temperature determinations for spots (on Sun-like stars especially) appear to be consistently lower than those obtained through lightcurve analysis. More accurate spectral determination of starspot temperatures is becoming possible by modelling spot occultations during planetary transits in spectroscopic observations (see Espinoza et al. 2019). Effective spot temperatures determined using MURaM simulations show good agreement with the solar case and the relationship between spot and photospheric temperature has a similar gradient to the Rackham et al. (2019) fit to the dwarf star data in Berdyugina (2005).
The 'facular only' case was investigated (see Fig. 6, bottom row), showing that the removal of the high-contrast spots results in a dramatic reduction in range variability at all coverage levels. On the Sun, high facular-to-spot coverage fractions Q obs > 40 are observed during low activity periods of the solar cycle. In some periods, the solar surface has been observed to be devoid of spots (while still populated with faculae). During activity maxima, Q obs decreases dramatically (Q obs ≈ 10) due to the increased presence of spots. We ascertained from varying Q that the range variability R var is a poor probe of facular presence. However, the mean brightness (and lightcurve shape) is significantly affected by faculae. Understanding how facular presence affects mean brightness is an important consideration when determining planet sizes via transit photometry. As would be expected, a linear change in mean brightness with spot and facular region coverage is seen when the facular-to-spot ratio is constant (see Fig. 11, right). For G and K stars with typically bright facular regions (regardless of facular field strength), the higher the facular-to-spot ratio, the greater the differential increase in brightness. The inverse is true for M stars with strong-field facular regions that appear dark relative to the photosphere at most limb distances.
Our simulation results for the G2 star reveal that mean brightness increases with coverage level for Q = 25 with strong-field facular regions, and Q = 40, 50 and 60 regardless of facular field strength. The Sun exhibits a strong activity-brightness correlation; at activity maximum, TSI (Total Solar Irradiance) is higher due to the facular brightening overcoming spot darkening. Comparisons between the chromospheric activity and photometric variability of solar analogs (Hall et al. 2009) have revealed that they typically exhibit at most a weak activity-brightness correlation. This could result from their higher average activity levels relative to the Sun, therefore having higher spot coverages and lower facular-to-spot ratios.
Comparing actress-simulated variability for a solar type star (viewed equator-on with weak-field facular regions, spot coverages uniformly distributed between 0.05 − 0.25% and Q = 25) and the 30-day variability measured over Cycle 23 (Basri et al. 2013) reveals that the range of variability and typical variability levels are of the same order (see Fig. 8, top). Throughout a typical solar cycle, there are more periods of relatively low sunspot area coverage than otherwise. This explains why the simulated distribution (with uniformly distributed coverages) has fewer low-variability occurrences than the observed distribution. The simulated distribution also has fewer high-variability occurrences: observed solar surface magnetic fields (and the emergent features) are rarely distributed uniformly in longitude, instead concentrating in 'clumps' and therefore amplifying the variability. This effect is explored through lightcurve modelling in Işık et al. (2020) and could be approximated by implementing 'active longitude' bands in the simulations.
We also compared G2 simulations (spot coverage range 1.0 − 20.0%, same simulation parameters as in the solar comparison) with the observed variability of the Kepler G-dwarf sample (Basri et al. 2013) (see Fig. 8, bottom). The observations have a much stronger positive skew than the simulated data, suggesting that there is a much higher proportion of relatively inactive stars in the sample than is reflected by the uniform spot coverage distribution used. In addition, the range of stellar inclinations in the G-dwarf sample will result in lower variability overall, as variability is maximised (in the Kepler band) when the star is viewed close to equator-on. The highvariability tail in the observed data cannot be explained by higher spot coverages alone; the results in Sect. 4.2 indicate that the increase in range variability with spot coverage lessens at higher coverage levels. Like the high-variability occurrences in the solar data, this could result from some longitudinal regions having high concentrations of magnetic activity.
A more complete comparison between simulated and observed data will be conducted in future: applying our understanding of the solar activity cycle (how feature coverage levels and emergence latitudes vary throughout), updating simulation parameters such as spot temperature and implementing active longitude bands. We also plan to investigate the wavelength dependence of activity-induced variability: solar facular contrasts are greatly enhanced in UV and are comparable to spot contrasts (see Yeo et al. 2013, and references therein), whereas contrasts of both the facular regions and spots are greatly diminished in the IR.
In this study, we consider stellar variability on rotational timescales only (days to months), where spots and faculae dominate. Another source of variability on these timescales is the evolution of magnetic features, which we do not consider at present (see Nèmec et al. 2020a,b for solar variability modelling with feature evolution). Variability on shorter timescales (minutes to hours) results mainly from granulation at the photosphere and pressure modes. Instrumental noise also contributes to observed variability on shorter timescales (Shapiro et al. 2017). Long-term variability (years to decades) is a consequence of the stellar activity cycle, throughout which typical levels of spot and facular surface coverage fluctuate (Solanki & Unruh 2013).

CONCLUSIONS
In this paper, we present actress, a geometrically accurate model for the stellar variability of late-type stars due to faculae and spots on rotational timescales, without spot evolution. The advantage of our modelling approach over previous works lies in the treatment of facular active regions (on the stellar photosphere). Intensity spectra calculated from MURaM 3D magnetoconvection simulations for G2, K0, M0 and M2 stars are used to derive limb-dependent intensities for the quiet photosphere and facular regions in a chosen spectral band. We provide the limb-dependent intensity coefficients used in this work for the Kepler band (TESS-band coefficients are also provided, see Appendix B).
We calculate single-phase Kepler-band lightcurves in Monte-Carlo simulations, investigating the sensitivity of the lightcurve variability proxy R var to changes in simulation parameters (spectral type, stellar inclination, feature distributions, spot temperature, facular field strength, spot area coverage and facular-to-spot coverage fraction). Our results show good qualitative agreement with prior studies. We attribute the quantitative differences to our accurate geometric setup and prescribed limb-dependent intensities. We found that, for a given spectral type and stellar inclination, spot temperature and spot area coverage have the largest effect on R var of all simulation parameters. This confirms the spot-dominated nature of variability in the Kepler band.
The inclusion of facular regions in our models (and varying facular simulation parameters) has a relatively small influence on R var . This does not indicate that faculae are unimportant, but rather that R var is not a good indicator of facular coverage. We found that faculae have a strong influence on mean brightness levels and on the lightcurve shape, potentially inhibiting planet size determination via transit photometry and stellar periodicity measurements respectively.
APPENDIX A: ROBUSTNESS OF R VAR While we do not include instrumental noise in our simulations, we evaluated the robustness of R var compared to the untrimmed lightcurve amplitude by using a Gaussian noise generator to add instrumental noise to a set of actress lightcurves. The fractional change in both variability proxies is plotted against noise level in Fig. A1.
The presence of noise typically increases both R var and the amplitude. The mean increase for R var at a noise level of 1 ppt (roughly comparable to the estimated precision for long-cadence observations of a Kepler magnitude K p = 15.5 target) is of the order of 2%, while the corresponding mean amplitude increase is roughly 25%. However, the range of changes in R var and amplitude are large. At the 95% confidence level, the range in R var change is [−4%, 15%] and the range in amplitude change is [−17%, 53%]. These results confirm the relative robustness of R var to noise. Figure A1. Plot of fractional change in variability proxy against noise level for both the untrimmed lightcurve amplitude (red) and R var (blue). Filled lines represent the means and shaded regions represent the 95% confidence intervals.

APPENDIX B: TESS-BAND INTENSITY COEFFICIENTS AND QUADRATIC EQUIVALENTS
We provide limb-dependent intensity coefficients for the quiet photosphere and facular regions in the TESS band (Table B1), for all available spectral types and magnetic field strengths. Variability results calculated in the TESS band are presented and compared with Kepler-band variability in Sect. 4.7. In addition, we provide equivalent Kepler and TESS-band fit coefficients for the quadratic limb darkening law (Tables B2 and B3) where I(1) is the intensity at disc centre and a and b are limb darkening coefficients. Fig. C1 shows a collection of spot temperature measurements from the literature. In principle, these can be used to re-determine the spot temperature dependence on spectral type. Spot temperatures collected here were determined using several different observational methods (lightcurve analysis, TiO-band modelling, Doppler imaging, line-depth ratios and transit modelling). These data were collected by Berdyugina (2005)    In Fig. C2, plots of mean R var calculated with spot temperatures from the linear fit to the combined data are shown alongside the variability results using Panja et al. (2020) spot contrasts, for all spectral types with weak-field facular regions. Variability levels are higher for all spectral types due to the larger spot temperature differences from the observationally derived fit, and the spectral type dependence is lost due to the greater discrepancy between the fit and the Panja et al. (2020) spot contrasts for later spectral types. This proves spot temperature to be a critical parameter for visible-band variability, further highlighting the importance of its accurate determination for stellar lightcurve modelling.  D1 shows how the latitude distributions of spots and facular regions change with spot coverage (or activity level) for the Monte-Carlo simulations described in Sect. 3.4. At low coverage levels (A spot 1.9% or A tot < 50%), the feature distributions are governed by the probability distributions shown in Fig. 4 (top panel), with spot and facular bands centred on ±16 • . The facular bands are broader than the spot bands. As coverage level increases beyond this point, the latitude constraint on faculae is removed, allowing high-latitude placement with a uniform distribution. The regions nearer the poles are populated before the mid-latitude ranges as there is less available surface area  Figure C2. Plot of mean range variability R var against spot area coverage A spot for all spectral types, with observationally-derived spot temperature contrasts (from the fit to all data in Fig. C1) at stellar inclination i = 90 • and weak-field facular regions. Variability results using Panja et al. (2020) spot contrasts are plotted for comparison (thin, dashed lines). Figure D1. Spot and facular region latitude distributions against spot coverage for a single Monte-Carlo simulation run. A feature density of 1.0 represents total coverage. at higher latitudes. At A spot 3.8%, the stellar surface is completely covered with spots and facular regions and the facular-to-spot coverage fraction decreases as more spots are added. We note that this appears to occur at a higher spot coverage in Fig. D1 due to the binning. The bin centres correspond to the sampling intervals for the lightcurve simulations.