Testing the Solar Activity Paradigm in the Context of Exoplanet Transits

The flux-transport model used here was developed by Schrijver (2001), with cycle modulation introduced by Schrijver & Title (2001) and with flux-decay properties as discussed by Schrijver et al. (2002). This same flux-transport code continues to be used in an assimilation mode for Solar and Heliospheric Observatory/MDI and Solar Dynamics Observatory/HMI magnetogram data starting in mid-1996, see, for example, Schrijver & DeRosa (2003) and Schrijver & Liu (2008).^{1 ,2}

Visualizations of the global field patterns are discussed below, but other realizations associated with two of the model runs here (identified as ${ \mathcal M }({ \mathcal A }=1)$ and ${ \mathcal M }({ \mathcal A }=30);$ see below in this section for their properties) can be found in Schrijver et al. (2003).

At each time step, a set of bipolar regions is randomly selected from a parent-distribution power-law frequency spectrum that is based on the work by Harvey & Zwaan (1993): they used magnetogram sequences throughout solar Cycle 21 (1975–1986) to determine times of maximum development of bipolar regions and with that determined the behavior of the flux-input spectrum with cycle phase. Their findings support a power-law distribution of bipolar-region sizes for which the power-law index is essentially unchanged through the cycle, so that only a multiplier is needed to describe the evolution through the activity cycle. The same spectrum is used throughout the simulated activity cycles, but modulated in time with a profile resembling a solar sunspot cycle, with successive cycles somewhat overlapping in early and late cycle phases. Stars of different activity levels are simulated by a multiplier on the cycle amplitude. As discussed below, a maximum flux for active regions is set as desired, below which the distribution is maintained the same in shape, except for an overall multiplicative factor.

These regions are then randomly distributed over a range of latitudes and with a range of tilt angles. The latitudes and tilt angles have a spread that decreases with increasing size of the region, while the mean latitude shifts from mid-latitudes toward the equator as a cycle progresses, with a 3 yr phase of overlap between successive cycles, all based on the average values and scatter about those derived from solar observations. Random values are generated, drawn from parent-distribution functions that are based on observed properties of solar bipolar regions, ranging from small ephemeral regions to large active regions.

The longitude of flux emergence for each bipole is in principle randomly drawn from a uniform distribution. However, this is modulated by where previously emerged regions continue to exist: on the Sun, bipolar regions have a strong tendency to emerge at locations where a previous generation emerged and often still exists. This nesting property is so pronounced that almost one in two active regions emerges inside another (Brouwer & Zwaan 1990; Harvey & Zwaan 1993). The nesting as seen for solar regions is modeled as described by Schrijver (2001). In short, the probability of emergence per unit area inside a magnetic plage region is set to 22 times that outside of a plage region (following Harvey & Zwaan 1993). Other than this hysteresis in flux emergence, no preferred longitudes are introduced.

The selected fluxes are distributed as flux concentrations over two adjacent patches of opposite polarity. Work by Schrijver & Harvey (1994) revealed that mature active regions have a mean flux density of 100–150 G (excluding spots and pores), regardless of the region's age or size. Here, a mean flux density is chosen of 180 G (as in the original model by Schrijver 2001) to accommodate a characteristic fraction of the flux contained in spots and pores. The regions are gradually introduced, with an equivalent rate of flux emergence of 5×10²¹ Mx day⁻¹ (based on data from Harvey & Zwaan 1993).

The flux-transport model uses a point-source approximation for the magnetic field. These sources can be thought of as moving along a pattern of vertices that are defined by the evolving supergranular network. The properties of that cellular pattern and the characteristic displacement velocity along its vertices determine the typical mean-free path for the flux concentrations. At each time step, the flux-transport model assumes that all flux within an area with a radius equal to the typical mean-free path length for the quiet-Sun network coagulate into a single concentration prior to the next time step. This concept was successfully tested for a range of solar activity levels. The derivation of the coagulation radius r_c = 4200 km is described in Section 4 of Schrijver (2001), based on the model by Schrijver et al. (1997).

In reality, the supergranular and larger-scale flows that transport magnetic flux are, of course, not strictly confined to the lanes between supergranules, while moreover the supergranular cell pattern itself is evolving. Hence, the flux concentrations in the model are abstractions of patches within the photospheric field for which the diffusion approximation can be applied on scales beyond the coagulation radius. On scales below r_c, granular convective motions interact with the magnetic flux through the laws of radiative magnetoconvection, reaching a quasi-equilibrium on a timescale of a few granular turnover timescales, well below the step time of 0.25 day used in the model. Note that the diameter of the coagulation patch of 2r_c = 8400 km is comparable to the box width of 9000 km used for the numerical models by Norris et al. (2017), which is used in Section 3.3 to model the photospheric brightness. Hence, Section 3.3 uses the fluxes Φ_i in the point sources of the surface flux-dispersal model, converted into a mean flux density of ${{\rm{\Phi }}}_{i}/(\pi {r}_{{\rm{c}}}^{2})$.

The flux concentrations initiate a grid-free random walk in which step lengths are a function of the absolute flux of the concentrations, again following observed trends that show larger concentrations to be less mobile, i.e., apparently more able to resist convective displacement.

These concentrations can collide to merge or (partially) cancel, and they can fragment, with collision cross sections and fractionation probabilities depending on the flux contained, as derived from solar observations.

Then, before starting on the next step, all concentrations are moved subject to the large-scale differential-rotation and meridional-flow profiles, after which a new set of bipolar regions is selected to match the average flux-input rate, which itself is modulated following the progression of an activity cycle shaped as the average sunspot cycle.

I use the following parameter settings (for details, see Schrijver & Title 2001): A time step of 6 hr, a flux-dispersal coefficient of D = 300 km s⁻² with a flux-dependent step size as described by Equations (A4) and (A5) in Schrijver (2001). The differential-rotation rate is set to that from Komm et al. (1993) and the meridional flow, tapered at high latitudes, like that of Van Ballegooijen et al. (1998), both as applied in the study by Schrijver & Title (2001). The half-life timescale on which flux concentrations "decay" (i.e., are randomly removed from the simulation) as introduced and described by Schrijver et al. (2002) is set to 5 yr. The model is set to use a fixed-amplitude cycle with equal strengths at each sunspot maximum, set to reach its peak 4 yr after sunspot minimum, with a cycle-to-cycle overlap period of 3 yr (see Schrijver et al. 2002). The duration of the full magnetic cycle is set to 21.9 yr, and its associated flux-injection profile in time and latitude is as in Equation (2) in Schrijver & Title (2001).

The model is run for five different sets of only two free parameters. The first, ${ \mathcal A }$, quantifies the amplitude of the stellar cycle; it is simply a multiplier on the frequency at which bipolar regions are inserted compared to a typical solar cycle, so that ${ \mathcal A }=1$ is used for a run to mimic the Sun (calibrated to Cycle 21), while ${ \mathcal A }\gt 1$ signifies a more-active star with a more pronounced cycle. The second parameter, Φ_max, is a threshold flux in bipolar regions at which the power-law frequency spectrum of fluxes is truncated, as discussed below. The parameters and results of the five runs are summarized in Table 2.

Table 2.  Model Properties and Modulation Estimates

Run	Φmax	$\langle | {fB}| \rangle$,	fd	Net Intensity			Rotation			Transit Residual
	(1022	$\langle | {fB}| {\rangle }_{\mathrm{fac}}$	(%)	Change δ* (%)			Modulation (%)			Ampl. Δ(λ) (ppm)
	Mx)	(G)		3870	6010	15975	3870	6010	15975	3870	6010	15975
${ \mathcal M }({ \mathcal A }=1)$	1.5	16, 15	0.06	+0.24	+0.04	+0.014	0.06	0.014	0.010	180	40	40
${ \mathcal M }{({ \mathcal A }=1)}^{+}$	30	60, 51	0.65	+0.55	−0.01	−0.00	0.38	0.31	0.12	420	180	70
${ \mathcal M }{({ \mathcal A }=10)}^{+}$	30	212, 132	5.9	−0.97	−1.60	−0.40	0.90	0.64	0.25	980	520	300
${ \mathcal M }({ \mathcal A }=30)$	1.5	135, 100	2.6	+0.26	−0.54	−0.12	0.30	0.22	0.08	1100	540	190
${ \mathcal M }{({ \mathcal A }=30)}^{+}$	30	479, 175	22.	−6.1	−5.3	−1.3	2.9	2.0	0.7	6300	3800	1400

Note. The table lists the following for a single, arbitrarily chosen step in the flux-transport model for an activity phase about 1.5 yr after cycle maximum, following the run identifier, which includes the cycle-strength factor ${ \mathcal A }$, and the maximum active-region flux per polarity Φ_max: globally averaged absolute magnetic flux density $\langle | {fB}| \rangle$ of all flux elements and $\langle | {fB}| {\rangle }_{\mathrm{fac}}$ for the faculae (with fluxes below 3×10²⁰ Mx); globally averaged areal filling factor f_d of dark pores and spots assuming an intrinsic field strength of 1.35 kG; the contrast δ_* of the total unocculted brightness of the active star and the unocculted inactive star from the perspective shown in Figure 1; peak-to-trough rotational modulation, and—in the final set of three columns—peak-to-trough amplitude Δ(λ) for the transit residuals, both as shown in Figures 1(a)–(d). Note that δ_*(λ) and Δ(λ) are for rough intercomparison only as they are computed for a single time step in the flux-dispersal model and for a transit across a single viewing angle of the stellar sphere.

Download table as:  ASCIITypeset image

A star like the present-day Sun is simulated as in the model developed by Schrijver (2001), here referred to as ${ \mathcal M }({ \mathcal A }=1)$. For this baseline solar model, the bipole-emergence rate is calibrated to the solar cycle, i.e., with a cycle-strength multiplier as in Schrijver & Title (2001) of ${ \mathcal A }=1$.

For stars like Kepler-17 and WASP-36, of solar spectral type but with rotation periods around 12 days, a coronal soft X-ray flux density is expected that is about 10 times higher than that of the Sun (e.g., Patten & Simon 1996; Mittag et al. 2018). As the soft X-ray flux density scales roughly linearly with the surface magnetic flux density (e.g., Schrijver & Title 2005), the mean surface magnetic flux density, $\langle | {fB}| \rangle$ (for filling factor f and intrinsic photospheric field strength B), should be roughly an order of magnitude higher than for run ${ \mathcal M }({ \mathcal A }=1)$, here particularly for the facular fields because spot fields appear to contribute little to the coronal brightness, so roughly ∼150 G. Reaching such levels of surface magnetic activity requires significant enhancements in the active-region injection frequency or mean active-region size, or both.

This estimate of the average surface magnetic flux density $\langle | {fB}| \rangle \sim 150$ G for stars of about half the solar rotation period is subject to a considerable uncertainty range. One discussion of this (with associated references) can be found in Vidotto et al. (2014): using the five power-law scalings from various sources in the literature that are listed in their Table 3 based on Zeeman broadening to quantify the surface magnetic field, one infers $\langle | {fB}| \rangle =40\mbox{--}160$ G starting from a solar value from run ${ \mathcal M }({ \mathcal A }=1)$ of ∼15 G. A value of ∼170 G results from a calibration from Rossby number via Zeeman–Doppler imaging to an equivalent magnetic flux density for a Zeeman broadening signal using Figure 1 and Equation (2) in See et al. (2019). Thus, the value of ∼150 G estimated above seems reasonable, although a substantial uncertainty should be allowed for.

Zoom In Zoom Out Reset image size

A much more active star, but with otherwise solar parameters, is simulated using ${ \mathcal A }=30$ in model ${ \mathcal M }({ \mathcal A }=30)$. The total flux on the stellar surface in the latter run peaks at about 7 times that of ${ \mathcal M }({ \mathcal A }=1)$: inserting more active regions with a given size spectrum leads to more flux cancellations, so that the total flux on the surface increases more slowly than linearly with the cycle strength (see, e.g., Schrijver 2001 for a discussion). Whereas it reaches a value of $\langle | {fB}| \rangle$ commensurate with that of a star with ${P}_{\mathrm{rot}}\sim 12$ days, the rotational amplitude is too low (see Section 4).

To increase the rotational modulation, larger regions could be introduced. This also strengthens the impact of activity nests, the tendency for active regions to emerge within existing other regions or in sites of an earlier such emergence. Thus, another run, ${ \mathcal M }{({ \mathcal A }=30)}^{+}$, also with ${ \mathcal A }=30$, is created but with a maximum active-region flux per polarity of Φ_max raised from 1.5×10²² Mx (as used in Schrijver 2001) to 3×10²³ Mx. The latter value allows regions as large as that expected to power the Carrington–Hodgson flare to emerge (based on a flux estimate by Aulanier et al. 2013). Such very large regions are rare on the Sun (the largest region observed to date, which occurred in 1947 April, is estimated to have had a flux at that level, see Taylor 1989; Schrijver et al. 2012), hence the choice of ${{\rm{\Phi }}}_{\max }=1.5\times {10}^{22}$ Mx for the solar activity level in run ${ \mathcal M }({ \mathcal A }=1)$, but the occurrence of very large flares on more-active stars suggests that larger regions may be more common on such stars (see, e.g., Maehara et al. 2012; Shibayama et al. 2013; Notsu et al. 2017).

The model runs described above suggested that stars like Kepler-17 and WASP-36 may be approximated by a run of intermediate ${ \mathcal A }$ and Φ_max, with parameters set as in run ${ \mathcal M }{({ \mathcal A }=10)}^{+}$. For comparison purposes, also a run ${ \mathcal M }{({ \mathcal A }=1)}^{+}$ is executed for a star of solar activity, but allowing for large active regions, setting Φ_max = 3 × 10²³ Mx.

Note that all runs for an active Sun-like star with ${ \mathcal A }=30$ have the same total amount of flux emerging per unit time in the range from $6\times {10}^{18}$ to 1.5×10²² Mx per polarity. The larger range of the power-law probability distribution for regions to be inserted into the stellar photosphere in runs with elevated Φ_max introduces additional large regions, leading to more flux on the stellar surface, while maintaining the frequencies of the more abundant smaller regions to match the shape of the size spectrum derived from solar observations (while allowing an overall multiplicative factor). Given the properties of the active-region frequency distribution, the rate of flux input for Φ_max = 3×10²³ Mx is increased by a factor of 1.7 relative to Φ_max = 1.5×10²² Mx. Larger regions disperse more slowly, which also leads to increased persistence of successive generations of flux emergence within them (in active-region nests), so that the total amount of flux in the photosphere increases more than the rate at which flux input increases if the spectrum is extended to higher fluxes (see Table 2).

To bootstrap the simulation runs relatively quickly into a state in which the polar caps are less sensitive to the initial state (an empty stellar surface), for all runs except ${ \mathcal M }{({ \mathcal A }=30)}^{+}$, half of all flux concentrations is randomly removed after the completion of the first 11 yr sunspot cycle (using the same procedure as in Schrijver & Title 2001). For these runs, the data shown here are taken from the third sunspot cycle, i.e., covering the years 22 to 33 for the simulations. Any effects of the initial state or the removal of half of all concentrations after the first 11 yr are moreover dampened by the half-life decay timescale of flux concentrations of 5 yr. For run ${ \mathcal M }{({ \mathcal A }=30)}^{+}$, when the code is tracking over 2×10⁵ elements near cycle maximum, the n² nature of the algorithm testing for collisions between flux concentrations leads to much longer computation time; because only the high-latitude fields require a multi-decade run, while these have little impact on the transit residuals being studied, for this run, data are shown in the same cycle phases as for the others, but from the first rather than third sunspot cycle.

Sample simulation results for the magnetic field distributions across the stellar surface are shown in Figures 1(a)–(d). As the focus of this study is the signal from exoplanetary transits, and as most planetary systems appear to have their normal vector rather well aligned with the stellar rotation axis, the perspective for visualizations of these simulations is chosen to lie in the plane of the stellar rotational equator.

Estimation of the brightness of facular patches on the faces of stars requires knowledge of their contrast as a function of viewing angle (or relative distance from disk center), wavelength, and corresponding characteristic absolute magnetic flux densities. Here, disk appearance, rotational modulation, and transit signals are simulated for three wavelengths for which Norris et al. (2017) show facular contrast versus limb distance (their Figure 8) and activity level: 3870, 6010, and 15975 Å. Their values are based on magnetoconvective calculations with radiative transfer using the MURaM code, validated successfully against solar observations³ . For values of μ between 0 and 0.2, close to the limb, for which they show no results, a smooth transition to zero contrast at the limb is imposed.

The background limb darkening of the quiet photosphere is set to that measured by Neckel & Labs (1994), except for the 15975 Å channel, for which I use the polynomial fit to the data presented in Figure 4 by Norris et al. (2017).

To estimate the average magnetic flux density to be used to compute facular brightenings, I assume that the flux in a model concentration is spread out over an area with a radius equal to the coagulation radius of r_c = 4200 km: concentrations that approach each other closer than that in the model are considered to merge. The flux threshold at which bright faculae transition to dark pores or spots is set to Φ_b = 3×10²⁰ Mx (see Table 4.1 in Schrijver & Zwaan 2000).

For spots, the average intrinsic field strength over the combined area of umbra and penumbra is set to 1.35 kG based on the study by Solanki & Schmidt (1993), who note that this value is similar to the field strength characteristic of smaller magnetic features; this value is therefore used as intrinsic field strength in the visualizations here for all flux concentrations. For the intensity contrast averaged over an entire pore or spot, I assume an average intensity contrast derived from the ratio of blackbody temperatures of the spot and photosphere, T_spot and T_phot, respectively, corresponding to the wavelength λ considered: ${c}_{{\rm{s}}}(\lambda )\,=\left(\exp ({hck}/(\lambda {T}_{\mathrm{phot}}))-1\right)/\left(\exp ({hck}/(\lambda {T}_{\mathrm{spot}}))-1\right)$, for h the Planck constant, c the speed of light in vacuum, and k the Boltzmann constant. Radiative transfer effects for different wavelengths are ignored in the present approximation.

For the Sun-like star simulated here, the photospheric temperature is set to T_phot = 5780 K. For the effective temperature of the spot area, including umbra and penumbra, I use T_spot = 5150 K as inferred from solar irradiance modeling by Fligge et al. (1998; see also, e.g., Fligge et al. 2001). More recent modeling may include umbral and penumbral components separately, such as in the study by Unruh et al. (2008), which uses a penumbral temperature of 5150 K and an umbral temperature of 4500 K. However, as the penumbral area dominates umbral areas by factors of typically 3:1 (Fligge et al. 1998) to 5:1 (Martinez Pillet et al. 1993), using a single, largely penumbral temperature as characteristic for the overall spot output suffices for the present purpose. With these values for T_phot and T_spot, the spot brightness relative to the photosphere at 3870 Å, 6010 Å, and 1.6 μm is 0.45, 0.60, and 0.79, respectively. The same relative limb-darkening curve is assumed as for the quiet photosphere at the wavelength under consideration.

Each flux concentration is mapped into the image pixel corresponding to its central location. For those with $| {{\rm{\Phi }}}_{i}| \leqslant {{\rm{\Phi }}}_{b}$, the effective radius A_eff is assumed to be r_c. For concentrations with $| {{\rm{\Phi }}}_{i}| \gt {{\rm{\Phi }}}_{b}$, an area of ${A}_{{\rm{\Phi }}}=| {{\rm{\Phi }}}_{i}| /{B}_{\mathrm{phot}}$, with B_phot = 1.35 kG is assumed to hold for the average field strength over the combined umbral and penumbral areas. For the concentrations for which A_Φ exceeds the unprojected surface area mapping into an image pixel, the area is approximated as a foreshortened circular disk (pixelation effects are corrected for by matching the total intensity in a disk to the area mapped onto an image array). For fractional pixels, and for concentrations for which the area A_eff is smaller than the solar surface area under an image pixel, the intensity in that pixel is modified using the fraction of the pixel area covered by the magnetic concentration, assuming the intensity of the complement of that pixel area is unaffected.

Simulated magnetograms are made by mapping the fluxes Φ_i of the flux concentrations onto an image, spread out over disks with radius A_Φ, and subsequently modified with a multiplicative factor for each pixel equivalent to assuming that the magnetic field is normal to the local surface.

Sample simulation results for the intensity images at three different wavelengths are shown in the right-hand three panels of Figures 1(a)–(d) (with doubled contrast to better show the faculae). The images in these figures were smoothed with a Gaussian with a FWHM value of 2 pixels. The images show that both the facular contrast and the spot contrast are largest for the shortest wavelengths, as expected from the facular contrasts and the wavelength-dependent contrast c_s(λ). And, also as expected, the faculae stand out most clearly toward the limb while spots and pores are most readily visible toward the central regions of the disks.

One property that jumps out in particular in the intensity images of Figures 1(c)–(d) is the abundance of dark features with diameters of order 8000 km. This is a consequence of the numerical prescription for the fractionation probability of flux concentration. This has a local minimum (and thus a corresponding maximum in the lifetime) around 10²¹ Mx. Consequently, clusters of such solar-like spots develop in the model. While this fractionation probability approximates the behavior of most sunspots well, it is as yet unknown how spots and spot clusters behave in very active stars or extremely large bipolar regions (if these indeed exist). Numerical experiments with radiative magnetoconvection by Beeck et al. (2015a, 2015b), however, suggest that in the presence of an increasingly strong field, (micro-)pores are a naturally occurring phenomenon for cool main-sequence stars when the mean flux density exceeds a few hundred Mx cm⁻². This may or may not be what happens on a very active star, but even if it does not happen, compact clusters of such spots in the model run here may effectively describe the appearance of any large spot in real life from the point of view of rotational modulation or transit signals for planets of a size comparable to Jupiter, as simulated here. Note that the images shown in Figures 1(a)–(d) have a best effective resolution of 5500 km at disk center (given the diameter of the rendered images of 512 pixels subjected to a Gaussian smoothing with FWHM of 2 pixels), to be compared to r_c = 4200 km.

Another feature that stands out in Figures 1(c)–(d) is the multitude of spots and pores that persist in the poleward arcs formed by decaying active regions, and also in flux dispersal toward the equator. In the models, these dark features are the result of the frequent collisions and temporary mergers between flux concentrations in regions of high average flux density. Such high flux concentrations on the Sun would correspond to pores and small spots.

First to discuss is the average relative offset ${\delta }_{* }\,=({I}_{* }(\infty ,\lambda )-{I}_{{\rm{q}}}(\lambda ))/{I}_{{\rm{q}}}(\lambda )$, where ${I}_{* }(\infty ,\lambda )$ is the brightness of the active star out of transit and I_q(λ) is that of the quiet-star (here quiet-Sun) photosphere. Values for δ_* are listed in Table 2. Note that these values should be compared across wavelength for a given set; substantial statistical fluctuations make these values less suitable for comparison between runs, although the trends are likely to stand. These values show that the tradeoff between bright faculae and dark spots is such that the overall brightness of a star with Sun-like activity (${ \mathcal M }({ \mathcal A }=1)$) is increased during cycle phases of enhanced activity relative to the reference quiet Sun. The averaged contrast over the three wavelengths, used as a rough proxy for the bolometric value, is ∼0.1%. This value compares well with the observed change in total solar irradiance over the cycle (e.g., Fröhlich 2016; Kopp 2016).

For the simulations of the more-active stars, the faculae play a relatively weaker role than starspots, resulting in a dimming of the overall stellar brightness. For ${ \mathcal M }{({ \mathcal A }=1)}^{+}$ this results in a near balance between spot darkening and facular brightening for the longer wavelengths, while for the most active star, simulated sunspot darkening outweighs facular contributions at all three model wavelengths. This trend is consistent with the empirical results of Lockwood et al. (1997) and Radick et al. (2018; see also Montet et al. 2017). In their long-term monitoring of a set of cool stars, they note that more-active stars tend to dim as chromospheric activity increases due to dynamo variability, suggesting a dominant role of starspots over faculae. In contrast, for many less-active stars—including the Sun—photospheric brightness increases with increasing activity, which is interpreted to be owing to a stronger role of faculae relative to that of spots, although the latter can cause dips in the brightness when crossing near central meridian due to solar rotation.

Figure 15 by Radick et al. (2018) suggests that spots and pores dominate in irradiance variations beyond a chromospheric activity level of $\mathrm{log}({R}_{{\rm{H}}+{\rm{K}}}^{{\prime} })\approx -4.75$, compared to a characteristic solar value of around −4.9. Converting the level of $\mathrm{log}({R}_{{\rm{H}}+{\rm{K}}}^{{\prime} })\approx -4.75$ into an S-index value (Noyes et al. 1984), and applying Equations (1) and (2) in Shapiro et al. (2014) yields a corresponding spot/pore filling factor of 1.2% and a facular filling factor of 6.8%, which with B = 1.35 kG yields $\langle | {fB}| \rangle \approx 100$ G for spots/pores and faculae combined. This appears compatible with the summary results in Table 2.

To compute the rotational modulation as shown by the solid curves in Figures 2–5, the perspective of the rendered photospheres is changed, but the magnetic field is held frozen in time (phase angle 0° corresponds to the images shown in Figures 1(a)–(d)). Dashed curves in these figures show the total magnetic flux on the "observer-facing side" of the star, rescaled to the same range for comparison, and inverted in sign because more flux is generally due to more active regions that put more spots and pores on the star, except in run ${ \mathcal M }({ \mathcal A }=1)$, where brightness and magnetic flux on a hemisphere are positively correlated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This single snapshot of rotational modulation clearly is only one of many instantiations of the field distribution, but gives an idea of the magnitude of the modulations that is sufficient for the present purpose. The amplitude of rotational modulation for the simulation of the Sun-like star (${ \mathcal M }({ \mathcal A }=1)$) at 4000 Å is ∼0.06%. This single example of rotational modulation lies only somewhat below the value of ≈0.1%–0.2% characteristic of the active Sun (Fröhlich 2016; Lee et al. 2016).

On the one hand, this single snapshot model is at least roughly of comparable magnitude to observations. On the other hand, the rotational modulation of more-active stars also comes out low if only the emergence frequency, but not the size spectrum, of active regions is changed. For ${ \mathcal M }({ \mathcal A }=30)$, the rotational modulation at visible wavelengths is ∼0.2%. This is considerably weaker than what is observed in stars of correspondingly high activity. One example is the rotational modulation observed for Kepler-17 (Table 1), which is about 1%–2%. Other examples of stars with comparable rotation periods can be found in Notsu et al. (2013) for a sample of stars that exhibit superflaring in the Kepler data: typical amplitudes of rotational brightness modulation for such stars fall in the range from 1% to 5%, with most clustering around 2%.

There are two ways by which the model's rotational modulation can be increased: raise the frequency of emergence of active regions, or allow for larger active regions to emerge. The frequency spectrum for emerging bipoles used in the present simulations is based on the work by Harvey & Zwaan (1993), and is thus characteristic of sunspot Cycle 21. Without guidance on, for example, preferred longitudes for flux emergence, statistically, one expects an increase in rotational modulation that scales as the square root of the number of regions (see also Rackham et al. 2018). Based on this, rather than modifying this frequency or introducing longitudinal patterns in an ad hoc fashion, the second possibility is explored. The largest active regions that emerged in the 846-day period of observations analyzed by Harvey & Zwaan (1993) had a total absolute magnetic flux per polarity of roughly 4×10²² Mx, with only a single occurrence of such a large region during the period of their analysis. This led Schrijver (2001) and Schrijver & Title (2001) to set Φ_max = 1.5×10²² Mx in their model runs to stay just below the most uncertain part of the active-region size spectrum with infrequent realizations on the Sun.

A later study by Zhang et al. (2010), spanning a longer period, showed regions with fluxes up to about 3×10²³ Mx per polarity. When the flux-emergence spectrum is extended to this flux for a star of otherwise solar properties (in run ${ \mathcal M }{({ \mathcal A }=1)}^{+}$), the rotational modulation at visible wavelengths reaches ∼0.3%–0.4%. This certainly lies in the range of solar rotational modulation, but now the added flux leads to surface-averaged flux densities that are too high (Table 2). Without experimenting further in the present context, I hypothesize that there may be a decline in active-region frequencies for region with fluxes approaching 3×10²³ Mx, dropping below the power law seen at lower fluxes. This would be in line with comments by Zhang et al. (2010), and with the occurrence of infrequent but large regions discussed by Schrijver et al. (2012), while occasionally resulting in larger rotational modulation and also regions large enough for infrequent powerful solar flares (Aulanier et al. 2013).

The result of allowing larger regions is also noticeable for the more-active stars: the rotational modulation found for ${ \mathcal M }{({ \mathcal A }=10)}^{+}$ and ${ \mathcal M }{({ \mathcal A }=30)}^{+}$ with ${{\rm{\Phi }}}_{\max }=3\times {10}^{23}$ Mx is significantly larger than for ${ \mathcal M }({ \mathcal A }=30)$, reaching roughly 1%–3% in the visible. Not only are these models in the range of the observed rotational modulation of Kepler-17 and other similar Kepler targets, but they also inject bipoles of sufficient magnitude to power the occasional superflares observed on such stars (e.g., Schrijver et al. 2012; Aulanier et al. 2013; Notsu et al. 2013).

Model ${ \mathcal M }{({ \mathcal A }=30)}^{+}$ has an average magnetic flux density of $\langle | {fB}| \rangle \sim 500$ G (Table 2). Notsu et al. (2015) show a trend of increasing $\langle | {fB}| \rangle$ versus the amplitude of the rotational brightness variation (their Figure 8(b)). They base this scaling on solar observations of the Ca ii infrared triplet compared to a magnetogram, and measure the stellar Ca ii infrared triplet to establish the mean stellar level of activity. If only the total flux in facular elements (with fluxes below 3×10²⁰ Mx) is taken to contribute in the Ca ii triplet, then for model values of $\langle | {fB}| {\rangle }_{\mathrm{fac}}\approx 180$ G the scaling found by Notsu et al. (2015) would map to brightness variation amplitudes by rotational modulation of ≈1.5%, which compares well with what is found here (Table 2).

The rotational modulation is the result of the inhomogeneous distribution of magnetic flux across the stellar surface, with the signal presenting a mean over entire hemispheres. The result is a relatively smooth curve as a function of the phase angle, as shown in Figures 2–5. The hemispheric smoothing can result in a simple modulation often at the apparent period of rotation (Figures 4 and 5), but for the less-active stars, rotation periods can be masked. The appearance of observed light curves of active stars may suggest the existence of only a few large spot groups, but the images in Figures 1(a)–(d) suggest that the real situation may be far more complex.

Figure 6 illustrates the rotational modulation over time by looking at the total amount of magnetic flux on the observer-facing hemisphere of the "rotating star": here, the magnetic evolution is shown for run ${ \mathcal M }{({ \mathcal A }=10)}^{+}$ over a total of 20 simulated rotation periods set to 12 days each, at the full model resolution of its 6 hr steps, centered on the time shown in Figure 1(c). The figure shows a clear long-term evolution in the total amount of flux, but for most of the period, there is a pronounced rotational modulation suggestive of a dominant active hemisphere where none is prescribed in the model.

Zoom In Zoom Out Reset image size

The average radius of the selected exoplanets in Table 1 is 0.14R_*, which is also characteristic of WASP-36b and Kepler-17b, which orbit the most Sun-like stars of the longest rotation periods in the sample listed. In the transit simulations discussed here, I therefore adopt a planetary radius of 0.14 R_⊙. Smaller virtual exoplanets can, of course, also be used, but that is an application that can await both advanced magnetoconvective models of stellar photospheres and very large aperture telescopes with increased S/N properties. One exception is made for comparison to a transit with a planetary radius of 0.02 R_⊙ discussed below.

Figures 1(a)–(d) show examples of transits paths across the model stars with a relative impact parameter for the transit of d_t = 0.25R_* (compare with the estimated values for the selected stellar sample in Table 1), for which Figures 7–10 show the transit residuals, calibrated relative to the inactive star and to zero out of transit:

$$\begin{eqnarray}&&{ \mathcal R }(t,\lambda )\equiv \displaystyle \frac{({I}_{* }(t,\lambda )-{I}_{{\rm{q}}}(t,\lambda ))}{{I}_{{\rm{q}}}(t,\lambda )}-\displaystyle \frac{({I}_{* }(\infty ,\lambda )-{I}_{{\rm{q}}}(\infty ,\lambda ))}{{I}_{{\rm{q}}}(\infty ,\lambda )},\end{eqnarray} \tag{ 1 }$$

where ${I}_{* }(t,\lambda )$ is the modeled transit intensity at time t (running in the figures from first contact with the image field of view to fourth contact) and wavelength λ, and ${I}_{* }(\infty ,\lambda )$ is the brightness of the active star out of transit; I_q are the same for the quiet, featureless photosphere. The curves are fairly smooth for two reasons: the random-walk diffusion of the magnetic flux leads to relatively smooth gradients in the surface field, with further smoothing imposed by the exoplanet diameter. But there are fine-scale variations that, if they are detectable in real observations, could be used to study structural details of stellar magnetic fields, provided that the effective telescope aperture is sufficiently large to keep exposures short while obtaining a high S/N.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 7–10 show only a single instance of an arbitrarily chosen step in the flux-transport simulations and a single viewing angle, while comparing the residuals for the three wavelengths used here. Figure 11 provides an impression of the full range of transit curves (top panel) and transit residuals (bottom panel) for run ${ \mathcal M }{({ \mathcal A }=10)}^{+}$ (compatible with a Sun-like star at P_rot ≈ 12 days) for only 3870 Å computed by rotating the star underneath a series of transits. The full range of transit residuals for this time step clearly shows the effect of limb brightening by faculae in the dips early and late in the transits. Moreover, the full range of residuals is more than double the amplitude of the single realization selected for Figure 9. For comparison, Figure 12 shows the simulated transit signals and residuals for a much smaller planet (${R}_{{\rm{p}}}/{R}_{* }=0.02$); note that for such a small planet, umbral–penumbral differences—here ignored—are in principle accessible, although this would require S/N levels that are currently unaccessible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All transit residuals reveal signatures of the facular brightening toward stellar limbs, recognized as dips near the beginnings and ends of the transit residuals (to be compared, e.g., to a similarly shaped transit light curve for Kepler-71 shown in Figure 4 of Zaleski et al. 2019). This is in effect a wavelength- and activity-dependent distortion of the limb-darkening curve (which should be anticipated to be asymmetric because the stellar activity along the transit path is expected to be generally asymmetric).

Next to be considered is the magnitude of the residuals: the peak-to-trough ranges in the transit residuals, Δ(λ) (summarized in Table 2), confirm the very weak transit signals expected for a star of solar activity, ${ \mathcal M }({ \mathcal A }=1)$, even when allowing for large active regions as in ${ \mathcal M }{({ \mathcal A }=1)}^{+}$. For the stars listed in Table 1 with rotation periods in the range of 8–12 days, transit-residual amplitudes reach values of up to 0.15%–0.4%. These values are much larger than those found for ${ \mathcal M }({ \mathcal A }=30)$, but are compatible with the model results for ${ \mathcal M }{({ \mathcal A }=10)}^{+}$ and ${ \mathcal M }{({ \mathcal A }=30)}^{+}$, i.e., for models in which much larger active regions are included in the emergence spectrum.

Figures 1(c) and (d) reveal the consequences of raising Φ_max to 3×10²³ Mx: the injection of very large regions leads to a larger length scale in the flux patterns (compared to solar patterns as in Figure 1(a)), and these large areas of high magnetic flux densities are conducive to forming or maintaining concentrations with high magnetic fluxes, many of which appear as spots and pores. At least that is what happens in the models based on the properties derived from solar observations.

The effective area of the occulting disk relative to that of the background star, based on the relative occultation depth $d(t,\lambda )=1-{I}_{* }(t,\lambda )/{I}_{* }(\infty ,\lambda )$ at any time t during the transit can be reformulated as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{p}}}^{2}(\lambda )}{{R}_{* }^{2}(\lambda )}=d(t,\lambda )\left(\displaystyle \frac{{I}_{* }(\infty ,\lambda )}{{I}_{{\rm{o}}}({\boldsymbol{r}}(t),\lambda )}\right),\end{eqnarray} \tag{ 2 }$$

where ${I}_{{\rm{o}}}({\boldsymbol{r}},\lambda )$ is the mean intensity of the patch occulted by the transiting exoplanet and ${I}_{* }(\infty ,\lambda )$ is the overall mean stellar intensity out of transit. Below, any wavelength dependence of R_* is ignored.

Determining the deceptively simple-looking ratio ${I}_{* }(\infty ,\lambda )/{I}_{{\rm{o}}}({\boldsymbol{r}}(t),\lambda )$ is a challenge: we do not know what is masked by the exoplanet, and what is important in the ratio is the difference between the surface structures on the entire observer-facing side of the star and that occulted by the planet, including limb darkening, which is relatively poorly known compared to that to which we have some observational access, namely along the transit path, and there is contaminated by surface activity.

Having observations at multiple wavelengths will help in constraining spot and plage contributions to the stellar brightness over the spectrum, but not without more study. For example, the observed intensities I_*(λ) are related to the intrinsic intensity of a non-active star, I_q(λ), the spot contrast and filling factor, c_s(λ) and f_s, and the facular brightness contrast across the surface, ${{ \mathcal C }}_{{\rm{f}}}({\boldsymbol{r}},\lambda )$, integrated over the observer-facing hemisphere, here with total area normalized to unity:

$$\begin{eqnarray}&&{I}_{\ast }(\lambda )={I}_{{\rm{q}}}(\lambda )\left(1+{f}_{{\rm{s}}}{c}_{{\rm{s}}}(\lambda )+\int {{ \mathcal C }}_{{\rm{f}}}({\bf{r}},\lambda ){\rm{d}}S\right).\end{eqnarray} \tag{ 3 }$$

If the facular contrasts would have been essentially multiplicatively scaling functions across wavelength, i.e., if ${{ \mathcal C }}_{{\rm{f}}}({\boldsymbol{r}},\lambda )\approx {C}_{{\rm{f}}}(r){c}_{{\rm{f}}}(\lambda )$ for a contrast C_f at some reference wavelength associated with faculae only, then

$$\begin{eqnarray}&&{I}_{\ast }(\lambda )\approx {I}_{{\rm{q}}}(\lambda )\left(1+{f}_{{\rm{s}}}{c}_{{\rm{s}}}(\lambda )+{c}_{{\rm{f}}}(\lambda )\int {C}_{{\rm{f}}}({\bf{r}}){\rm{d}}S\right),\end{eqnarray} \tag{ 4 }$$

with the integral, independent of wavelength, a constant to be determined. Then the problem of determining the unknowns would have been a straightforward inversion problem. But the fact that the limb-brightening curves for faculae modeled by Norris et al. (2017) do not transform into one another by simple multiplicative factors implies that in order to fully disentangle the spot and facular components, the spatial distribution of these components across the stellar disk needs to be reasonably well known. For some purposes and some wavelengths, the approximation in Equation (4) may suffice, but determining this is beyond the scope of this paper.

An estimate of the relative impact of the effect of the activity in the transit chord on the radius of the exoplanet and any surrounding atmosphere can be obtained by writing ${I}_{{\rm{o}}}({\boldsymbol{r}}(t),\lambda )={I}_{{\rm{q}}}({\boldsymbol{r}}(t),\lambda )(1+{\delta }_{{\rm{o}}}({\boldsymbol{r}}(t),\lambda ))$, where I_q(${\boldsymbol{r}})$(t), λ) is the quiet, limb-darkened photosphere behind the occulting planet and ${\delta }_{{\rm{o}}}({\bf{r}}(t),\lambda )$ is the (spot plus faculae) relative difference in intensity of the occulted patch due to stellar activity. Similarly, we can write ${I}_{* }(\infty ,\lambda )\,={I}_{{\rm{q}}}(\infty ,\lambda )(1+{\delta }_{* }(\lambda ))$. For relatively small perturbations, so if ${\delta }_{{\rm{o}}}({\bf{r}}(t),\lambda )\ll 1$ and δ_*(λ) ≪ 1, then Equation (2) transforms into

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{R}_{{\rm{p}}}(\lambda ) & \approx & {R}_{\ast }(\lambda )\,{d}^{\displaystyle \frac{1}{2}}(t,\lambda ){\left(\displaystyle \frac{{I}_{{\rm{q}}}(\infty ,\lambda )}{{I}_{{\rm{q}}}({\boldsymbol{r}}(t),\lambda )}\right)}^{\displaystyle \frac{1}{2}}\\ & & \times \,\left(1+\displaystyle \frac{1}{2}{\delta }_{\ast }(\lambda )-\displaystyle \frac{1}{2}{\delta }_{{\rm{o}}}(r(t),\lambda )\right).\end{array}\end{array}\end{eqnarray} \tag{ 5 }$$

The final part of this equation essentially quantifies the "transit light source effect" related to stellar magnetic activity (e.g., Rackham et al. 2018, 2019b).

Ignoring measurement uncertainties and assuming the quiet-star photospheric brightness, including limb darkening, is known, the final expression between parentheses in Equation (5) quantifies the effect of stellar activity on the estimated planetary radius. Figure 13 shows an example of the overall impact on estimated values of R_p(λ) for a transiting exoplanet with radius 1.4R_Jup for ${ \mathcal M }{({ \mathcal A }=10)}^{+}$ by the last term in the final parenthetical expression of Equation (5): in the near-IR, the impact is of order 2000–4000 km, while differences across wavelengh from the blue to the near-IR amount to at least 1000 km even when times of minimal spot coverage are selected during the transit (thin gray lines in Figure 13).

Zoom In Zoom Out Reset image size

If the wavelength-dependent impact of faculae and spots on overall brightness were known, i.e., if we knew δ_*(λ)—which we do in this case for the models—then this reduces the impact on the estimated planet radius in this example (thick black lines in Figure 13), depending on whether the transit covers sufficiently inactive photospheric regions. If, however, one were to use the amplitude of the rotational modulation along with the assumption that the asymmetry between hemispheres is entirely attributable to a single dominant spot group, and were to ignore facular contributions, then the estimated values for δ_*(λ) would be substantially smaller than the true values, leading to wavelength-dependent errors equivalent to several 1000 km. For the other models in Table 2, these correction terms impact radius estimates from some 500 km in stars with solar-like activity (${ \mathcal A }=1$) to of order 5000 km at the high-activity end (${ \mathcal A }=30$).