Forward modelling of brightness variations in Sun-like stars II. Light curves and variability

Context. The amplitude and morphology of light curves of Sun-like stars change substantially with increasing rotation rate: brightness variations are amplified and become more regular. This has not been explained so far. Aims. We develop a modelling approach for calculating brightness variations of stars with various rotation rates and use it to explain the observed trends in stellar photometric variability. Methods. We combined numerical simulations of magnetic flux emergence and transport with a model for stellar brightness variability to calculate synthetic light curves of stars as observed by the Kepler telescope. We computed the distribution of the magnetic flux on the stellar surface for various rotation rates and degrees of active-region nesting (i


Introduction
Planet hunting missions such as the Convection, Rotation and planetary Transits (CoRoT Baglin et al. 2006;Bordé et al. 2003), Kepler (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014) allow studying stellar brightness variations caused by transits of magnetic features as stars rotate.Such brightness variations were discovered for the Sun almost half a century ago (Willson et al. 1981;Willson & Hudson 1981).Since then the models of solar brightness variations have matured and are now not only capable of accurately reproducing most of the available measurements (see Solanki et al. 2013;Ermolli et al. 2013, for reviews) but they also provide a starting point for explaining a plethora of stellar photometric data (see, e.g.Lagrange et al. 2010;Meunier & Lagrange 2013;Meunier et al. 2015;Borgniet et al. 2015;Nèmec et al. 2020b).
Recently Nèmec et al. (2020b) (hereafter N20b) have combined the Spectral And Total Irradiance Reconstruction model (SATIRE, Fligge et al. 2000;Krivova et al. 2003) together with a surface flux transport model (SFTM, Cameron et al. 2010) to compute the power spectra of solar brightness variations as they would be measured at different inclinations, i.e. the angle between solar rotation axis and direction to the observer.These calculations helped to remove a number of important observational biases when comparing solar variability to that of other stars (Nèmec et al. 2020a;Reinhold et al. 2020).Notably, by employing the N20b model and using the approach developed by Witzke et al. (2018Witzke et al. ( , 2020) ) to extend it to stars with non-solar metallicities, Reinhold et al. (2021) have found that rotation periods of a majority of the G-dwarfs with near-solar age remain undetected.These results provided an explanation for the discrepancy between the predictions of the number of Sun-like rotators in the Kepler field and the actual number of detected ones (see van Saders et al. 2019) .
In this work, we make an additional important extension of the solar paradigm for modelling variability of stars rotating faster than the Sun.Namely, we combine the solar variability model of N20b, which was extensively tested against solar irradiance measurements, with the modelling framework for computing the surface distribution of magnetic flux on stars with solar fundamental parameters but various rotation rates developed by Işık et al. (2018) (hereafter Paper I).The Flux Emergence And Transport (FEAT) model presented in Paper I involves physicsbased calculations of the emergence latitudes and tilt angles of bipolar magnetic regions (BMRs) for given stellar rotation rates, and the subsequent modelling of the evolution of the radial magnetic flux at the photosphere via an SFTM.The FEAT model is self-consistently able to reproduce the observations of polar spots that appear on stars with rotation periods below about 3 days (see, e.g., Jeffers et al. 2002;Marsden et al. 2004;Järvinen et al. 2006;Waite et al. 2015).Recently, the FEAT model was successfully applied to the young solar analogue EK Dra (S , enavcı et al. 2021), to explain the Doppler images that indicated near-polar spots and extended spot patterns towards low latitudes.In the present work, we extend the model of Paper I to calculate brightness variations of stars with a variety of rotation rates observed at various inclination.We also allow for different degrees and modes of nesting of magnetic features on their surfaces (i.e. the tendency of active regions to emerge in the vicinity of recently emerged regions).We compare our results to the observational trends found by McQuillan et al. (2014) and Santos et al. (2021) from Kepler data and propose a possible explanation of these trends.
In Sec. 2 we briefly describe the FEAT model (see Işık et al. 2018, for a more detailed description) and explain how this model is extended for calculating stellar brightness variability.In Sec. 3 we discuss the resulting light curves (LCs).These LCs are then used to calculate the amplitude of the variability, which we compare to observations in Sec. 4. We discuss our findings within the frameworks of various other recent studies in Sec. 5, before we present our conclusions in Sec. 6.

Model
Our model consists of two building blocks: calculations of the surface distribution of magnetic features and the subsequent calculations of their effect on stellar brightness.

Flux emergence and transport
To simulate the magnetic flux emergence on other stars, we utilise the FEAT model (Paper I).
In essence, the FEAT model extends the pattern of emergence and evolution of the magnetic fields observed on the surface of the Sun to stars rotating faster than the Sun (and, thus, more active than the Sun).This is done in five steps (see Paper I, Fig. B1): Step 1.We adopt the synthetic record of solar active-region emergences during cycle 22 from Jiang et al. (2011).While this approach cannot be used to reconstruct the solar irradiance on a specific day, N20b have shown that Jiang et al. (2011) records allow reproducing the overall pattern of solar variability.Cycle 22 was chosen, as it represents a cycle of moderate to strong activity level, making it suitable for modelling the most active stars.
Step 2. We define the time-dependent emergence rate of BMRs on a star as S (t) = S (t) • s, where S and S are stellar and solar emergence rates, respectively, and s is a scaling factor.To reflect the observed rotation-activity relation, we followed Paper I and took s = ω ≡ Ω /Ω , where Ω is the rotation rate of a star and Ω is the solar rotation rate.
Step 3. The resulting input record of emergences is mapped down to the base of the convection zone using thin flux tube simulations (Schüssler et al. 1996;Işık et al. 2018) for the solar rotation rate Ω .
Step 4. The record of emergences at the base of the convective zone is mapped back to the surface, but in contrast to Step 3 using thin flux tube simulations for a star with a given rotation rate Ω .These simulations follow the rise of the flux tubes throughout the convection zone up to the surface of the star, where they emerge in the form of a loop with two footpoints of opposite polarity (i.e.BMRs).An important feature of these simulations is that they account for the Coriolis effect that pushes rising flux tubes towards higher latitudes for higher rotation rates Schüssler & Solanki (1992).
Step 5. Assuming that the size distribution of BMRs does not change with the rotation rate, we modify the emergence locations of flux loops, to simulate different degrees and modes of active-region nesting (see a similar approach by Işık et al. 2020), motivated by the observed activity complexes on the Sun.The details of this step will be described later in this Section.
Step 6.In this last step, the calculated locations of emergence and the tilt angles of BMRs are fed into the surface flux transport model (SFTM).The SFTM describes the passive transport of magnetic fields of the BMRs on the surface of stars, by taking into account both large-scale flows (meridional flow and differential rotation) and the diffusion of the fields (Cameron et al. 2010;Işık et al. 2018).
In the present work, we employ the steps outlined to obtain the surface distribution for stars with 1, 2, 4, and 8Ω .The emergence patterns of the BMRs (e.g. the input record) for 1Ω and 8Ω are shown in Fig. 1.The main purpose of this study is to model stellar variability on the rotational timescale.We follow the approach of Paper I and take the solar temporal profile of activity (even though it is unlikely to be observed on the faster rotating and more active stars, we discuss this choice later in this paper), but we model the light curves only over the four-year window centred at the activity maximum (marked in Fig. 1).Our light curves (and the resulting variability amplitudes) correspond to a representative stellar activity level scaled from the solar maximum.Hence they are largely unaffected by the temporal profile of the underlying activity cycle.
In Fig. 2 we present snapshots of the magnetic field distribution on the surface of stars with 1Ω (top row) and 8Ω (bottom row) for various inclinations.The figure nicely demonstrates the difference in the emergence frequency, as the 8Ω star is covered by more BMRs.Also the difference in the latitudinal distribution is remarkable.This is a result of Step 4 of the FEAT approach, which in turn leads to the formation of strong polar fields for the fast rotator depicted here.We also show snapshots of the magnetic field distribution of the surface of stars with Ω and 8Ω for various inclinations in Fig. 2. Additionally, the emergence pattern of BMRs is altered by introducing active-region nesting.This effect is observed on the Sun, albeit to a small degree (see, e.g., Pojoga & Cudnik 2002, who estimated that 40-50% of solar active regions can be associated to nests) and it was recently suggested that nesting can be substantially stronger on highly variable stars with near solar rotation rates (Işık et al. 2020).Following Işık et al. (2020), we introduce two modes of nesting (Step 5): the active-longitude (AL) and the free nesting (FN) modes.In both modes, we define a degree of nesting, p, which gives the probability of a given BMR to be part of a nest.In the AL mode, those nests are centered around two fixed longitudes with a separation of 180 degrees.If the BMR is drawn to belong to a nest, then it is shifted to one of the two ALs (with equal probability).Its new longitude is defined by a 1D normal distribution around the AL (with the standard deviation 10 • ), whereas the latitude is kept unchanged.This ensures that the latitudes of emergence still follow the general trends of the solar butterfly diagram.In the FN mode, nesting occurs around central latitudes and longitudes, which are randomly picked from the unaltered (i.e.non-nested) emergence record.If the BMR is drawn to belong to a nest then its emergence location is moved to a random location drawn from a 2D normal distribution with the mean at the nest centre having standard deviations of 2 • in latitude and 3 • in longitude.For more details on the nesting definition and choice of the parameters we refer to Işık et al. 2018 for the FN employed here and Işık et al. 2020 for AL.
We show a comparison between the two modes of nesting for 1Ω in Fig. 3.In order to distinguish more easily between the two modes, we show both time-latitude (left panels) and timelongitude diagrams, with the colourbar indicating the longitude and latitude each.The AL nesting is easier to identify in the timelongitude diagrams, whereas the FN nesting is easier to observe in a time-latitude graph.

Defining the area coverages
The output of the FEAT model consists of full stellar surface maps (360 • by 180 • ) of magnetic fields, with a resolution of 1 • × 1 • per pixel.As our brightness calculations rely on the area coverages of spots and faculae, we have to first convert the magnetic field maps into the surface distributions of spots and faculae.N20b did this by following the evolution of sunspots after their emergence to calculate the coverage of the solar disk by spots.While this approach proved itself to be accurate for calculating solar variability (see also Dasi-Espuig et al. 2014) it did not account for spots which appeared due to the superposition of the magnetic flux from different active regions.In other words, the main assumption of this approach was that all spots on the stellar surface have emerged as spots within the corresponding active region.While this is a good assumption for the present Sun it prohibits the formation of polar spots, which are observed on young and rapidly rotating G stars (Jeffers et al. 2002;Marsden et al. 2004;Järvinen et al. 2008;Waite et al. 2015;S , enavcı et al. 2021).These polar spots most likely form via flux superposition.We therefore take a slightly different approach than in N20b.We acknowledge that, apart from the intrinsic mechanism in the FEAT model it is theoretically also possible to form polar spots by overly increasing the rate of BMR emergence (Schrijver & Title 2001), or meridional flow rate (Holzwarth et al. 2006), but leaving the emergence latitudes solar-like.
To calculate the spot area coverage of a given pixel of the synthetic magnetogram returned by FEAT we define two thresholds: a lower cut-off, B min , and an upper saturation level, B max .The spot coverage of a given pixel is related to the field in the pixel following where α m,n s is the spot filling factor of a pixel with coordinates m and n of the magnetogram, and |B mn | is the absolute value of the field in this pixel.
In order to calculate the faculae area coverage in each pixel, we follow an approach similar to N20b by setting a saturation threshold B sat (see also papers describing SATIRE, e.g., Fligge et al. 2000;Krivova et al. 2003).If we find that a pixel is covered already partially by spots, we disregard this pixel for the facular masking.If a given pixel is spot-free, we calculate the faculae area coverage following: The values of the parameters are selected such that the current model returns the same rotational variability (see Sect. 2.4) for the solar case (i.e.Ω = Ω ) as the N20b model.

Calculating the brightness variations
Using Eq. 1 and 2 we obtain maps of area coverages per pixel of the visible stellar disc.What part of the full surface (hence the disc) of a star is visible depends on the inclination between the stellar rotation axis and the line-of-sight to the observer.The spectral irradiance, S (t, λ) (i.e., the spectral stellar flux, normalized to 1 AU), where t is the time and λ the wavelength, is then calculated by summing up the intensities weighted by the corresponding area coverages by magnetic features of a pixel as given by S (t, λ) = S q (λ) + mn k I k mn (λ) − I q mn (λ) α k mn (t) ∆Ω mn . (3) The summation is done over the pixels of the maps of the visible 2D stellar disc and the m and n indices are the pixel coordinates (longitude and latitude, respectively), α k mn is the pixel (m,n) coverage by magnetic feature k (in the present work faculae, umbra or penumbra), ∆Ω mn is the solid angle of the area on the stellar disc corresponding to one pixel, as seen from the distance of 1 AU.I k mn is the intensity spectrum of magnetic feature k observed at the location corresponding to pixel (m,n).We use the values computed by Unruh et al. (1999) with the radiative transfer code ATLAS9 (Kurucz 1992;Castelli & Kurucz 1994).
S q is the quiet-star irradiance, defined as S q (λ w ) = mn I q mn (λ w )∆Ω mn . (4) Note that the solid angles of the pixels, as well as the corresponding intensity values depend on the vantage point.Hence S (t, λ) is sensitive to the stellar inclination.The calculations presented in this work are performed in the Kepler passband, following where λ 1 and λ 2 are the blue and red threshold wavelengths of the filter passband, R(λ) is the response function of the filter and S (λ, t) is the spectral irradiance at a given wavelength and time t, h is the Planck constant, and c is the speed of light.

Defining the parameters of the model
As mentioned previously, we fix the parameters of the model introduced in Sect.2.2 to return the same level of rotational variability, represented via R var , of the solar case as N20b.N20b showed that their model is able to reproduce the observed brightness variations of the Sun.The N20b model and the model we develop here are both based on the SFTM, with similar underlying statistical BMR emergence records.We therefore use the N20b model as our reference in the present work.For this, we considered the four-year interval (indicated by vertical dashed black lines in Fig. 1) around the maximum of the synthetic cycle.
We then split the time series into 90-day segments, detrended the new time series by their mean value, and calculated the difference between the extrema in each of the segments using the approach outlined above and that of N20b.We note that we directly consider the difference between the extrema instead of the differences between the 95th and 5th percentiles of sorted flux values, as is usually done in the literature with the more noisy Kepler measurements (see, e.g., Basri et al. 2013).
We show a comparison between the R var values for the solar case as returned by N20b and the present model in Fig. 4. The best-values for the parameters B min , B max , and B sat were found to be 60, 700 and 250 G respectively and were chosen, as they resulted in a slope of the linear regression close to unity (1.029) and a high r 2 value (0.957).The mean rotational variability in the present model is comparable to that of the N20b model (1459.6 to 1457.71 ppm).However, we note that the threshold approach used in this model favours spots over faculae.While this is not necessarily accurate for the solar rotator in the present work, it leads to the formation of polar spots and spot domination of the stellar variability for the faster rotating stars.
We present maps for different inclinations and nesting degrees of the spot and facula distributions (following the description in Sect.2.2 and the choice of the parameters presented just above) for a star with 8Ω as returned by FEAT with various nesting degrees (see Sect. 2.1) in Fig. 5 and Fig. 6, respectively.Clearly visible in Fig. 5 are the polar spots.

Light curves
Using the model described in Sect.2.1, we generated a synthetic 11-year-long cycle for each rotation rate considered in this work.We focus on the four years around the activity maximum for two reasons.Firstly, we aim to explain the upper envelope of the variability distribution of stars as a function of the rotation period (McQuillan et al. 2014).Secondly, we model a single activity cycle, without overlapping cycles preceding and following it.Modelling cycle overlap is beyond the present scope, and it would affect the activity level around cycle minima, which we thus exclude from the analysis.
We note that the LCs shown in this section are 180-dayslong snippets of the full LCs and were simply chosen for demonstration.
Firstly, we consider the LCs of stars that are observed equator-on (i = 90 • ). Figure 7 displays the detrended LCs for the non-nested case (p = 0), for different rotation rates.One can see that the faster the star rotates, the higher its amplitude of variability (mainly due to the shorter rotation periods -the lifetimes of the active regions change far less) which is a consequence of the activity-rotation scaling (see Step 2 from Sect.2.1).The small bumps in the LCs that are mostly seen in the case for the Ω = 1Ω are a result of firstly having all BMRs emerge at their maximum size in the SFTM.Since the emergence rate is generally low and not many other BMRs are present, this affects the LCs more severely in this low activity case.Secondly, and lastly, due to the nature of the threshold approach outlined in Sect.2.2 flux superposition (in case of same polarity encounters) and cancellation (in case of opposite polarity encounters) leads to the possibility of flux switching from being associated with faculae to spots (or spots to faculae) rapidly.Similar effects for the solar rotator will be visible in the following plots as well for the very same reason.
Figure 7 shows that not only amplitude but also the shape of the LCs strongly depends on the stellar rotation rate.For the case of the solar rotation (Ω = 1Ω , top panel of Fig. 7), most of the individual dips in the LCs correspond to transits of different active regions (since active regions evolve on timescale shorter than the solar rotation).As active regions emerge randomly in time, the LCs appear quite irregular.In contrast, the LCs of the more rapidly rotating stars show gradually more regular patterns in brightness variations (see the lower three panels in Fig. 7).This is because active regions on such stars can survive several rotation periods.Furthermore, the large amount of BMR emergences over mid-to high latitudes with large tilt angles leads to the formation of polar spots at about 4Ω , with prominent polar spot caps being present for the 8Ω rotator (see Fig. 5).The formation of polar spots for stars with those rotation rates are consistent with Paper I and Doppler-imaging observations.We note that the polar spots in our simulations turn out to be non-axisymmetric unipolar caps.Their overall structure is rather stable, because their decay is compensated by the magnetic flux coming from the new emergences, as long as the activity level and the BMR polarity orientations are sustained.As a next step, we consider models with active-region nesting.
We first consider the effect of the active-longitude nesting (see Sect. 2.1). Figure 8 shows LCs synthesised with AL nesting of 70% (i.e.p = 0.7).The overall shape of the LC for 1Ω is still rather irregular, compared to the faster rotators, due to the low emergence frequency.With increasing rotation rate, dips related to BMR transits not only occur at a separation of the rotation period, but also at half of the rotation period.In addition to a change in the morphology of the LCs, the variations are amplified with respect to the corresponding non-nested case at each rotation rate (black solid lines versus coloured lines).
Figure 9 gives the most extreme case of AL nesting, where we assume that all BMRs emerge in one of the two active longitudes (e.g.p =1).Clearly, both dips (at one and half-rotational period interval) occur in all of the cases and the LC amplitudes are further augmented.The LCs in Fig. 9 additionally show that the two peaks have almost the same amplitude for the faster rotators.In case of the solar rotator and the low emergence frequency, temporarily asymmetries between the two AL can arise in case of the emergence of a large BMR.However, with increasing emergence frequency, the two ALs become more and more symmetric and hence the amplitude of the two dips due to the ALs as the star rotates are, to a large extent, comparable.
We now focus on the behaviour of the LCs in the FN mode.In Fig. 10 we present the LCs produced with p = 0.7.One can see that the amplitudes of the variability increase for all four shown rotation periods.The light curves also appear more regular than those calculated with p = 0. We increase the nesting even further to p = 0.99 in Fig. 11.The change in the LCs with respect to the non-nested case is remarkable.The amplitude of the LCs is enhanced for all cases and the runs with p = 0.99 exhibit regular patterns, even in the solar case (top panel of Fig. 10).Interestingly, in the displayed LCs, not only dips with separation of the rotation period, but also half-of the rotation period (e.g.most prominently seen in the 4Ω star) appear.This is interesting, as the dips at half of the rotation period, appear and disappear and are not constant, compared to the AL case.We will discuss this further in Section 5. Next, we consider the inclination effect on the LCs.For demonstration, we limit ourselves to the non-nested cases with different rotation rates.In Fig. 12 we show the time-span of 0-90 days from Fig. 7 for inclinations of 90 • , 60 • and 30 • .In the given timespan, for 1Ω , the amplitude of the variability de-   creases with decreasing inclination.The shape of the transits also change.For 2 and 4Ω , the LC amplitudes decreases for the inclinations shown here as well.Interestingly, the situation changes for the 8Ω case.The amplitude increases from i = 90 • to 60 • and then decreases from i = 60 • to i = 30 • .Also the amplitude of variability observed at i = 30 • is larger than that at i = 90 • .These inclination dependencies can be explained with the help of the magnetic field maps in Fig. 5.For the case of 1Ω all regions emerge within ± 30 • around the equator.The Coriolis effect gets stronger with increasing rotation rate, so that for the 8Ω star, the BMRs can emerge at latitudes up to ± 70 • , while a latitudinal belt free of active regions opens up around the equator between ±20 • latitudes (see Sec. 5).For i = 90 • , the  high-latitude spots appear close to the limb, where their effect on the brightness is significantly reduced due to the foreshortening.With decreasing inclination, the majority of spots shift towards the centre of the visual disc and their effect on brightness increases, so that at intermediate inclinations, the variability reaches a maximum (see Figs. 13 panel a), before it starts decreasing again, as the spots move towards the limb of the visible disc once more.

Comparison to observations
In the following, we compare our model to observational records of Sun-like stars obtained by the Kepler telescope.Since our model builds on the solar paradigm, we select stars with nearsolar effective temperatures between 5500-6000 K and surface gravity log g ≥ 4. We express the variability through the quantity R var , first introduced by Basri et al. (2010Basri et al. ( , 2011)).We use a slightly modified version of the range R var : we compute the difference between the 95th and 5th percentile of the sorted differential intensities.Even for the latest Kepler data reduction (DR25) it occurs that some quarters still include instrumental systematics that might influence the range.We then compute the median absolute deviation (MAD) of all R var values, and remove those quarters, that deviate more than 6x the MAD.Afterwards the median is taken of all remaining quarters, which we call R var .Since the SFTM returns instantaneous values with 6-hour cadences, we take every 12th data point (the Kepler long cadence is ≈30 min) and compute R var for this unbinned time series.All computations are based on the latest Kepler data release (DR25) using the PDC-MAP light curves.We also cross-checked our calculated variabilities represented through R var , with the metric R per used by McQuillan et al. (2014) and found very good agreement.
Before computing R var of the models, we need to add noise to the light curves.We follow the same strategy as described in detail in Reinhold et al. (2020).Here, we multiply the noise with a factor of √ 3h/30min = √ 6 to account for the different time bins.For each inclination, 1000 noise realizations are considered, and the mean and standard deviations are computed.This has been done only for the 1Ω p = 0 case, as we found that the noise level is significantly lower than the actual stellar variability in all other cases.
In Figs. 13 and 14 we show the comparison between the calculated R var values of our simulated stars for various degrees of nesting in the AL and the FN mode, respectively, and the two samples of Kepler stars shown as grey ( McQ14) and black (S21) dots.Evidently, if we do not include any nesting (p = 0), our calculated variabilities underestimate the bulk of the observed variability amplitudes, especially for the faster rotators.With increasing nesting level (Fig. 13 and 14 b, c, d, and e), R var increases and the values move towards the upper edge of the distribution.Interestingly, p = 0.99 in the FN mode (Fig. 14 e) overestimates the variability of the solar rotators but leads to similar variability values as the upper envelope of the variability distributions of the faster rotating stars.We note that, while the numbers of BMR emerging is the same for a given rotation rate between the non-nested and nested cases, with higher nesting degree, the spot disc area coverage increases due to the formation of spots by flux superposition.As a consequence, the spot area coverage is not preserved in contrast to the approach presented in Işık et al. (2020) and the nesting has a stronger effect on variability in our model than in Işık et al. (2020).We will elaborate on this point further in Sect. 5.
For a more detailed comparison between the simulations and the observations, we bin the distribution of observed variabilities.Namely, we compare the variabilities returned by our model for a star rotating X times faster than the Sun, with variabilities of a sample of stars with periods in the range [23/X, 27/X] days from the Kepler samples.This comparison is shown in Fig. 15.The histograms in grey and black display the range of variabilities within each of the rotation period bins.We note that the number of stars within each rotation bins decreases from panel a to d. Similarly to Figs. 13 and 14, Fig. 15 shows that the calculations with p = 0 clearly underestimate the variability for the 2Ω and 4Ω sub-samples (while the small number of Kepler stars in the 8Ω rotation bin makes interpretation of panel d rather difficult).
For the stars with near-solar rotation periods (i.e the 1Ω case), there is a substantial difference between the stars in McQ14 and S21.The sample of S21 contains many more stars with variabilities lower than the solar variability at the maximum of cycle 22 (blue curve in the left panel of Fig. 15), while both samples roughly contain the same number of stars substantially more variable than the Sun (i.e. the high-variability tail).This result again raises the question of whether the Sun could also become as variable as those stars in the high-variability tail (Reinhold et al. 2020).Moreover, it agrees with the conclusions drawn in Reinhold et al. (2020Reinhold et al. ( , 2021) ) that the solar variability is not unexpectedly low but that the rotation periods of many stars with At the same time, the p = 0.99 FN calculations in Fig. 15 lie towards the upper bound of the variability distribution of the faster rotators, while highly overestimating the variability of stars with near-solar rotation periods.One can see that different nesting modes can lead to similar values of median R var , especially if the inclination of the stellar rotation axis is not known.This degeneracy might be lifted if alternative metrics are used.We discuss further metrics of characterising stellar variability including the morphology of the LCs in Sect. 5.

Discussion
The premise of this work was to extend the solar paradigm to model the distribution of magnetic features on stars rotating faster than the Sun to use these to consecutively calculate the stellar light curves and the amplitude of the variability.Figure 15 shows that while our calculations with a rather high degree of nesting can reproduce the bulk of variabilities in the Kepler sample they do not catch the maximum of the variability distribution.This might be because we have used the emergence frequency of BMRs for solar cycle 22 as reference cycle and scaled the emergence frequency as a function of the rotation period based on that cycle (i.e. a star rotating twice as fast as the Sun exhibits twice as many BMR emergences).First, solar cycle 22 does not represent the maximum level of activity the Sun is capable of (see, e.g.Usoskin 2017).Second, the activity-rotation scaling itself is rather approximate (see Işık et al. 2018, for more details).We refer here to Işık et al. (2020), who considered the effect of activity increase on the variability with a considerably simpler model limited to stars with near-solar rotation rates.
Empirical studies suggest that faster rotating stars have activity cycles shorter than the present Sun (Böhm-Vitense 2007).In spite of this, as the main purpose of this study is to model stellar variability on the rotational timescale, we follow the approach of Paper I and assume that the solar temporal profile of activity remains unchanged for the faster rotators, for simplicity.However, we only consider the maximum level of activity (in a four-year window during activity maximum, see Fig. 1).Thus, our light curves (and the resulting peak-to-peak variability) correspond to the maximum of the stellar activity cycle and are largely unaffected by the temporal profile of the underlying activity cycle.
Another parameter space that we have not taken into account in this study is the effect of stellar fundamental parameters (e.g. the effective temperature and the stellar metallicity).While we limited our sample of Kepler stars to a temperature range of 5500-6000K, it contains stars with a rather broad range of metallicity values.Witzke et al. (2020) have shown that changing the metallicity for a star with solar level of activity enhances the possibility of recovering its rotation period, as the star moves out of the compensation regime between facular and spot contribution on the rotational timescale.It is less intuitive, however, how the change in metallicity will affect the calculations presented in this paper.According to our calculations, the rotational variability in the rapid rotators is primarily driven by spots, yet the spot com-ponent is less affected by the change in the metallicity than the facular component (see Witzke et al. 2018, for a more coherent discussion).However, metallicity might have an impact on the activity of a star and on the surface distribution of magnetic features (Amard & Matt 2020;Amard et al. 2020;See et al. 2021) -an effect that we deem to be outside the scope of the present study.
The latitude of emergence calculated by flux-tube simulations (see Sect. 2 and Paper I) puts a well-defined lower limit for the latitude of emergence, which becomes more visible for fast rotators (see Fig. 1).This is a consequence of the inward-directed Coriolis force in the rotating frame, consistently acting on rising flux tubes having a pro-grade azimuthal flow (i.e., they rotate faster than their locality).This leads to a well-defined minimum latitude of emergence, corresponding to the minimum non-zero latitude of injection near the equator at the base of the convection zone.Whether such a latitudinal gap around the equator occurs on rapid rotators is unknown (see S , enavcı et al. 2021, for a discussion), but one would expect that stochastic effects (e.g.convection) can induce scatter around the minimum latitude of emergence.We reckon, however, that this additional scatter will not affect the photometric variability to a large degree.
When comparing different rows in Fig. 5, one interesting feature becomes apparent: the spot area coverage is larger for the free nesting mode than for the non-nested case and active longitude nesting.This is because the proximity of neighbouring magnetic flux elements leads to a high possibility of same-polarity encounters.Magnetic flux, which accumulates this way, leads to the formation of spots.We note that the possibility of spontaneous spot formation via flux superposition has been detected in numerical simulations (Kitiashvili et al. 2010), but so far such a formation has not been observed on the Sun, probably due to its relatively low activity level and a small degree of nesting of solar magnetic features.
Figures 7-11 show that nesting affects not only the amplitude of the variability, but also the morphology of the LCs (see also discussion in Sect.3), in parallel with the computations by Işık et al. (2020).For example, in the case of AL nesting, dips in the LCs mainly occur each half of the stellar rotation period.Basri & Nguyen (2018) have introduced the "single/double ratio" (SDR) metric, which provides information about the ratio of time a star spends in single-or double-dip modes (i.e. its LC shows one or two peaks per period, respectively).The SDR was proposed to be an effective metric for characterizing stellar LCs.SDR as well as other metrics might be useful for testing the models and constraining the distribution of stellar magnetic features.Indeed, Fig. 15 shows that both AL-and FN-type modes of nesting can result in the same amplitude of variability, albeit at different nesting degrees.Any metric, such as R var , which is used in the present work, do not take into account LC morphology, as they only measure the peak-to-peak variability.Using metrics sensitive to the morphology of the LCs will help to distinguish between various modes of nesting.This will be addressed in a forthcoming publication.
Finally, our calculations are based on the assumption that the size distribution of emerging spots is the same on solar-like stars and the Sun, i.e. it does not depend on the rotation rate and stellar activity level.Specifically, the source term in Jiang et al. (2011) used in this study (see Step 1 in Sect.2.1) is based on the sunspot size distribution during solar cycle 22.This neglects that the spot size distribution might depend on the activity level (see, e.g., Solanki & Unruh 2004;Krivova et al. 2021).We note that increase of spot sizes in the source term of SFTM model would simultaneously amplify the variability of a star and make its LC more regular.Thus, the change of distribution of emerging spots might be another mechanism capable of explaining Kepler observations (together with nesting investigated in this study).Since the lifetime of sunspots depends on their sizes, a change in the size distribution of spots would also lead to a change in their lifetimes, which would have a very direct effect on the variability amplitude and the LC statistics, with longer-living spots making the LC more regular.In that sense extending the lifetime of a spot should have a similar effect as stronger nesting.This is an important parameter whose influence will be studied in a future publication.

Conclusions
We coupled the model for the emergence and surface transport of magnetic flux in Sun-like stars (Işık et al. 2018, Paper I) with a model for calculating stellar brightness variations (partly based on the approach presented in Nèmec et al. 2020b).This allowed us to compute light curves of stars with rotation rates between 1 and 8Ω as they would be observed by Kepler at different inclination angles.Following up on the findings of Işık et al. (2018) and Işık et al. (2020), we investigated the impact of active-region nesting on the light curves and, in particular, the amplitude of the variability.
We compared the output of our model to the observed variabilities of Kepler stars in the temperature range 5500-6000K.In particular, we aimed at explaining the dependence of the amplitude of the variability on the rotation rate.Recently, Işık et al. (2020) showed that the model without nesting underestimates the variability of stars with known near-solar rotation periods (see also Reinhold et al. 2020).We found that the same is true for stars rotating faster than the Sun.Our runs without nesting dramatically underestimate the stellar variability at all rotation rates.
We showed that the observed dependence of Kepler variabilities on the rotation period, for stars with detected rotation periods, can be explained by an increase of nesting degree with the rotation rate, in parallel with increasing activity level.As both modes of nesting used in this work led to similar levels of variability for stars with different rotation rates, we plan to further investigate the use of metrics that consider LC morphologies (instead of the peak-to-peak variability), to retrieve more information regarding the surface magnetic activity of stars.Additionally, we plan to include the effects of different rotation-activity relationships, cycle length and of stellar fundamental parameters (i.e.effective temperature and metallicities), on the variability.The applications of the FEAT model extend beyond stellar photometric variability, which is presented in this work.The model has been adapted to study the astrometric jitter introduced by stellar magnetic activity (Sowmya et al. 2021(Sowmya et al. , 2022) ) and Doppler imaging (S , enavcı et al. 2021).It can also be used to study the magnetic contamination of high and low-resolution transmission spectra (see, e.g.Rackham et al. 2018Rackham et al. , 2019;;Dravins et al. 2021;Rackham et al. 2022).

Fig. 1 .
Fig. 1.Butterfly diagrams of BMR emergence for 1Ω left panel and 8Ω right panel.The colour-bar gives the longitudes of emergence and the vertical dashed lines indicate the 4 years of the 11 year cycle considered in our brightness variation calculations.No amount or nesting is included here.

Fig. 2 .
Fig. 2. Snapshot of magnetic field distribution for 1Ω (top row) and 8Ω (bottom row) as returned by the FEAT model at various inclination angles around maximum of the activity cycle.First column shows i = 90 • , second column i = 57 • , third column i = 30 • , and forth column i = 0 • .

Fig. 3 .
Fig. 3. Butterfly diagrams of BMR emergence for 1Ω , with p = 0.7 in the different nesting modes considered in this work.Top row: AL nesting, bottom row: FN nesting.Left panels are time-latitude diagrams with the colourbar indicating the longitudes, right panels are timelongitude diagrams, with the colourbar indicating the latitude.

Fig. 4 .
Fig. 4. Comparison of R var (in ppm) as returned by the N20b model and the present model.All values are given in ppm.Black solid line gives the linear regression, whereas the blue dashed line is the 1-to-1 correspondence between the two models

Fig. 7 .
Fig. 7. Synthetic light curves (LCs) for stars with different rotation rates as they would be observed in the Kepler passband at an inclination of i = 90 • .Shown are not-nested cases (p = 0) with rotation rate values of 1Ω (blue), 2Ω (orange), 4Ω (green), and 8Ω (red).

Fig. 8 .
Fig. 8. Similar to Fig. 7, where the black curves represent the calculations with p = 0 for each rotation rate and the coloured curves represent those with added AL-type nesting at p = 0.7.
2 using the updated fundamental parameters of Kepler stars from Berger et al. (2020).Next, we consider stars with known rotational periods using the data sets of McQuillan et al. (2014) (hereafter McQ14) and Santos et al. (2021) (hereafter S21).These constraints lead to a sample of 6,228 stars when using McQ14 rotation periods and 11,493 stars when using S21 rotation periods.

Fig. 13 .
Fig. 13.Comparison of R var as a function of the rotation period between stars with effective temperatures with 5500-6000 K and log g > 4.2 with detection rotation periods from McQ14 (grey dots) and S21 (black dots) and the modelled stars.Each panel includes different nesting probabilities p in the form of active-longitude (AL) nesting.The different colours indicate the inclination of the modelled stars.

Fig. 15 .
Fig. 15.Dependence of R var on the inclination.Each panel represents different realisations of the BMR emergence for a given rotation rate.The grey histograms are drawn from the McQ14 stars and the black histograms from the S21 stars, with the limitations of [23/X, 27/X] days for a given rotation rate X in units of the solar rotation (i.e. for X=2, meaning stars with 2Ω , the stars have rotation periods between[11.5,13.5]days).