Issue 
A&A
Volume 620, December 2018



Article Number  A177  
Number of page(s)  14  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201833393  
Published online  14 December 2018 
Forward modelling of brightness variations in Sunlike stars
I. Emergence and surface transport of magnetic flux
^{1}
MaxPlanckInstitut für Sonnensystemforschung,
JustusvonLiebigWeg 3,
37077
Göttingen,
Germany
email: isik@mps.mpg.de; solanki@mps.mpg.de
^{2}
Feza Gürsey Center for Physics and Mathematics, Boğaziçi University,
Kuleli 34684,
Istanbul,
Turkey
^{3}
School of Space Research, Kyung Hee University,
Yongin,
GyeonggiDo,
446701,
Republic of Korea
Received:
8
May
2018
Accepted:
15
October
2018
Context. The latitudinal distribution of starspots deviates from the solar pattern with increasing rotation rate. Numerical simulations of magnetic flux emergence and transport can help model the observed stellar activity patterns and the associated brightness variations.
Aims. We set up a composite model for the processes of flux emergence and transport on Sunlike stars to simulate stellar brightness variations for various levels of magnetic activity and rotation rates.
Methods. Assuming that the distribution of magnetic flux at the base of the convection zone follows solar scaling relations, we calculate the emergence latitudes and tilt angles of bipolar regions at the surface for various rotation rates, using thinfluxtube simulations. Taking these two quantities as input to a surface flux transport (SFT) model, we simulate the diffusiveadvective evolution of the radial field at the stellar surface, including effects of active region nesting.
Results. As the rotation rate increases, (1) magnetic flux emerges at higher latitudes and an inactive gap opens around the equator, reaching a halfwidth of 20° for 8 Ω_{⊙}; and (2) the tilt angles of freshly emerged bipolar regions show stronger variations with latitude. Polar spots can form at 8 Ω_{⊙} by accumulation of followerpolarity flux from decaying bipolar regions. From 4 Ω_{⊙} to 8 Ω_{⊙}, the maximum spot coverage changes from 3 to 20%, respectively, compared to 0.4% in the solar model. Nesting of activity can lead to strongly nonaxisymmetric spot distributions.
Conclusions. On Sunlike stars rotating at 8 Ω_{⊙} (P_{rot} ≃ 3 days), polar spots can form, owing to higher levels of flux emergence rate and tilt angles. Defining spots by a threshold field strength yields global spot coverages that are roughly consistent with stellar observations.
Key words: stars: activity / stars: solartype / starspots / stars: magnetic field / methods: numerical / magnetohydrodynamics (MHD)
© ESO 2018
1 Introduction
The advent of spaceborne photometry, which is designed primarily to detect extrasolar planets, has made shortterm stellar activity signals measurable (Koch et al. 2010). Brightness variations in Sunlike stars on timescales ranging from hours to decades are driven by convective flows and magnetic fields at the surface, while the observed patterns of variability also depend on the stellar rotation and related quantities (Aigrain et al. 2004; Shapiro et al. 2014, 2017).
Dark spots and bright faculae observed on the Sun are formed by relatively strong magnetic flux concentrations threading the photosphere. To investigate radiative variations driven by magnetic activity and rotation on Sunlike stars, it is important to numerically simulate the distribution and the evolution of surface magnetic fields at large scales.
The latitudinal distribution of starspots is one of the extensively studied features of magnetic activity on rapidly rotating cool stars, via Doppler and ZeemanDoppler imaging (Strassmeier 2009; Donati & Landstreet 2009). For Sunlike stars in particular, the mean latitude of cool spots and magnetic fields has been observed to increase with the rotation rate. Spots near or at the rotational poles are observed on rapidly rotating Sunlike stars with rotation periods below about 3 days (Jeffers et al. 2002; Marsden et al. 2006; Järvinen et al. 2007; Waite et al. 2015).
As a possible explanation for polar spots, Schüssler & Solanki (1992) suggested a mechanism in which the action of the axially directed Coriolis force in the rotating frame of a rising flux tube becomes comparable to or even dominates the outward buoyancy force for sufficiently fast rotation. Schüssler et al. (1996) demonstrated how faster rotation deflects rising tubes towards high latitudes on Sunlike stars, and Granzer et al. (2000) obtained latitudinal distributions for a range of stellar masses.
A surface mechanism which can also contribute to the formation of polar spots was suggested by Schrijver & Title (2001) using a randomwalk model of flux dispersal and transport. In this model, a strong highlatitude concentration of magnetic fields of a single polarity was formed by the trailing polarities of bipolar regions when the flux emergence rate was 30 times the solar value. A ringlike structure was formed around a polar spot of opposite polarity, which had a much longer decay time than what a diffusion model would infer.
Another attempt was made by Işık et al. (2007b), who used a diffusive surface flux transport (SFT) model subject to observed stellar surface differential rotation estimates and stellar radii. They suggested that the long lifetimes of polar spots can be explained by the midlatitude emergence of a large bipolar region with a large tilt angle, presumably owing to strong action of the Coriolis effect. In both Schrijver & Title (2001) and Işık et al. (2007b), a largely unipolar polar spot was produced by the diffusion of flux from the followerpolarities of bipolar magnetic regions (BMRs) at low or medium latitudes.
Holzwarth et al. (2006) also used an SFT model to show that a sufficiently fast poleward meridional flow at the surface can also lead to polar spots, but with intermingled polarities. In an attempt to reproduce the flipflop activity variations that have been reported for some active stars (Berdyugina 2005), Elstner & Korhonen (2005) applied a nonaxisymmetric meanfield dynamo model, which produced an asymmetric distribution of magnetic field concentrations near the rotational poles.
Schüssler et al. (1996) modelled the latitudinal distribution of emerging flux on solarmass stars, and Işık et al. (2007a, 2011) studied the combined effects of the dynamo, emergence, and SFT processes on stellar cycles. It is known that the average inclination of BMRs with respect to the eastwest direction (the tilt angle) is crucial in maintaining the solar activity cycle (Baumann et al. 2004; Jiang et al. 2010). However, a systematic modelling of Sunlike stars with physically consistent computation of tilt angles and the subsequent SFT here and after has not yet been undertaken.
To better understand the implications of the rotationactivity relation not too far from the solar regime and to provide a testing platform of forward modelling for photometric reconstructions, here we construct a modelling framework. It incorporates the observed properties of the solar cycle, the effects of rotation on rising flux tubes in the convection zone, and SFT. In this paper, we present the method and discuss the impact of increasing rotation rate and activity on the surface patterns of magnetic flux. A following paper will be devoted to syntheses of light curves based on the magnetic flux emergence (Sect. 2.2) and transport models (Sect. 2.3) presented here. The purpose of the model presented here is restricted to reproducing brightness variations of Sunlike stars in highprecision data such as from Kepler (Koch et al. 2010) and TESS (Ricker et al. 2015), particularly from very short timescales up to several stellar rotations. We do not, however, rule out a later extension to also model spectroscopic (e.g. radial velocity) and spectropolarimetric changes by including horizontal magnetic fields. In Sect. 3, we present results from our scaled stellar models, namely, temporal and latitudinal patterns of activity. The limitations of the model and the relevance of our results with respect to the observations of active Sunlike stars is discussed in Sect. 4.
2 Model setup
To model the emergence and the evolution of the surface magnetic field, we first generated an 11yr semisynthetic sunspot group record based on the statistics of solar cycle 22 (Sect. 2.1.1). With increasing stellar rotation rate, we linearly scaled the flux emergence rate, in agreement with the observed rotationactivity relation. Stronger solar cycles are known to show a tendency for higher mean emergence latitudes (Solanki et al. 2008; Jiang et al. 2011). We therefore used an empirical relation to extrapolate the mean latitude of emergence to higher activity levels (Sect. 2.1.2). We assumed that the resulting distribution represents the butterfly diagram of fluxtube eruptions at the base of the convection zone of a G2Vtype star with a given activity level. Using fluxtube simulations, we then computed the emergence latitudes and tilt angles of rising loops for a given stellar rotation rate (Sect. 2.2). A schematic presentation of the algorithm is given in detail in Appendix B.
The latitudinal distribution of emerging flux is therefore determined by two independent effects in our approach: (i) the base latitudinal distribution changes with the flux emergence rate (extrapolation of the empirical solar relation between the activity level and the mean latitude); and (ii) the rotation rate affects the action of the Coriolis force, which tends to deflect rising flux tubes poleward.
The resulting starspot emergence records are then used as input to the SFT model to obtain the evolution of surface magnetic flux (Sect. 2.3). In the following subsections we describe the various steps in this chain in greater detail.
2.1 Synthetic starspot group records
2.1.1 The input solar cycle
Using the RGO/SOON records between 1700 and 2010, Jiang et al. (2011, hereafter JCSS11) found statistical relationships between various characteristics of sunspot groups. They generated a semisynthetic spot group record (SGR) based on these relationships. We adopted the solar SGR generation procedures described in JCSS11 and Cameron et al. (2016, hereafter CJS16; see their Sect. 2.2) to simulate a relatively strong, artificial cycle lasting about 11 yr. We describe the details of this procedure in Appendix A.
In summary, the synthetic cycle module provides the following quantities related to sunspot group emergence as a function of time (cycle phase): the total number of sunspot groups (Eq. (A.1)); the emergence times and latitudes, which are determined by the mean latitude of sunspot groups (Eqs. (A.2) and (A.3)); the width of the latitudinal spread (Eq. (A.4)); random longitudes (nesting is introduced in Sect. 2.1.3); sunspot group areas drawn from a composite empirical distribution (Eq. (A.5)).
Instead of using the tilt angles synthesised within the JCSS11 model from empirical relationships, we use the results of the fluxtube simulations to be described in Sect. 2.2. For the solar case, the results are similar to the tilt angle relationship assumed by JCSS11, while for more rapidly rotating stars very different tilt angle distributions are obtained.
Having defined the base model of the solar cycle, we now describe how we determine the flux emergence rate and the mean latitude of flux eruption, for more active Sunlike stars.
2.1.2 Scaling the cycle on more active stars
To scale the stellar cycle amplitude, S_{⋆}, we define , where S is in units of the maximum of the annual running mean of the sunspot group number, S_{⊙} is set to 156, which is the observed value for Cycle 22 (see Appendix A.1), and is then a relative measure of the emergence frequency (or rate) of BMRs; for simplicity, this quantity will be called the emergence frequency or flux emergence rate (in solar units) throughout the rest of the paper, as they are proportional quantities. To evaluate , we focused on two extreme cases: (i) a constant , to isolate the rotational effects on the (otherwise solar) emergence pattern; and (ii) a linear scaling , following the observed proportionality Bf = 50v_{eq} between the field strength (Bf, where f is the filling factor) and the equatorial rotational velocity v_{eq} of Sunlike stars, estimated by Reiners (2012). To first order, Bf is proportional to the total magnetic flux on the stellar disc. The scaling therefore provides a rough representation of the rotationactivity relationship of G dwarfs.
The effect that stronger solar cycles start at higher latitudes has already been considered by JCSS11. We extrapolated this effect to higher activity levels () and expressed the mean latitude of flux eruptions from the base of the stellar convection zone, , by modifying Eq. (A.3), as (1)
where we assumed k = 0.022 for , as in Eq. (A.3). We chose k = 0.014 for to limit the maximum latitude of the input cycle to 73°. This roughly corresponds to the highlatitude edge of the main region of instability of flux tubes (see Sect. 2.2.1).
Rather than assuming the surface emergence record described above to represent the distribution at the base of the convection zone, we corrected the eruption latitudes at the base of the convection zone for the weak poleward deflection of rising flux tubes at the solar rotation rate (see Appendix B).
2.1.3 Nesting of active regions
It is known that sunspot groups tend to emerge within “nests” of activity (e.g. Castenmiller et al. 1986). To simulate the observed clumping of activeregion emergence, we modified the longitudes and latitudes obtained in the previous steps by setting a probability that a given BMR would be part of a nest, which we call nesting probability (Appendix C). When this effect is included, the resulting cycle variation of loworder multipoles such as the equatorial dipole as well as the open flux better represent the observed variations.
Figure 1 shows the resulting solar reference butterfly diagram with , and a nesting probability of 70%. The only significant difference from the synthetic cycle22 butterfly diagram of the JCSS11 model is the clustered emergence pattern.
Fig. 1 Semisynthetic emergence record (SGR) for solar cycle 22, using emergence latitudes and tilt angles resulting from the remapping procedure for the solar rotation and flux emergence rate . The colour bar shows the tilt angles with a rather high saturation level, so that the tilt angles can be compared with faster rotating cases. The plot on the left shows the histogram of emergence latitudes. 

Open with DEXTER 
2.2 The rise of flux tubes in the convection zone
We model the latitudinal distribution of emerging magnetic flux on Sunlike stars using numerical simulations of buoyantly rising magnetic flux tubes. We adopt the thin flux tube approximation (Spruit 1981) and model flux tubes leading to sunspot groups as initially toroidal flux rings in mechanical equilibrium (Caligari et al. 1998) in the form given by FerrizMas & Schüssler (1995) and Caligari et al. (1995). At the equilibrium state, the flux tube is assumed to lie in the stably stratified overshoot region at the base of a solarlike, nonlocalmixinglength convection zone model, which is adopted from Skaley & Stix (1991).
We take the angular rotation rate Ω as a function of radius r and latitude λ, (2)
where Ω_{⋆} is the equatorial rotation rate at the stellar surface, c_{1} = 0.0876, c_{2} = 0.0535, and c_{3} = −0.0182. The base of the convective overshoot region in the stratification model is at r_{0} = 0.724 R_{⊙}; it is taken as the centre of the tachocline here. Furthermore, d = 0.075R_{⊙} is the width of the error function defining the thickness of the tachocline and the constants (c_{1}, c_{2}, c_{3}, d) are chosen such that Eq. (2) closely mimics helioseismic inversions of solar internal rotation (Schou et al. 1998). Observational studies indicate that surface differential rotation increases rather weakly with the rotation rate as , where n was estimated to be 0.15 by Barnes et al. (2005), and 0.2 for G stars by Balona & Abedigamba (2016). For simplicity, we set Δ Ω_{⋆} = ΔΩ_{⊙} in both the radial and latitudinal directions. The factor and the term in Eq. (2) account for keeping the same (solar) differential rotation rate in both the radius and latitude, as Ω_{⋆} increases.
2.2.1 Initial properties of flux tubes
Following FerrizMas & Schüssler (1995), we solve the sixthorder dispersion relation for linear perturbations of a toroidal flux ring in mechanical equilibrium, taking the thermodynamical quantities from the stratification model and the rotation profile from Eq. (2).
Figure 2 shows the linear stability diagrams for toroidal flux tubes in the middle of the overshoot region (r_{mid} = 0.728R_{⊙}) as a function of the initial latitude λ_{0} and the field strength B_{0}, for different rotation rates and the same Sunlike stratification. In light and darkshaded regions, the fastestgrowing wave mode has an azimuthal wavenumberof m = 1 and m =2, respectively.The radial location is about 5000 km beneath the base of the convection zone (the term “base” signifies the depth at which the convective heat flux changes its sign). The red dots in Fig. 2 show λ_{0} and B_{0} of flux tubes chosen for the nonlinear simulations (Sect. 2.2.2) with 2° latitudinal steps, all corresponding to a characteristic linear growth time of 50 days. This growth time ensures that the tubes are sufficiently close to the onset of instability and that the Joy’s law resulting from the simulations for matches well with the solar observations. We did not consider the islands of instability because (1) the corresponding growth time is not reached in the lower latitude island; and (2) the highlatitude island is not reached by the input butterfly diagram in any case, except for . To be conservative, we preferred to limit the simulations to the main region of instability by decreasing the factor k in Eq. (1) for .
The maximum latitude at the base is 73° for (Sect. 2.1.2). To cover the entire latitude range in the input solar cycle for all flux emergence rates, we set up fluxtube simulations with initial latitudes at the base of the convection zone up to λ_{0} = 73° for , to move from step I to II. This is the same maximum latitude as for the input model at the surface because the flux tubes rise almost radially for , especially at high latitudes (Fig. 3).
The stability diagrams in Fig. 2 show that the onset of instability is shifted towards higher field strengths for higher rotation rates, by a factor of about 3 for , compared to the solar value, owing to the enhanced Coriolis force component directed towards the rotation axis, which has a stabilising effect. The dynamics of the tube is governed predominantly by the buoyancy, curvature, and Coriolis forces in the rotating frame. The tube radius is relevant only for the drag force, which only weakly affects the resulting emergence properties. For all the simulations, we set the initial crosssectional radius to 2000 km, about 3.6% of the local pressure scale height. With this radius and the initial field strengths corresponding to a growth time of 50 days (Fig. 2), the magnetic flux within a tube is typically about 10^{22} Mx for and 3 × 10^{22} Mx for ^{1}. The fluxtube rise simulations serve only to obtain emergence latitudes and tilt angles, which are only slightly affected by the initial tube radius via the drag force.
Fig. 2 Initial latitudes and field strengths of flux tubes (red diamonds) chosen for fluxtube rise computations, plotted over linear stability diagrams for flux tubes at the middle of the overshoot region, for (panel a), (panel b), (panel c), and (panel d), with differential rotation ΔΩ_{⋆} = ΔΩ_{⊙}. The contour lines denote the growth time in days. The linearly stable regime is shown in white, and the unstable regime is shaded with light grey where the fastestgrowing mode is m = 1, and with darkgrey for m = 2. The red diamonds show the initial tube parameters chosen for the numerical simulations, at a growth time of 50 days. 

Open with DEXTER 
2.2.2 Simulations of flux tube emergence
To model the rise of flux tubes through the convection zone, we carried out simulations starting from the initial parameters mentioned above (Sect. 2.2.1, see Fig. 2). We used the code developed by MorenoInsertis (1986) and extended to threedimensional (3D) geometry by Caligari et al. (1995). It solves the dynamical evolution of a one dimensional (1D), initially toroidal ring of mass elements embedded in a 1D background stratification (Skaley & Stix 1991) in the 3D Lagrangian frame. The mass elements are subject to body forces including the drag force, Lorentz force, and buoyancy, as well as the pseudoforces induced by the Coriolis and centrifugal effects. The evolution of the tube is considered as an isentropic process. The magnetic flux and the integrity of the tube with its closed structure are also conserved. We chose 10^{3} mass elements for the flux tube, which is initially in mechanical equilibrium and subject to a linear combination of azimuthally periodic spatial perturbations with wave numbers in the range 1 ≤ m ≤ 5. Their magnitudes are 10^{−4} in units of the local pressure scale height, in each of the three dimensions. Unstable tubes experience magnetic buoyancy instability and develop loops that rise up to a heliocentric radial distance of about 0.98R_{⊙}. At this point, the simulation halts owing to the ambient pressure scale height becoming comparable with the tube diameter, violating the thintube criterion. We roughly define this stage as the “emergence” of the loop, though the loop is still under the surface. In general, the fastestgrowing azimuthal wave mode in the nonlinear simulations is consistent with the prediction of linear stability analysis (Fig. 2). When m = 2, two buoyant loops form with a 180degree phase difference, and one of them emerges before the other. The simulations are stopped at this point due to the thinfluxtube criterion, so it is not possible to track the other emergence at the opposite longitudinal hemisphere. Therefore, our simulations may be somewhat underestimating the amount of magnetic flux that emerges at the stellar surface when m = 2 is the dominant mode of instability.
The initial (λ_{0}, B_{0}) determined by the growth time of 50 days and the initial radial location are set such that the loops yield a realistic set of emergence latitudes and tilt angles for the Sun. The time from the initial state to the emergence, that is, including the development of the instability, ranges from a few hundred to a thousand days and increases with Ω_{⋆}. We measurethe emergence latitude from the apex of the tube at the end of the simulation. We determine the tilt angle by using the longitudes and latitudes of the preceding and follower legs of the flux loop at 0.97R_{⊙}. The results are roughly consistent with the tilt angles obtained from the slope of the tangent vector at the apex.
Figure 3 shows, as a function of the emergence latitude λ_{e}, the tilt angle α (this dependence is called Joy’s law in solar physics) and the latitudinal deflection λ_{e} − λ_{0} as a function of the initial latitude, for different rotation rates. Simulations were made for the full range of latitudes in the case , including latitudes where no emergence occurs in (smaller diamonds in Fig. 3). This additional latitude range was needed because we scaled the mean latitude of activity with (Eq. (1)), in accordance with the empirical solar model extrapolated to higher flux emergence rates. In general, both λ_{e} and α increase with the rotation rate, owing to enhanced Coriolis force components in the rotating frame.
There are a few abrupt changes in α(λ_{e}), which are also visible in λ_{e} −λ_{0}. To understand the origin of such features, we first eliminated the possibility of a numerical resolution issue. For this, we set up tubes with higher numbers of mass elements (up to 4000), but these runs converged to very similar values of emergence latitudes and tilt angles. Following additional simulations with slightly different initial field strengths, we found that these jumps were robust features. They are shaped by physical effects, that is, they represent regimes where different forces and/or different wave modes become dominant. For instance, the local peak in the tilt angle for at λ_{e} ≃ 15° occurs at the transition of the fastestgrowing mode from m = 1 to m = 2 (see Fig. 2a).
To reveal the nature of the largest jump for , we show in Fig. 4 the emergence phases of flux tubes starting at λ_{0} = 45°, 46°, and 47°, that is, roughly where the jump in the tilt angle at about λ_{e} = 50° occurs (see Fig. 3). The plots clearly depict the transition from the case when a highly tilted part of the tube emerges earlier (Fig. 4a), to the case when a much less tilted part emerges earlier (Fig. 4c). The radial and azimuthal projections of the tubes mark this transition clearly. The relative phase speeds of the two competing loops vary with λ_{0}, such that the lowtilt loop intrudes the hightilt loop above a certain initial latitude of about 46°. The transition from partial to full intrusion is responsible for the lowtilt plateau in Fig. 3a, for λ_{e} between about 50° and 65°. The two loops fully merge beyond λ_{e} = 70°, where moderate tilts of about 22° are reached (see Fig. 3a).
As an independent test, thin flux tube simulations for have been kindly made by M. Weber (priv. comm.) using the code developed by Fan et al. (1993). The only two differences with our setup were that she assumed rigid rotation for the stellar interior and started the tubes in the lower convection zone, where the superadiabaticity was positive. Nevertheless, the resulting latitude dependence of the tilt angle and the latitudinal deflection turned out to be qualitatively similar to our case for . The distribution and amplitude of the abrupt variations in the tilt angle roughly agreed with each other in these two sets of independent simulations.
It is quite possible that both loops forming out of a single flux tube eventually emerge, producing active regions with quite different tilts. Hence, regions with systematically different tilts can coexist at nearly the same latitude on rapidly rotating stars. However, since the thin tube approximation does not allow us to continue the computations further, we have simply used the tilts of the region emerging first.
Fig. 3 Panel a: tilt angle of emerging flux loops vs. their emergence latitude for (black diamonds), (blue asterisks), (green squares), and (red crosses). Panel b: latitudinal deflections (poleward is positive) as a function of the initial latitude at the base of the convection zone. The initial λ_{0} and B_{0} are as inFig. 2. Emergence latitudes for which no activeregionsized BMR emerges for are shown by small grey diamonds (see text). 

Open with DEXTER 
Fig. 4 Geometry of three emerging flux loops with initial latitudes 45° (left panels),46° (middle panels), and 47° (right panels) at the base of the convection zone, for . Upper panels: inner sphere is drawn at 0.72R_{⊙}, and the outer sphere at 0.97R_{⊙}. The parts of the tube that are beneath the outer layer are shaded in grey, whereas the emerged parts are brighter. For each mass element of the tube, the colours denote the crosssectional tube diameter (the redder the larger). Lower panels: latitudinal and radial projections of the tubes, where each mass element is represented by a dot. The horizontal line on the radial profile corresponds to the location of the outer sphere (0.97R_{⊙}), where the tilt angle (α) is measured from the footpoint locations. The red arrows denote the apex of each tube. 

Open with DEXTER 
2.3 Surface flux transport
We used thesolar and stellar SGRs described in Sect. 2.1 as input to the SFT model (see Jiang et al. 2014, for a review), for which we employed the code developed by Baumann et al. (2004). The code solves, with oneday steps, the magnetic induction equation at the solar/stellar surface, where the field is assumed to be purely radial. This allows us to consider the field as scalar, with a sign representing the magnetic polarity. This is a reasonable assumption for the kiloGauss fields on the Sun found in active region plage, solar network, and sunspot umbrae, though the geometry of the highly inclined fields in penumbrae is not taken into account. Because the purpose of this series of studies is to simulate brightness variations, we did not attempt to model the horizontal components of the magnetic field, for which the relationship with brightness variations in Sunlike stars is unclear. Exceptions are spot penumbrae, which can be treated as having a homogeneous brightness at the effective spatial resolution that can be reached in stellar observations. In the current endeavour to model brightness variations, we implicitly assume that larger scale horizontal fields are transients that occur only during flux emergence, and neglect their signature in the brightness distribution on a stellar disc.
2.3.1 Properties of bipolar magnetic regions
In our SFT model, the distribution of the field on the solar surface is represented in terms of spherical harmonic functions, with a maximum degree of l = 64, corresponding to a resolution at the level of supergranular cells on the Sun. The freshly emerged BMRs are defined as two circular regions of opposite polarity, with a fixed upper limit of the field strength, B_{max}. The interpolarity distance (ranging between 3 and 10°) controls thesize of each BMR, as the characteristic radius of each polarity is fixed at 4°. We adopt B_{max} = 374 G, which was determined by Cameron et al. (2010), who matched the variation of the total unsigned magnetic flux from an SFT simulation to magnetographic observations of the Sun from Mount Wilson and Wilcox Solar Observatories.
2.3.2 Surface flows
The surface fields are subject to differential rotation and poleward meridional flow, which are assumed stationary, and follow the same profiles as in Baumann et al. (2004). The latitudinal shear considered in the SFT model is very similar to Eq. (2) at r= R_{⊙}. The difference does not affect the resulting surface flux distributions. The meridional flow reaches a poleward speed of 11 m s^{−1} at midlatitudes and ceases at λ = ±75°. For the effects of smaller scale flows (supergranulation), the SFT model considers the diffusion term in the induction equation, with a horizontal turbulent diffusivity of 250 km^{2} s^{−1} (Cameron et al. 2010) and a radial diffusivity of 25 km^{2} s^{−1} (Baumann et al. 2006).
2.3.3 Initial magnetic field
As the initial condition, we assume a dipolar field reaching G at the poles, whichtakes into account the possibility of stronger initial axial dipole moments for higher levels of flux emergence rate, .
The resulting time series of surface maps of the radial magnetic field are represented in the socalled Carrington frame, where the latitudes of about ± 5° are at rest. We assumed that the SFT parameters were invariant for all sets of , except for the initial field condition.
2.4 Summary of assumptions
Here we summarise our assumptions when modelling faster rotating, more active suns.
 1.
The timelatitude distribution of flux tube eruptions at the base of the convection zone follows the statistical properties of the solar butterfly diagram (Sect. 2.1.1), with the cycle duration set to 11 yr.
 (a)
The mean latitude of eruptions at the base follows the same linear scaling with the cycle strength, as observed on the solar surface (Sect. 2.1.2).
 (b)
The initial field strengths of erupting flux tubes correspond to a constant linear growth time (50 days) of the magnetic buoyancy instability (Fig. 2) at the middle of the convective overshoot layer near the base (Sect. 2.2.1).
 (a)
 2.
The number of eruptions throughout the 11yr cycle scales linearly with the rotation rate, , except for the comparison case (Sect. 2.1.2).
 3.
The stratification and the differential rotation (ΔΩ) in the convection zone are kept the same as in the solar model (Sect. 2.2).
 4.
The SFT coefficients and the largescale flows are the same as in the solar model (Sects. 2.3.1 and 2.3.2).
 5.
The initial surface field is dipolar with a peak unsigned strength (at the rotational poles) equal to G (Sect. 2.3.3).
3 Results
3.1 Butterfly diagrams
Figure 5 shows the butterfly diagrams for six sets of . For 2, 4, and 8, we consider either a solar (, left panels) or a scaled stellar (, right panels) flux emergence rate (Sect. 2.1.2).
For (panels a, c, and e), the systematic increase of both the mean emergence latitude and the tilt angle with is evident, owing to stronger Coriolis acceleration of rising tubes. For (panels b, d, and f), the imposed increase in leads to a further increase in the mean latitude of activity, owing not only to the Coriolis effect, but also to the solarlike scaling of the mean latitude (Eq. (1)), as can be seen in the associated latitude histograms. In addition, the gap of inactivity around the equator widens with increasing . However, the gaps are not very differentfrom each other in cases and 8 (see also Fig. 3b for comparison). The colours represent the tilt angles, as in Fig. 1. The abrupt changes in the tilt angle are visible at some latitudes (cf. Fig. 3a), especially for .
3.2 Variation of surface magnetic activity
The synthetic emergence records presented in Sect. 3.1 are now provided as input to the SFT model. In the first set of simulations, we keep the solar flux emergence rate () and increase the rotation rate (Sect. 3.2.1). In the second set, we increase both quantities with (Sect. 3.2.2).
3.2.1 Solar flux emergence rate ( = 1)
Figure 6 shows the 27day averages of the unsigned flux over the stellar surface and the unsigned mean polar field strength (averaged over both polar caps, which are defined for λ  > 75°) using the same colours as in Fig. 3. Taking λ_{e} and α for with the solar flux emergence rate (, Fig. 6a) leads to a total flux variation which increases only weakly with . This enhancement of the flux is due to the systematic increase in the average tilt angles of emerging BMRs, which weakens flux cancellation over the polarity inversion line of each BMR. While the mean tilt angle ⟨α ⟩ increases from about 5° to 30°, the latitudinal separation between the polarities increases with sin ⟨α⟩. The poleward deflection of rising loops indirectly contributes to the increase (with ) of flux for : the emergence latitudes become more confined to a range where the latitudinal shear is stronger than at the solar activity belts. The strong differential rotation at the activity belt tends to disperse the opposite polarities at an even higher rate, immediately following any BMR emergence. This decreases the cancellation between the opposite polarities within the same BMR. Such selfcancellation accounts for a significant fraction of the initial flux decrease of a BMR.
The isolated effect of faster rotation for (Fig. 5a,c, and e) is more conspicuous for the polar field amplitudes (Fig. 6c), for which the increasing tilt angle is the major contributor. This leads to a larger latitudinal separation between the polarities of a BMR, allowing themeridional flow to transport a larger amount of the following polarity to the pole, while the leading polarity has more time to cancel with the field from other BMRs. The increasing latitudinal deflection produces the opposite effect for the polar field than for the total magnetic flux: by shifting both polarities of each BMR towards higher latitudes, it decreases the efficiency of crossequatorial flux cancellation. However, the effect remains sufficiently weak for up to , for which most BMRs still emerge well below the latitude of the fastest meridional flow (~ 37°). Another consequence of the increase in ⟨α⟩ is that the polar field reverses its polarity and reaches its peak value earlier for increasing .
3.2.2 Flux emergence rate scaled with rotation ( = )
For the set (Figs. 5b,d, and f), the amplitudes of both the total flux and the polar field increase more strongly with than in the case of (Figs. 6b,d). Still, the total flux and the polar field increase at a rate lower than the scaling itself. The polar field reversals occur earlier owing to increasing average tilt angles, but there is no difference in the reversal time from to (8,8). The reason is that a significant amount of BMRs in the case (8,8) emerge within the lowtilt plateau above λ ~ 50° (see Figs. 3a and 5f). In this case, does not provide sufficient contribution to the global (axial) dipole moment to reverse it earlier than in the case (4,4). We note that in (8,8) the initial polar field that has to be reversed is stronger by a factor of two, relative to (4,4), owing to our assumption G.
The rotational dependencies of the total surface flux and the polar flux are shown in Fig. 6e in terms of the time integral of both quantities. The dependence is almost linear for both and , with similar slopes. In the latter case the timeintegrated total flux scales as , with c ≈ 0.9. This is expected from the linearity of the SFT process (Schrijver 2001). The nonlinear dependence of the BMR dipole moments on does not lead to a significant deviation here. Eighttimesfaster rotation leads to a change in the cumulative polar flux of about 2 × 10^{24} Mx, which is about 8% of the corresponding change in the cumulative surface flux. It should be noted, however, that the cumulative polar flux scales nonlinearly with , at a gradually increasing rate, owing to rotationally induced effects. The contributors to the positive correlation between the two quantities are now both the increased latitudinal separation between polarities (increasing mean tilt angle) and the imposed dependence .
Fig. 5 As in Fig. 1, but now for various sets of as indicated at the top of each panel. Left panels: for (solar emergence rate). Right panels: for . 

Open with DEXTER 
Fig. 6 Total unsigned surface flux for (panel a) and (panel b), and the average polar flux density for latitudes poleward of ± 75°, for (panel c) and (panel d). The colours denote , as in Fig. 3. Panel e: comparison of the time integrals of polar flux based on the average polar field and the total unsigned flux, for (diamonds) and (asterisks). 

Open with DEXTER 
3.2.3 Magnetic butterfly diagrams for =
Figure 7 shows the longitudinal averages of the surface magnetic field for models, as afunction of time. One can directly see the effect of increased tilt angles (i.e. latitudinal separation of the preceding and the follower polarities) on the quicker formation of a stronger polar field. For the more rapidly rotating stars, the increase in both the flux emergence rate and the tilt angles cause the time at which the polar field reaches its maximum value to converge towards the time of the maximum total surface flux.
Figure 8 shows the unsigned radial field at the surface, averaged over longitude, for the cases shown in Fig. 7. In this way, it is easier to compare the distribution of the activity belts with that of polar fields, because the arithmetic cancellation of oppositepolarity fields is avoided during longitudinal averaging. It is evident that midlatitude activity strengthens considerably for and 8. Above , the polar fields start to become comparable with the midlatitude activity, and they also become dominant for an increasingly larger fraction of the cycle.
Fig. 7 Timelatitude diagrams of azimuthally averaged surface magnetic field for (panel a), (panel b), (panel c), and (panel d). For all cases, . We note thedifferent saturation levels of the colour scale. 

Open with DEXTER 
Fig. 8 Timelatitude diagrams of azimuthally averaged B_{r} for (panel a), (panel b), (panel c), and (panel d). for all cases.Here, the saturation levels are all 100 G. 

Open with DEXTER 
3.3 The distribution and coverage of starspots ( = )
To define spot areas from surface magnetic maps, we followed the simplest approach of setting a threshold to the magnetic field strength. All pixels with a field strength above the threshold are considered to belong to spots. We determine the threshold value from the condition that the time average of the spot coverage through the cycle, , roughly corresponds to the cycleaveraged umbral spot coverage observed on the Sun, which is about 0.002. Using this criterion, we have found a threshold of 187 G, corresponding to about 50% of B_{max} defined in Sect. 2.3. We have then used this threshold to determine the latitudinal distribution of the spot coverage for all four rotation rates with .
Figures 9a–d show timeaveraged latitudinal profiles of spot occupancy as a fraction of the stellar surface area. The area fractions were calculated by counting the pixels above the threshold, taking into account their areas on a spherical grid. In most observational studies, a factor of cos λ is used when estimating the fractional spot area per latitude bin. This means that our profiles are comparable with such latitudinal profiles presented in the literature (e.g. Järvinen et al. 2007; Waite et al. 2017). When the spot occupancy would be given as a fraction of the latitudinal band area, however, the polar spot of the case would lead to a much larger coverage near the pole than at midlatitudes. The spot coverage over the whole stellar surface, averaged over the whole activity cycle (darker curves), , increases from the solar value of 0.2–0.4% for , 1.2% for , and 10% for . For comparison, oneyear averages centred at the activity maximum are also given (lighter colours). There is a marked tendency for the mean latitude and the latitudinal spread of starspots to increase with and . The mean latitude at each hemisphere increases, owing to the Coriolis deflection of rising flux and the enhanced poleward transport of highly tilted source regions. The increase in the latitudinal spread is led by (i) higher flux emergence rate, (ii) the scaling of mean latitude with (Eq. (1)), and (iii) weaker flux cancellation between opposite polarities of BMRs, owing to a higher average tilt angle. It is evident from Fig. 9d () that, through the flux transport at sufficiently high and , spotted regions are formed at even higher latitudes than their emergence latitudes. Such starspots are formed by signed magnetic flux being transported towards higher latitudes and concentrated there. Figure 9e shows the variation of the surface coverage of starspots for all the four cases above. The cycle means are about onethird of the annual means around maxima, except for the case (1,1), for which the ratio is about 0.5.
Figure 10 shows snapshot surface distributions of signed field from the maximum phases of activity, corresponding to each of the considered rotation rates, for the case . The polar fields become stronger for higher , while the activity belts move towards higher latitudes. A strong unipolar polar cap forms for , but it is surrounded by patches of oppositepolarity field, owing to the leading polarities of freshly emerged, tilted BMRs. A conspicuous feature is visible for : d B∕dλ changes its sign at λ ≃ 75°, where the poleward meridional flow piles up the field diffusing from lower latitudes. In the later stages of the cycle, this ringlike structure diffuses to form a circular region peaking at the rotational pole.
In Fig. 11 we show the unsigned magnetic field strength, filtered with the threshold of 187 G. These snapshots can be seen as approximations of intensity images by omitting the limbdarkening and facular brightening effects. The strengthening of the polar field and of midlatitude activity with increasing is visible. Based on the spot threshold field strength, the polar spot forms only for , and it shows the same ringshaped pattern as in Fig. 10 with occasional plumes in the direction of the equator. Here, the solar case exhibits too large “sunspots”, though the total area fraction is solarlike. This is because the resolution of the SFT model is nothigh enough to resolve typical sunspots. In addition, the saturation level for the unsigned field is only 800 G, which is well below the average field strength of sunspot umbrae. Consequently, the simulated stellar images presented here are meant to represent a mediumresolution picture of spots on Sunlike stars. The resolution is therefore between those of solar fulldisc whitelight images and those produced by Doppler imaging (e.g. Solanki & Unruh 2004). A more rigourous treatment of spot areas is beyond the scope of this study, and will be employed in the next paper for proper modelling of brightness variations.
It is known that faster rotating G stars generally show stronger brightness variability (McQuillan et al. 2014). In this context, nesting of active regions can have an influence on the variability amplitudes on timescales comparable with the rotation period. This is because the surface distribution of spots would become less homogeneous in longitude (rotational phase), when spots tend to emerge within nests. As a visual demonstration of the effect of nesting, we show in Fig. 12 poleon views (i = 0) of our case at three different phases of the activity cycle, for random (unnested) and strongly nested cases, where the nesting probability was chosen to be p = 0.7 (same as in all the previous figures). We also display the signed field distributions corresponding to the nested case, for the overall field geometry and strength to be evaluated. Without forced nesting, the spot distribution appears more axisymmetric. With nesting, the highly nonaxisymmetric spot distribution is likely to induce largeramplitude brightness variations on rotational timescales.
Fig. 9 Timeaveraged latitudinal distributions of the fraction of surface area covered by starspots for (panel a), (panel b), (panel c), and (panel d), where , given at the top of each panel. The colours correspond to the cases in Figs. 3 and 6; lighter ones are averages over a oneyear window centred at the activity maximum and darker curves represent cycle averages. Timeaveraged surface coverages for maxima and wholecycles are given inside each plot. Panel e: variation of the total spot coverage for each case. 

Open with DEXTER 
Fig. 10 Snapshots of magnetic field strength from the runs for (panels a–d), corresponding to the cycle maximum in each case. The inclination angle of the rotation axis with respect to the line of sight is 30°. The latitudinal circles are drawn at λ = 37.5°, where the poleward flow speed has a maximum, and at λ = 75°, above whichit is assumed to be zero. 

Open with DEXTER 
Fig. 11 As in Fig. 10, but for the unsigned magnetic field strength with an inclination of 60°. The maps show the distribution of starspots, which are defined as all pixels above the threshold of 187 G. 

Open with DEXTER 
Fig. 12 Poleon views of the radial field, representing spot distributions for unnested (left column) and nested (middle column) cases for . The right column shows the corresponding signed magnetic field strength for the nested case, with a colour saturation at ± 870 G. Dotted circles represent the latitudes at 30° and 60°. We note that the spots are defined above 187 G. 

Open with DEXTER 
4 Discussion
We have developed a twopart model, to provide timeresolved maps of the radial magnetic field on Sunlike stars with rotation rates in the range Ω_{⊙}≤Ω ≤ 8Ω_{⊙}, which corresponds to an (sidereal) equatorial rotation period range from 25 to 3 days. The platform developed here will be used in the forward modelling of brightness variations on timescales covering active region evolution, stellar rotation, and the activity cycle (the second paper in this series). It also has the potential to be used in synthesising spectra covering photospheric lines used in Doppler imaging.
The thin fluxtube simulations successfully model the basic dynamical aspects of the emergence of largescale flux loops in the case of the Sun, despite several idealisations involved (Caligari et al. 1995). Distributions of tilt angles in faster rotating suns have not been well investigated so far, despite their potential effects on the distribution and evolution of surface magnetic flux. We have demonstrated in this study that these effects can be significant. Though we focused on the photospheric distribution of largescale radial fields, the dynamical effects of increasing rotation rate on the emerging flux would certainly have implications for coronal magnetic structure (Gibb et al. 2016).
The model provides clues about how patterns of stellar activity are likely to change with increasing rotation and fluxemergence rates. We assumed that the timelatitude pattern of flux eruptions at the base follows the solar butterflydiagram trends. Because currently there is no empirical evidence favouring any specific pattern for the internal toroidal field for , we preferred this simple and conservative approach to also test the applicability of the solar paradigm. Therefore, the only change in the dynamo with increasing rotation rate is that it produces more toroidal flux, which reaches higher latitudes at the base of the convection zone.
We extrapolated the empirical relation for the mean latitude of the solar butterfly pattern to higher levels of activity, as if a solar cycle had an amplitude times its Cycle22 level. This was done to obtain the base distributions of such very strong cycles before they rise to the surface, for 1Ω_{⊙} (step II in Appendix B). This step is necessary because observations show that sunspots appear at higher average latitudes during stronger cycles (Solanki et al. 2008; Jiang et al. 2011). However, it is also possible that in real stars the deviation from a solarlike butterfly diagram at the base of the convection zone can be significant, especially for more rapid rotators. In future work, extrapolation of composite or BabcockLeighton type dynamo models to faster rotating Sunlike stars can be employed (Işık et al. 2011; Karak et al. 2014) to obtain an estimate of the butterfly diagram of the internal toroidal field more consistently (see also Warnecke 2018).
Our simulations show that the polar field exhibits the following trends with increasing : (i) it strengthens relative to lower latitude activity, and (ii) it reverses its polarity increasingly earlier with respect to the activity maximum (Fig. 6d), in spite of the fact that the initial polar field was scaled with the rotation rate. The first trend results from two effects, namely, higher tilt angles and stronger activity. The second tendency can lead to an earlier amplification of the toroidal field by the action of differential rotation on the poloidal field. In addition to meridional circulation and dynamo effects (Jouve et al. 2010; Işık et al. 2011), it can thus contribute to shortening the cycle period for more active stars.
There are two competing effects in our models which determine the polar field amplitude: the increasing tilt angles (with rotation rate) tend to amplify the developing polar fields, whereas the gap of inactivity opening around the equator (with faster rotation) leads to weaker crossequatorial precedingpolarity flux cancellation, which limits the growth of a strong polar cap. We speculate that there is a critical rotation rate at which the two following timescales become comparable: the timescale for the magnetic flux of each polarity within a BMR to diffuse and cancel each other, and the timescale for the precedingpolarity flux from a typical BMR to diffuse and cancel with the precedingpolarity flux from a corresponding BMR from the opposite hemisphere. Beyond such a critical rotation rate, polar fluxes would be saturated, provided that we are using the same solar transport parameters. Such high rotation rates and activity levels are beyond the scope of this paper.
Another effect that can hinder the formation of circumpolar spots is related to the dynamo process. If the cycle period decreases and cycle overlap increases significantly with the rotation rate, then the polar fields resulting from SFT may not reach sufficient strengths to form spots before the subsequent cycle peaks (see the Sunlike model with P_{rot} = 2 days in Işık et al. 2011). However, circumpolar spots are indeed observed on some Doppler images of rapidly rotating Sunlike stars. Future observations and modelling of cycles on Sunlike stars should address this issue.
Our case can be compared to the Sunlike (G1.5V) star EK Draconis, which has . This star also exhibits a polar spot and midlatitude activity in several Doppler and ZeemanDoppler images (Järvinen et al. 2007; Waite et al. 2017). This qualitative agreement is gratifying. Furthermore, its mean umbral spot coverage is estimated to be in the range 0.25–0.40 based on TiOband observations by O’Neal et al. (2004), whereas our model for (8,8) gives a cycle mean of 0.07 and a maximum of 0.20. We consider these values to be in reasonable agreement given the rather simple thresholding used to determine spot areas, and also that we scale the amount of magnetic flux linearly with the rotation rate.
We applied nesting of BMR emergence as observed on the Sun to more active stars. The resulting SFT simulations led to substantial rotational asymmetry in the starspot distributions on active stars. If the observed starspots are in fact lowresolution manifestations of such active nests (Özavcı et al. 2018) then the observed sizes and lifetimes of these structures may not be indicative of the intrinsic sizes and lifetimes of their constituent spots (see Işık et al. 2007b, for solutions of monolithic vs. clustered unipolar spot diffusion).
ZeemanDoppler imaging studies show that the azimuthal component of the magnetic field strengthens and becomes comparable to or even exceeds the radial field for very active Sunlike stars. Our SFT model considers only the radial component of the field. Inclusion of the horizontal components (e.g. Gibb et al. 2016; Lehmann et al. 2017) represents a considerable extension of our current understanding, and they are beyond the present scope of modelling brightness variations of moderately active Sunlike stars. These components will therefore be introduced in a future study.
Our assumption that the meridional flow and differential rotation profiles are identical to those on the Sun has an effect on the resulting latitudinal distribution of the magnetic field. For instance, we expect that a meridional flow reaching up to latitudes higher than our assumed 75° would lead to a field peaking at the poles earlier in a given cycle. However, such details should not dramatically change the production of a polar spot. Deviations of these profiles from the solar patterns would affect the evolution of flux in various latitudinal zones. Alternative profiles from theoretical models of flow fields as a function of the rotation rate can be used in our model (Küker et al. 2011; Kitchatinov & Olemskoy 2012).
Finally, we note that it might be possible to reproduce the observed surface patterns of activity starting from different combinations of dynamogenerated toroidal field patterns, initial fluxtube conditions, and largescale flows in the convective envelope. However, our aim has not been to seek the parameters and functions which give the best fits to the observed patterns, but to stick to the solar paradigm and test its validity for more active and more rapidly rotating configurations. In qualitative terms, our models are more or less consistent with the overall observed patterns of surface inhomogeneities in longitude and latitude, as well as surface fractions. We reserve quantitativecomparisons with the observed stellar brightness variations for the followup studies, in which we shall estimate fractions of spot and facular areas from the modelled surface magnetic fields, and generate light curves.
We planto include further features into the model in forthcoming studies, such as inflows towards the activity belts, active longitudes, and random scatter in the tilt angle, which will be relevant for multiplecycle models. Such parametrisations can also be implemented as part of an extrapolation of a solar fluxtransport dynamo model.
5 Conclusion
Under the assumption that the toroidal field at the base of the convection zone follows the extrapolated solar patterns, we have developed a numerical platform combining two main physical effects which are likely responsible for the observed variety of activity patterns on Sunlike stars: (i) the rise of flux tubes under 3D effects of the relevant largescale forces, and (ii) the surface transport of emerging largescale magnetic fields, under the effects of largescale surface flows. We find that for the final modelled distribution of magnetic fields on stellar photospheres, the following proposed processes play important roles: (i) the deviation from radial fluxtube rise (Schüssler & Solanki 1992), (ii) the evolution of the field on the stellar surface (Schrijver & Title 2001), and (iii) increasing BMR tilt angles (Işık et al. 2007b, 2011). We also find that the onset of polar spot formation on Sunlike stars can occur by accretion of trailing polarity flux from BMRs, between four and eight times the solar rotation rate and the flux emergence rate. Our results also show that nesting of emerging bipoles can have substantial effects on the surface distribution of starspots.
The models developed here can be used in forward modelling of brightness variations in magnetically active Sunlike stars. At a later stage, we plan to extend them to be helpful in understanding other observations of active Sunlike stars.
Acknowledgements
We thank Jie Jiang and Robert Cameron for sharing their code which generates semisynthetic solar emergence records, and Maria Weber for confirming our fluxtube emergence results using a different code. We also thank the referee for the comments that helped to improve this manuscript. E.I. acknowledges support by the Young Scientist Award Programme BAGEP2016 of the Science Academy, Turkey. This work has been partially supported by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea. A.S. acknowledges funding from the European Research Council under the European Union Horizon 2020 research and innovation programme (grant agreement No. 715947).
Appendix A Synthetic input solar cycle
A.1 Emergence frequency
To model the temporal profile of the monthly group sunspot number, R_{G}, we follow CJS16 and use the function devised by Hathaway et al. (1994) with additional Gaussian random noise, Δ R_{G}: (A.1)
where t is the time in months and a is the amplitude. We fit the first function on the RHS of Eq. (A.1) to sunspot group data for cycle 22 from RGO, which yields a = 0.00336. The quantities b(a) and c control the length and the asymmetry of the cycle, respectively. For b, we adopt Eq. (4) of Hathaway et al. (1994) and c = 0.71. The standard deviation of the probability distribution of ΔR_{G}(t) around R_{G}(t) is approximated by 0.5R_{G}(t). This level was estimated by measuring the deviations of the observed R_{G} of Cycles 21–23 from the corresponding fits of the first term on the RHS of Eq. (A.1). Finally, the number of BMRs emerging in a given month is obtained as R_{G}∕2.75, based on a calibration carried out by CJS16, through fits to the monthly number of groups observed for Cycles 21–24.
A.2 Sunspot group latitudes
Next, we set up a synthetic butterfly diagram in the input model, following the procedure of JCSS11 (also followed by CJS16). We take the average latitude of spot groups at a given temporal phase bin i of the cycle to be described by the quadratic function (A.2)
where 1 ≤ i ≤ 30 (the cycle is split into 30 temporal bins), is the average latitude oversolar Cycles 12–20, and ⟨λ⟩ is the average latitude of sunspot groups over the cycle. This was obtained by JCSS11 for a given cycle using a linear fit to data from all the available cycles, as (A.3)
where S_{⊙} is the cycle amplitude. To representthe Sun, we set S_{⊙} = 156 from the maximum of the twelvemonth running mean of the observed R_{G} of Cycle 22 (see Table 1 of JCSS11). The second term on the RHS of Eq. (A.3) models the observation that stronger solar cycles have a higher mean latitude of emergence (Solanki et al. 2008). For the latitudinal spread around ⟨λ ⟩, the model assumes a Gaussian distribution with a standard deviation σ^{i} that varies with the cycle phase, following the quadratic function (A.4)
where 1 ≤ i ≤ 30 as in Eq. (A.2).
A.3 Sunspot group areas
Sunspot group areas (A) are randomly picked from either a power law or a lognormal distribution, depending on the size (in millionths of a solar hemisphere), μH, as given by (A.5)
To include the dependence of the mean group area on the cycle phase, we adopt the relation given by JCSS11, but with 1 ≤ i ≤ 20, (A.6)
Appendix B The link between the base and the surface
To synthesise starspot emergence records, we follow the algorithm outlined below. The steps I–III correspond to those in Fig. B.1.
 I.
Generate a random realisation of the spot group record (SGR) with a certain flux emergence rate (see Sects. 2.1.1–2.1.2). The mean latitude of this SGR depends on (Eq. (1)).
 II.
Map the latitudes obtained in step I down to the base of the convection zone, by interpolating the fluxtube rise table for the solar rotation rate, (Sect. 2.2.2).
 III.
Take the time series resulting from step II as the base latitudes of the flux tubes leading to spot groups, for a given rotation rate, . Generate the emergence latitudes and tilt angles for , by interpolating within the fluxtube rise table (Sect. 2.2.2).
 IV.
To simulate the effect of active region nesting, modify the longitudes and latitudes obtained in the previous steps, using a probabilistic approach (Appendix C).
Fig. B.1 A schematic representation of the mappings between the surface and the base of the convection zone in the model. Ω denotes the equatorial rotation rate in solar units. The roman numerals refer to the steps of the procedure. 

Open with DEXTER 
Appendix C Nests of activity
To simulate the observed tendency of flux emergence in the vicinity of recent emergence, we modified the longitudes and latitudes of BMRs in our standard starspot record, using a probabilistic approach. We first set a generic probability 0 < p < 1 for each BMR with coordinates (λ, ϕ) to be part of a nest. In this study, we chose a rather high probability of p = 0.7 to clearly demonstrate the effects of nesting. This value has made the cycle variation of the equatorial dipole moment much more similar to the observed variations (Wang 2014, see Fig. 3) in comparison to the unnested case.
If a uniformly chosen random number E_{i} ~ U[0, 1] is less than p, then the coordinates of the ith BMR of the unnested record is considered as a potential nest centre, . If the next random number E_{i+1} < p, then the (i + 1)th BMR belongs to the nest of the ith BMR. The coordinates of such a BMR are then replaced by the coordinates , which are drawn from normal distributions centred around with widths 2° in latitude and 3° in longitude, close to the empirical values obtained by Pojoga & Cudnik (2002) for Cycle 23. The procedure is continued iteratively until for the (i + k)th BMR E_{i+k} ≥ p holds. Such a BMR is considered as an isolated spot group, while its original coordinates are kept unchanged.
References
 Aigrain, S., Favata, F., & Gilmore, G. 2004, A&A, 414, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balona, L. A.,& Abedigamba, O. P. 2016, MNRAS, 461, 497 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J.F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Baumann, I., Schmitt, D., Schüssler, M., & Solanki, S. K. 2004, A&A, 426, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baumann, I., Schmitt, D., & Schüssler, M. 2006, A&A, 446, 307 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdyugina, S. V. 2005, Liv. Rev. Sol. Phys., 2 [Google Scholar]
 Caligari, P., MorenoInsertis, F., & Schussler, M. 1995, ApJ, 441, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Caligari, P., Schüssler, M., & MorenoInsertis, F. 1998, ApJ, 502, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R. H., Jiang, J., Schmitt, D., & Schüssler, M. 2010, ApJ, 719, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R. H., Jiang, J., & Schüssler, M. 2016, ApJ, 823, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Castenmiller, M. J. M., Zwaan, C., & van der Zalm, E. B. J. 1986, Sol. Phys., 105, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., & Landstreet, J. D. 2009, ARA&A, 47, 333 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Elstner, D., & Korhonen, H. 2005, Astron. Nachr., 326, 278 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., Fisher, G. H., & Deluca, E. E. 1993, ApJ, 405, 390 [NASA ADS] [CrossRef] [Google Scholar]
 FerrizMas, A., & Schüssler, M. 1995, Geophys. Astrophys. Fluid Dyn., 81, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Gibb, G. P. S., Mackay, D. H., Jardine, M. M., & Yeates, A. R. 2016, MNRAS, 456, 3624 [NASA ADS] [CrossRef] [Google Scholar]
 Granzer, T., Schüssler, M., Caligari, P., & Strassmeier, K. G. 2000, A&A, 355, 1087 [NASA ADS] [Google Scholar]
 Hathaway, D. H., Wilson, R. M., & Reichmann, E. J. 1994, Sol. Phys., 151, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Holzwarth, V., Mackay, D. H., & Jardine, M. 2006, MNRAS, 369, 1703 [NASA ADS] [CrossRef] [Google Scholar]
 Işık, E., Schmitt, D., & Schüssler, M. 2007a, Astron. Nachr., 328, 1111 [NASA ADS] [CrossRef] [Google Scholar]
 Işık, E., Schüssler, M., & Solanki, S. K. 2007b, A&A, 464, 1049 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Işık, E., Schmitt, D., & Schüssler, M. 2011, A&A, 528, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Järvinen, S. P., Berdyugina, S. V., Korhonen, H., Ilyin, I., & Tuominen, I. 2007, A&A, 472, 887 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jeffers, S. V., Barnes, J. R., & Collier Cameron, A. 2002, MNRAS, 331, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Işık, E., Cameron, R. H., Schmitt, D., & Schüssler, M. 2010, ApJ, 717, 597 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, J., Cameron, R. H., Schmitt, D., & Schüssler, M. 2011, A&A, 528, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jiang, J., Hathaway, D. H., Cameron, R. H., et al. 2014, Space Sci. Rev., 186, 491 [NASA ADS] [CrossRef] [Google Scholar]
 Jouve, L., Brown, B. P., & Brun, A. S. 2010, A&A, 509, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karak, B. B., Kitchatinov, L. L., & Choudhuri, A. R. 2014, ApJ, 791, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2012, MNRAS, 423, 3344 [NASA ADS] [CrossRef] [Google Scholar]
 Koch, D. G., Borucki, W. J., Basri, G., et al. 2010, ApJ, 713, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Küker, M., Rüdiger, G., & Kitchatinov, L. L. 2011, A&A, 530, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehmann, L. T., Jardine, M. M., Vidotto, A. A., et al. 2017, MNRAS, 466, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, S. C., Donati, J.F., Semel, M., Petit, P., & Carter, B. D. 2006, MNRAS, 370, 468 [NASA ADS] [CrossRef] [Google Scholar]
 McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 MorenoInsertis, F. 1986, A&A, 166, 291 [NASA ADS] [Google Scholar]
 O’Neal, D., Neff, J. E., Saar, S. H., & Cuntz, M. 2004, AJ, 128, 1802 [NASA ADS] [CrossRef] [Google Scholar]
 Özavcı, I., Şenavcı, H. V., Işık, E., et al. 2018, MNRAS, 474, 5534 [NASA ADS] [CrossRef] [Google Scholar]
 Pojoga, S., & Cudnik, B. 2002, Sol. Phys., 208, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Reiners, A. 2012, Liv. Rev. Sol. Phys., 9, 1 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telescopes Instrum. Syst., 1, 014003 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J. 2001, ApJ, 547, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Title, A. M. 2001, ApJ, 551, 1099 [NASA ADS] [CrossRef] [Google Scholar]
 Schüssler, M., & Solanki, S. K. 1992, A&A, 264, L13 [NASA ADS] [Google Scholar]
 Schüssler, M., Caligari, P., FerrizMas, A., Solanki, S. K., & Stix, M. 1996, A&A, 314, 503 [NASA ADS] [Google Scholar]
 Shapiro, A. I., Solanki, S. K., Krivova, N. A., et al. 2014, A&A, 569, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shapiro, A. I., Solanki, S. K., Krivova, N. A., et al. 2017, Nat. Astron., 1, 612 [NASA ADS] [CrossRef] [Google Scholar]
 Skaley, D., & Stix, M. 1991, A&A, 241, 227 [NASA ADS] [Google Scholar]
 Solanki, S. K., & Unruh, Y. C. 2004, MNRAS, 348, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Solanki, S. K., Wenzler, T., & Schmitt, D. 2008, A&A, 483, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Spruit, H. C. 1981, A&A, 98, 155 [NASA ADS] [Google Scholar]
 Strassmeier, K. G. 2009, A&ARv, 17, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Waite, I. A., Marsden, S. C., Carter, B. D., et al. 2015, MNRAS, 449, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Waite, I. A., Marsden, S. C., Carter, B. D., et al. 2017, MNRAS, 465, 2076 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M. 2014, Space Sci. Rev., 186, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J. 2018, A&A, 616, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
We note that the fluxes of individual flux tubes with the chosen B_{0} and λ_{0} (Sect. 2.2) are in the same range with the BMRs used in the SFT (Sect. 2.3). Whenever a surface source is introduced in the SFT simulation, however, its flux is determined only by the empirically synthesised size distribution.
All Figures
Fig. 1 Semisynthetic emergence record (SGR) for solar cycle 22, using emergence latitudes and tilt angles resulting from the remapping procedure for the solar rotation and flux emergence rate . The colour bar shows the tilt angles with a rather high saturation level, so that the tilt angles can be compared with faster rotating cases. The plot on the left shows the histogram of emergence latitudes. 

Open with DEXTER  
In the text 
Fig. 2 Initial latitudes and field strengths of flux tubes (red diamonds) chosen for fluxtube rise computations, plotted over linear stability diagrams for flux tubes at the middle of the overshoot region, for (panel a), (panel b), (panel c), and (panel d), with differential rotation ΔΩ_{⋆} = ΔΩ_{⊙}. The contour lines denote the growth time in days. The linearly stable regime is shown in white, and the unstable regime is shaded with light grey where the fastestgrowing mode is m = 1, and with darkgrey for m = 2. The red diamonds show the initial tube parameters chosen for the numerical simulations, at a growth time of 50 days. 

Open with DEXTER  
In the text 
Fig. 3 Panel a: tilt angle of emerging flux loops vs. their emergence latitude for (black diamonds), (blue asterisks), (green squares), and (red crosses). Panel b: latitudinal deflections (poleward is positive) as a function of the initial latitude at the base of the convection zone. The initial λ_{0} and B_{0} are as inFig. 2. Emergence latitudes for which no activeregionsized BMR emerges for are shown by small grey diamonds (see text). 

Open with DEXTER  
In the text 
Fig. 4 Geometry of three emerging flux loops with initial latitudes 45° (left panels),46° (middle panels), and 47° (right panels) at the base of the convection zone, for . Upper panels: inner sphere is drawn at 0.72R_{⊙}, and the outer sphere at 0.97R_{⊙}. The parts of the tube that are beneath the outer layer are shaded in grey, whereas the emerged parts are brighter. For each mass element of the tube, the colours denote the crosssectional tube diameter (the redder the larger). Lower panels: latitudinal and radial projections of the tubes, where each mass element is represented by a dot. The horizontal line on the radial profile corresponds to the location of the outer sphere (0.97R_{⊙}), where the tilt angle (α) is measured from the footpoint locations. The red arrows denote the apex of each tube. 

Open with DEXTER  
In the text 
Fig. 5 As in Fig. 1, but now for various sets of as indicated at the top of each panel. Left panels: for (solar emergence rate). Right panels: for . 

Open with DEXTER  
In the text 
Fig. 6 Total unsigned surface flux for (panel a) and (panel b), and the average polar flux density for latitudes poleward of ± 75°, for (panel c) and (panel d). The colours denote , as in Fig. 3. Panel e: comparison of the time integrals of polar flux based on the average polar field and the total unsigned flux, for (diamonds) and (asterisks). 

Open with DEXTER  
In the text 
Fig. 7 Timelatitude diagrams of azimuthally averaged surface magnetic field for (panel a), (panel b), (panel c), and (panel d). For all cases, . We note thedifferent saturation levels of the colour scale. 

Open with DEXTER  
In the text 
Fig. 8 Timelatitude diagrams of azimuthally averaged B_{r} for (panel a), (panel b), (panel c), and (panel d). for all cases.Here, the saturation levels are all 100 G. 

Open with DEXTER  
In the text 
Fig. 9 Timeaveraged latitudinal distributions of the fraction of surface area covered by starspots for (panel a), (panel b), (panel c), and (panel d), where , given at the top of each panel. The colours correspond to the cases in Figs. 3 and 6; lighter ones are averages over a oneyear window centred at the activity maximum and darker curves represent cycle averages. Timeaveraged surface coverages for maxima and wholecycles are given inside each plot. Panel e: variation of the total spot coverage for each case. 

Open with DEXTER  
In the text 
Fig. 10 Snapshots of magnetic field strength from the runs for (panels a–d), corresponding to the cycle maximum in each case. The inclination angle of the rotation axis with respect to the line of sight is 30°. The latitudinal circles are drawn at λ = 37.5°, where the poleward flow speed has a maximum, and at λ = 75°, above whichit is assumed to be zero. 

Open with DEXTER  
In the text 
Fig. 11 As in Fig. 10, but for the unsigned magnetic field strength with an inclination of 60°. The maps show the distribution of starspots, which are defined as all pixels above the threshold of 187 G. 

Open with DEXTER  
In the text 
Fig. 12 Poleon views of the radial field, representing spot distributions for unnested (left column) and nested (middle column) cases for . The right column shows the corresponding signed magnetic field strength for the nested case, with a colour saturation at ± 870 G. Dotted circles represent the latitudes at 30° and 60°. We note that the spots are defined above 187 G. 

Open with DEXTER  
In the text 
Fig. B.1 A schematic representation of the mappings between the surface and the base of the convection zone in the model. Ω denotes the equatorial rotation rate in solar units. The roman numerals refer to the steps of the procedure. 

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