© K. L. Yeo et al. 2021

1. Introduction

The global solar dynamo produces the cyclic emergence of bipolar magnetic regions (BMRs) on the solar surface which underlies the 11-year activity cycle (Charbonneau 2020). Emergent BMRs range from the larger sunspot-bearing magnetic bipoles, which go on to form active regions, to the smaller spot-free and relatively short-lived (timescale of hours) ephemeral regions (ERs) (Harvey 1993, 2000; Hagenaar 2001; Hagenaar et al. 2003; Thornton & Parnell 2011). How the larger BMRs relate to the sunspots they enclose and how ER emergence relates to the number of sunspots present are of interest as they connect directly to how solar magnetism and activity relate to sunspot prevalence, with implications for the study of a wide range of solar phenomena.

The large-scale photospheric magnetic field evolves as BMRs emerge, their magnetic flux is transported and dispersed by systematic and turbulent motions, and magnetic elements of opposite polarity approach and cancel each other (van Driel-Gesztelyi & Green 2015). This process can be described via the surface flux transport equation (Leighton 1964; DeVore et al. 1984) or, at a simpler level, by a closed system of ordinary differential equations (Solanki et al. 2000, 2002). Models have been developed to reproduce the past and forecast the future evolution of the large-scale photospheric magnetic field by solving the transport equation for a given set of initial conditions and details of BMR emergence, aptly termed surface flux transport models (SFTMs) (Jiang et al. 2014). Some studies have simulated the evolution over the past few centuries by prescribing to SFTMs the details of BMR emergence over this period, as determined from sunspot number or sunspot area records (e.g. Wang et al. 2005; Cameron et al. 2010; Jiang et al. 2011a,b). Sunspot number and sunspot area records are employed as they are the only direct observations of solar surface magnetism to go this far back in time. Otherwise, SFTMs usually make use of solar magnetograms (e.g. Whitbread et al. 2018; Virtanen et al. 2018; Yeates 2020) or Ca II K spectroheliograms (e.g. Virtanen et al. 2019) instead. While magnetogram and Ca II K spectroheliogram observations do not extend as far back in time as sunspot observations, the details of BMR emergence can be established more directly and accurately from such data.

Since the large-scale photospheric magnetic field is the observable part of the global dynamo and determines the structure of the coronal magnetic field, SFTMs represent a powerful tool in the study of these solar phenomena (Cameron et al. 2012; Mackay & Yeates 2012; Jiang et al. 2013). Models of solar irradiance variability are often employed to generate the historical solar forcing input required by climate simulations (Haigh 2007) and to investigate how the brightness variability of the Sun compares to that of other cool stars (Shapiro et al. 2014, 2015, 2016; Nèmec et al. 2020a,b). Through their incorporation into models of solar irradiance variability, discussed next, SFTMs are also relevant to the study of the Earth’s climate and of the solar-stellar connection.

Solar irradiance variability on timescales greater than a day is believed to be mainly driven by the variation in the prevalence and distribution of dark sunspots and bright faculae and network on the solar disc as the photospheric magnetic field evolves (Solanki et al. 2013; Yeo et al. 2017). Accordingly, solar irradiance models describe the variation in this quantity due to the intensity deficit from sunspots and the excess from faculae and network, and this information is determined from observations of solar magnetism (Domingo et al. 2009; Yeo et al. 2014a).

The variation in solar irradiance over the past few centuries, a period of particular interest to climate studies, is often modelled from sunspot number and sunspot area records (e.g. Lean 2000; Krivova et al. 2007, 2010; Dasi-Espuig et al. 2014, 2016; Coddington et al. 2016; Wu et al. 2018). Modelling solar irradiance variability with sunspot observations means inferring not just the effect of sunspots, but also the effect of faculae and network on solar irradiance from information about sunspots alone. This can be aided by SFTMs aimed at reproducing the evolution of the large-scale photospheric magnetic field from the same sunspot observations. For example, Coddington et al. (2016) incorporated into their solar irradiance reconstruction the decadal variability in the SFTM reported by Wang et al. (2005). In another instance, the solar irradiance reconstruction by Dasi-Espuig et al. (2014, 2016) is based on synthetic full-disc magnetograms generated from sunspot area and sunspot number records using the SFTM of Cameron et al. (2010) and Jiang et al. (2011a,b) (hereafter CJSS).

Various studies found the facular area to increase with increasing sunspot area at a diminishing rate such that it appears to gradually saturate (Foukal 1993, 1996, 1998; Chapman et al. 1997; Shapiro et al. 2014). Solanki & Fligge (1998) and Yeo et al. (2020) (hereafter YSK), compared various disc-integrated quantities that are sensitive to the amount of faculae and network present to sunspot area and/or sunspot number and returned analogous results. As first pointed out by Foukal (1993), these results suggest that the balance between the competing effects of bright faculae and network and of dark sunspots on solar irradiance might shift towards the latter as we go from low to high activity levels. This is believed to extend to cool stars in general, such that as we go from less to more active stars, we transit from a faculae-dominated regime where facular and network brightening dominates starspot darkening to a spot-dominated regime where the converse is true. This explains why less active stars appear to become brighter with increasing activity, while more active stars exhibit the opposite trend (Lockwood et al. 2007; Hall et al. 2009; Shapiro et al. 2014; Radick et al. 2018). The balance between bright faculae and network and dark sunspots on the solar disc, and how that might change with the activity level and the inclination of the Sun to the observer, is of interest for the relevance to the understanding of the brightness variability of the Sun and other cool stars. In recent studies Nèmec et al. (2020a,b) extended the model of solar irradiance variability by Dasi-Espuig et al. (2014, 2016), which incorporates the SFTM by CJSS, to probe the inclination-dependence of this balance.

In order to reproduce the evolution of the large-scale photospheric magnetic field over a given period of time, SFTMs require as input the physical properties of each BMR that had emerged over this time interval, including the timing, location, and amount of magnetic flux at emergence. However, deriving this information from observations of sunspots is less than straightforward. As ERs do not contain or emerge synchronously with sunspots and are difficult to observe due to their relatively small sizes and short lifetimes, the relationship between ER emergence and sunspot prevalence is not exactly known. As for sunspot-bearing BMRs, it is less than straightforward to infer the details of their emergence from sunspot number records as these are essentially sunspot and sunspot group counts, and therefore do not contain information on individual sunspot groups. While sunspot area records do tabulate the location and area of individual sunspot groups, how the amount of magnetic flux in a given newly emerged sunspot-bearing BMR, Φ_BMR, depends on the area of the enclosed sunspots, A_BMR,S, is not well constrained. By the term ‘newly emerged’ we refer to the point where the sunspot-bearing BMR has fully emerged but not yet started to decay. This paper focuses on the issue of the relationship between Φ_BMR and A_BMR,S. For the sake of brevity, hereafter the term BMR denotes sunspot-bearing BMRs alone, excluding ERs.

The difficulty in establishing BMR emergence from sunspot observations leaves SFTMs to make varying assumptions. The Wang et al. (2005) model, based on the Hoyt & Schatten (1998) sunspot group number record, accommodated the fact that there is no information on individual sunspot groups in such data by assigning the same amount of magnetic flux to each of the emergent BMRs in a given solar cycle. In the CJSS model, detailed in Sect. 3.1, the A_BMR,S of each BMR is taken directly from sunspot area records. We let Δβ denote the separation between the two poles of a BMR. The model combines an empirical relationship between Δβ and A_BMR,S and the assumption that Φ_BMR is linearly proportional to (Δβ)² to determine Φ_BMR from A_BMR,S.

Few studies have attempted to constrain the relationship between Φ_BMR and A_BMR,S. From an examination of the faculae in active regions and newly emerged BMRs, Schrijver & Harvey (1994) postulated that Φ_BMR scales linearly with the total area extended by the faculae and sunspots in each BMR, A_BMR. This can be combined with the empirical relationship between facular area and sunspot area (Foukal 1993, 1996, 1998; Chapman et al. 1997; Shapiro et al. 2014) to roughly relate Φ_BMR to A_BMR,S. In a more recent study, Lemerle et al. (2015) compared the database of newly emerged BMRs compiled by Wang & Sheeley (1989) to the USAF-NOAA sunspot area record and found Φ_BMR to conform to a power-law relationship with A_BMR,S with a power-law exponent of 0.71. The analyses of Schrijver & Harvey (1994) and Lemerle et al. (2015) (hereafter SH and LCC), were restricted by certain limitations of the observations available to them, discussed in Sect. 2.3. It would be straightforward to derive the relationship between Φ_BMR and A_BMR,S from an extensive and reliable record of the magnetic and physical properties of newly emerged BMRs. There is, however, no such record currently available, as also discussed in Sect. 2.3.

In this study we aim to establish the relationship between the amount of magnetic flux in each newly emerged BMR, Φ_BMR, and the area of the enclosed sunspots, A_BMR,S. We make use of the empirical model relating the disc-integrated facular and network magnetic flux, ΣΦ_FN, to the total surface coverage by sunspots, ΣA_S presented by YSK (Yeo et al. 2020), which builds on the earlier work by Preminger & Walton (2005, 2006a,b, 2007), abbreviated here as PW. As we demonstrate, the set-up of the YSK model allows us to establish from it an empirical relationship between Φ_BMR and A_BMR,S while circumventing the lack of a proper database of newly emerged BMRs. In the following, we describe the YSK model (Sect. 2.1), derive the empirical relationship between Φ_BMR and A_BMR,S (Sect. 2.2), and compare it to the SH and LCC results (Sect. 2.3). Then, we apply each of these constraints on Φ_BMR to the SFTM by CJSS (Sect. 3.1) and examine the effect on the model output (Sect. 3.2). Making use of the surface flux transport simulation based on the empirical relationship between Φ_BMR and A_BMR,S derived here, we will look at the balance between the amount of faculae and network and the amount of sunspots on the solar disc, and how that changes with solar activity level (Sect. 4). Finally, we provide a summary of the study in Sect. 5.

2. The relationship between the amount of magnetic flux in newly emerged bipolar magnetic regions and the area of the enclosed sunspots

2.1. Yeo et al. (2020) model

In this section we derive an empirical relationship between the amount of magnetic flux in each newly emerged bipolar magnetic region or BMR, Φ_BMR, and the area of the enclosed sunspots, A_BMR,S. For this purpose we employ the empirical model of the relationship between the disc-integrated facular and network magnetic flux, ΣΦ_FN, and the total surface coverage by sunspots, ΣA_S, recently reported by YSK (Yeo et al. 2020). The YSK model is given by

$$\begin{array}{c}\hfill \mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}={\left(\mathrm{\Sigma }{A}_{\mathrm{S}}\right)}^{h_1}\otimes H+h_2,\end{array}$$(1)

where h₁ and h₂ are fit parameters. The convolution of (ΣA_S)^h₁ with H, a finite impulse response (FIR) filter, describes the variation in ΣΦ_FN associated with sunspot-bearing BMRs. The h₂ term represents network magnetic flux that did not originate from sunspot-bearing BMRs; it is assumed in the model to be constant. The FIR filter is given by

$$\begin{array}{c}\hfill H(t)=h_{3}max[exp(-\frac{|t|}{h_4})cos\left(\frac{2\pi t}{t_{\odot }}\right),0],\end{array}$$(2)

where −0.25 t_⊙ ≤ t ≤ 14.25 t_⊙, t_⊙ = 26.24 days denotes the synodic rotation period of the Sun, and h₃ and h₄ are fit parameters.

Yeo et al. (2020) extended the composite time series of daily ΣA_S by Balmaceda et al. (2009) (Fig. 1a) and of daily ΣΦ_FN by Yeo et al. (2014b) (Fig. 1b) to June 2017 and constrained the fit parameters with the result. Yeo et al. (2014b) had derived daily ΣΦ_FN back to 1974 by taking daily full-disc magnetograms from the Kitt Peak Vacuum Telescope (Livingston et al. 1976; Jones et al. 1992), SoHO/MDI (Scherrer et al. 1995), and SDO/HMI (Scherrer et al. 2012), accounting for the systematic differences between the observations from the various telescopes due to instrumental factors and isolating the faculae and network present in each magnetogram. YSK could not extend the Balmaceda et al. (2009) time series beyond June 2017 as it is based on USAF-NOAA sunspot areas for the period of 1986 on, and these measurements are available only up to this time. For ΣA_S in millionths of the solar hemisphere (μHem) and ΣΦ_FN in units of Weber (Wb), h₁ = 0.843, h₂ = 1.82 × 10¹⁴ Wb, h₃ = 1.99 × 10¹¹ Wb, and h₄ = 33.5 days. The robustness of the model is visibly evident in Fig. 1b; the model reconstruction of ΣΦ_FN from ΣA_S (red) is very similar to the Yeo et al. (2014b) time series (black), and the 365-day running mean of the difference between the two (blue) is close to null and does not indicate any clear systematic departures. The Pearson’s correlation coefficient, R, between the two time series is 0.950, meaning the model reproduces about 90% of the observed variability. We refer the reader to Yeo et al. (2020) for a detailed discussion of the uncertainties in the model and how it compares to earlier models.

Fig. 1. Observed ΣAS and ΣΦFN, YSK model (Eq. (1)) reconstruction of ΣΦFN from ΣAS and H term in this model (Eq. (2)). (a) Composite time series of the total surface coverage by sunspots, ΣAS, by Balmaceda et al. (2009). (b) Composite time series of the disc-integrated facular and network magnetic flux, ΣΦFN, by Yeo et al. (2014b) (black), the reconstruction of this quantity from ΣAS with the YSK model (red), and the 365-day running mean of the difference between measurement and model (blue). The Balmaceda et al. (2009) time series goes back to 1874 but the range of the plot is confined to 1974 to 2017, the period where it overlaps with the Yeo et al. (2014b) time series. The gaps in the Yeo et al. (2014b) time series correspond to periods where there no suitable magnetogram observations with which to determine ΣΦFN were available. (c) The H term, which describes the response of ΣΦFN to the appearance of a sunspot per unit sunspot area. The dotted line follows the exponential envelope of H and t⊙ denotes the synodic rotation period of the Sun (see Sect. 2.1 for details).

The FIR filter H (Eq. (2)) is based on the earlier findings of PW (Preminger & Walton 2005, 2006a,b, 2007). Sunspots form some time after the BMRs they belong to have started emerging on the solar surface; they are short-lived compared to active regions and their decay products, such that the latter persists long after their sunspots have decayed (van Driel-Gesztelyi & Green 2015). Consequently, the response of any quantity that is sensitive to the presence of BMRs, active regions, and their decay products, including ΣΦ_FN, to the appearance of a sunspot similarly starts some time before it appears and persists for a period after it dissipates. PW deconvolved ΣA_S from the average unsigned magnetic flux density over the solar disc, solar irradiance, and various chromospheric and coronal indices, deriving FIR filters that elucidate the response of each quantity to the appearance of a sunspot and how this response varies with time due to the phenomenon just described. Preminger & Walton (2007) captured the key features of these empirical FIR filters with a model FIR filter (Eq. (2) in their paper), which YSK modified slightly to more accurately reflect the form of the empirical FIR filters, yielding H. Depicted in Fig. 1c, H describes a series of sinusoidal lobes, 0.5 t_⊙ wide and t_⊙ apart, of diminishing amplitude, as delineated by the exponential envelope of time constant h₄ (dotted line). The first lobe has an amplitude of h₃ and is centred on t = 0. The sequence of lobes is truncated at t = 14.25 t_⊙. The FIR filter charts the amount of facular and network magnetic flux associated with a sunspot-bearing BMR, per unit sunspot area, from 0.25 t_⊙ before it fully emerges to 14.25 t_⊙ after. The rising edge of the first lobe corresponds to the BMR emerging, the peak of the first lobe to full emergence, and the falling edge and the subsequent lobes to its decay and modulation with solar rotation.

If the amount of facular and network magnetic flux associated with a sunspot-bearing BMR scales linearly with the area of the enclosed sunspots, then the variation in ΣΦ_FN due to sunspot-bearing BMRs would be given by ΣA_S ⊗ H. However, as noted in the introduction, various studies have found the facular area to increase with increasing sunspot area at a diminishing rate, such that it is at a lower level than if the two quantities are linearly related (Foukal 1993, 1996, 1998; Chapman et al. 1997; Shapiro et al. 2014). This implies that the larger an active region, the higher (lower) the proportion of the total area occupied by sunspots (faculae). Based on this observation, YSK argued that it is unlikely that the amount of facular and network magnetic flux associated with a sunspot-bearing BMR conforms to a linear relationship with the area of the enclosed sunspots. They introduced the fit parameter h₁, such that the variation in ΣΦ_FN is given by (ΣA_S)^h₁ ⊗ H instead. This and the modifications introduced to the model FIR filter reported by PW produced a closer agreement between modelled and measured ΣΦ_FN, detailed in Yeo et al. (2020).

After the time of the YSK study, the ΣA_S composite by Balmaceda et al. (2009) was revised by Mandal et al. (2020). Both the Balmaceda et al. (2009) and Mandal et al. (2020) time series are based on, up to 1976, Royal Greenwich Observatory (RGO) sunspot area measurements. They are therefore exactly identical to one another up to this time, and differ only in the observations used to extend the RGO time series to more recent times. Mandal et al. (2020) noted that these two time series differ from one another by about 6% between 1977 and 1985, and that there are otherwise no noticeable systematic differences between the two. The revision of the YSK model with the Mandal et al. (2020) time series is outside the scope of the current study. In any case, we do not expect the result of applying the Mandal et al. (2020) time series to the YSK model to differ significantly from what was obtained with the Balmaceda et al. (2009) time series as input. The disparity between these two time series between 1977 and 1985 is relatively small, especially in comparison to the random error in sunspot area measurements (20% to 50%, Sofia et al. 1982). In addition, less than 10% of the data points in the daily ΣΦ_FN composite by Yeo et al. (2014b) lie in this period.

2.2. Derivation

The total magnetic flux in each newly emerged BMR is the sum of the magnetic flux in its faculae and in its sunspots, or

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BMR}}={\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{F}}+{\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{S}}.\end{array}$$(3)

Here we establish how Φ_BMR,F and Φ_BMR,S relate to A_BMR,S in turn. To derive the relationship between Φ_BMR,F and A_BMR,S, let us isolate from the YSK model the contribution to ΣΦ_FN by the faculae in newly emerged BMRs, ΣΦ_BMR,F. As noted in Sect. 2.1, the peak of the first lobe of the FIR filter, which is at t = 0, corresponds to the response of ΣΦ_FN to a sunspot-bearing BMR when it has fully emerged and has not started to decay. Taking Eq. (1), we isolate ΣΦ_BMR,F by replacing H with H(t = 0) = h₃ and factoring out the basal level, h₂, arriving at

$$\begin{array}{c}\hfill \mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{F}}=h_3{\left(\mathrm{\Sigma }{A}_{\mathrm{S}}\right)}^{h_1}.\end{array}$$(4)

The PW and YSK approach, by modelling the sunspot-related variation in a given disc-integrated quantity as the convolution of ΣA_S with a FIR filter and with the various quantities at daily cadence, treats the appearance of sunspots as one-day events. In other words, it implicitly assumes that each sunspot exists for just one day, such that the total surface coverage by sunspots, ΣA_S, is treated in the model as the total surface coverage by sunspots in newly emerged BMRs, ΣA_BMR,S. Replacing ΣA_S in Eq. (4) with ΣA_BMR,S, we get

$$\begin{array}{c}\hfill \mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{F}}=h_3{\left(\mathrm{\Sigma }{A}_{\mathrm{BMR},\mathrm{S}}\right)}^{h_1}.\end{array}$$(5)

Assuming this relationship applies to individual BMRs, then

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{F}}& =h_3{A}_{\mathrm{BMR},\mathrm{S}}^{h_1}\hfill \\ \hfill & =1.99\times {10}^{11}{A}_{\mathrm{BMR},\mathrm{S}}^{0.843}\mathrm{Wb}.\hfill \end{array}$$(6)

Cho et al. (2015) examined how the amount of magnetic flux in sunspots and pores varies with their size. The authors noted that for fully formed sunspots, the amount of magnetic flux appears to scale linearly with its area, confirming the early result by Sheeley (1966). Assuming the mean magnetic flux density over sunspots is 1350 G or 4.11 × 10¹¹ Wb μHem⁻¹ (Solanki & Schmidt 1993; Lites et al. 1993; Keppens & Martinez Pillet 1996), this means

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{S}}=4.11\times {10}^{11}{A}_{\mathrm{BMR},\mathrm{S}}\mathrm{Wb}.\end{array}$$(7)

Taking Eqs. (6) and (7) into Eq. (3), we get

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\mathrm{BMR}}& =4.11\times {10}^{11}{A}_{\mathrm{BMR},\mathrm{S}}+1.99\times {10}^{11}{A}_{\mathrm{BMR},\mathrm{S}}^{0.843}\mathrm{Wb}\hfill \\ \hfill & =4.11\times {10}^{11}(A_{\mathrm{BMR},\mathrm{S}}+0.485A_{\mathrm{BMR},\mathrm{S}}^{0.843})\mathrm{Wb}\hfill \\ \hfill & =4.11\times {10}^{11}{\stackrel{\sim }{A}}_{\mathrm{BMR}}\mathrm{Wb},\hfill \end{array}$$(8)

which relates the amount of magnetic flux in newly emerged BMRs to the area of the enclosed sunspots. The ${\stackrel{\sim }{A}}_{\text{BMR}}$ term, representing $A_{\text{BMR},\text{S}}+0.485A_{\text{BMR},\text{S}}^{0.843}$, denotes for a given BMR the area of the sunspot that would contain the same amount of magnetic flux (compare Eqs. (7) and (8)). It is worth noting here that the YSK model (Sect. 2.1) and the analysis based on this model that we have just presented pertains only to sunspot-bearing BMRs and does not apply to ephemeral regions.

In Fig. 2 we chart Φ_BMR,S, Φ_BMR (dashed and solid black lines, top panel), and the ratio of the two, $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ (black line, bottom panel) as a function of A_BMR,S. As is visibly evident, the empirical relationship between Φ_BMR and A_BMR,S presented in Eq. (8) indicates that most of the magnetic flux in newly emerged BMRs is in the sunspots they enclose. The quantity $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ increases gradually from about 0.7 at A_BMR,S = 1 μHem⁻¹ to 0.9 at 6000 μHem⁻¹ (left axis, Fig. 2b). Schrijver (1987) established an empirical relationship between the proportion of the magnetic flux in active regions contained in their sunspots and the proportion of the area of active regions occupied by their sunspots. Assuming it applies to newly emerged BMRs, then

$$\begin{array}{c}\hfill \frac{{\mathrm{\Phi }}_{\mathrm{BMR},\mathrm{S}}}{{\mathrm{\Phi }}_{\mathrm{BMR}}}=5.1\frac{A_{\mathrm{BMR},\mathrm{S}}}{A_{\mathrm{BMR}}},\end{array}$$(9)

Fig. 2. Relationship between BMR magnetic flux and sunspot area. (a) Amount of magnetic flux in newly emerged BMRs, ΦBMR, as a function of the area of the enclosed sunspots, ABMR,S, as indicated by the empirical relationship between ΦBMR and ABMR,S presented in this study (black solid line, Eq. (8)) and the constraints on ΦBMR reported by SH (blue, Eq. (10)) and LCC (red, Eq. (13)). For the estimate of ΦBMR based on the SH result, the lower and upper bounds, corresponding to k = 1 and k = 2, are given. The dashed line follows the amount of magnetic flux in newly emerged BMRs contained in the embedded sunspots, ΦBMR,S (Eq. (7)). (b) Corresponding plot of $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$. To an approximation, $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ is directly proportional to the fraction of the area of each BMR occupied by their sunspots, $\frac{A_{\text{BMR,S}}}{A_{\text{BMR}}}$, indicated on the right axis (Eq. (9)) (see Sects. 2.2 and 2.3 for details).

meaning the trend in $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ with A_BMR,S (left axis, Fig. 2b) is accompanied by $\frac{A_{\text{BMR,S}}}{A_{\text{BMR}}}$ rising in proportion from about 0.14 to 0.18 (right axis). In other words, the smallest newly emerged BMRs are about 70% sunspots by magnetic flux and 14% sunspots by area, and this increases steadily to about 90% and 18%, respectively for the largest newly emerged BMRs.

We note the following. Firstly, active regions and newly emerged BMRs are not equivalent to one another, meaning the application of the findings of Schrijver (1987) to the latter (Eq. (9)) is not strictly correct. Secondly, some of the smallest BMRs might contain pores instead of fully developed sunspots (e.g. Švanda et al. 2016), such that the amount of sunspot magnetic flux they enclose might be smaller than indicated by Eq. (7). Nonetheless, the assertion here that the sunspots in a given newly emerged BMR can enclose most of its magnetic flux while occupying only a small proportion of its area is reasonable and within expectation given that sunspots and pores are magnetically denser than faculae.

2.3. Comparison to Schrijver & Harvey (1994) and Lemerle et al. (2015)

In this section we compare the empirical relationship between Φ_BMR and A_BMR,S, derived here (Eq. (8)) to the constraints on Φ_BMR reported by SH (Schrijver & Harvey 1994) and by LCC (Lemerle et al. 2015). Taking Kitt Peak magnetograms, SH isolated the faculae in active regions and newly emerged BMRs and found the facular magnetic flux in individual regions to scale linearly with the area extended by these features. Based on this result and assuming an average BMR magnetic flux density of 100 G or 1.52 × 10¹⁰ Wb μHem⁻¹, SH estimated that

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BMR}}=1.52\times {10}^{10}k{A}_{\mathrm{BMR}}\mathrm{Wb},\end{array}$$(10)

where k is a free parameter that the authors did not constrain. SH interpreted the quantity k − 1 for a given newly emerged BMR as the sum of Φ_BMR,S, and the magnetic flux that dissipated before the BMR had fully formed as a proportion of Φ_BMR. Since this ratio is unlikely to exceed unity, at least not by a wide margin, k most probably lies between one and two. Instead of establishing the relationship between Φ_BMR and A_BMR directly, SH had estimated it from the relationship between facular magnetic flux and area. This was necessary because, as noted in precursor studies to the SH study conducted by the same authors (Schrijver 1987; Harvey & Zwaan 1993), sunspot magnetic flux cannot be taken directly, at least not reliably, from Kitt Peak magnetograms.

Ideally, to relate Φ_BMR to A_BMR,S using the relationship between Φ_BMR and A_BMR reported by SH (Eq. (10)), it should be combined with the relationship between A_BMR and A_BMR,S. Unfortunately, the relationship between A_BMR and A_BMR,S is not known. Instead, let us assume the empirical relationship between the total surface coverage by faculae and by sunspots reported by Chapman et al. (1997) (hereafter CCD) can be taken as the relationship between the area of the faculae and of the sunspots in individual newly emerged BMRs, such that

$$\begin{array}{c}\hfill A_{\mathrm{BMR}}=414+22A_{\mathrm{BMR},\mathrm{S}}-0.0036A_{\mathrm{BMR},\mathrm{S}}^2.\end{array}$$(11)

Inserting Eq. (11) into Eq. (10), we obtain

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BMR}}=1.52\times {10}^{10}k[414+22A_{\mathrm{BMR},\mathrm{S}}-0.0036A_{\mathrm{BMR},\mathrm{S}}^2]\mathrm{Wb}.\end{array}$$(12)

The estimate of Φ_BMR from this expression at k = 1 and k = 2 represents the lower and upper bound of this quantity based on the combination of the SH and CCD findings.

Lemerle et al. (2015) compared the Φ_BMR measurements in the database of newly emerged BMRs compiled by Wang & Sheeley (1989) to USAF-NOAA sunspot areas. They could not match each newly emerged BMR in the Wang & Sheeley (1989) database to the corresponding sunspot group in the USAF-NOAA record, which excluded a direct comparison and compelled them to compare the distribution of values in the two records instead. From this, LCC estimated that

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{BMR}}=1.12\times {10}^{12}{A}_{\mathrm{BMR},\mathrm{S}}^{0.714}\mathrm{Wb}.\end{array}$$(13)

Taking Eqs. (12) and (13), we recalculate Φ_BMR and $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$, drawn in Fig. 2 (blue and red lines) along the values from the empirical relationship between Φ_BMR and A_BMR,S derived in this study (black solid lines).

Figure 2a shows, at lower levels of A_BMR,S, that the estimates of Φ_BMR based on the SH (blue) and LCC results (red) are higher than the current estimate (black solid line). However, Φ_BMR increases with increasing A_BMR,S at a slower rate in the former, such that at higher levels of A_BMR,S the estimates of Φ_BMR based on the SH and LCC results drop below the current estimate and also below Φ_BMR,S (dashed line). This translates into the analogous behaviour of $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ illustrated in Fig. 2b. At low A_BMR,S the estimates of $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ based on the SH (blue) and LCC results (red) are lower than the current estimate (black), but they rise more rapidly with increasing A_BMR,S such that they eventually exceed not just the current estimate, but also unity. The amount of magnetic flux in a given BMR has to be greater than the amount carried by the sunspots it encloses, meaning Φ_BMR going below Φ_BMR,S or, equivalently, $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ going above unity is unphysical. The observation here that the estimates of Φ_BMR based on the SH and LCC results drop below Φ_BMR,S towards higher A_BMR,S is unlikely to be from any uncertainty in the latter. The margin is so wide that even if the mean magnetic flux density over sunspots is only one-quarter the value 1350 G we adopt here, Φ_BMR,S will still be greater than the estimates of Φ_BMR based on the SH and LCC results for the largest BMRs. The constraint on Φ_BMR reported by SH, at least when combined with the findings of CCD as we do here (Eq. (11)), and that by LCC appear to grossly underestimate the amount of magnetic flux in larger BMRs.

The discrepancies between the results reported in this study and earlier by SH and LCC that we have just highlighted could have arisen from the following. SH examined active regions and newly emerged BMRs together even though the two are not equivalent. Furthermore, SH had to constrain Φ_BMR indirectly as sunspot magnetic flux cannot be determined directly and reliably from the magnetogram observations they utilised. How the SH results compare to the LCC and current results could have also been affected by the fact that we had to combine the SH result with the findings of CCD in order to make this comparison, and the latter does not strictly apply to newly emerged BMRs.

As for the LCC study, the authors could not compare the Wang & Sheeley (1989) and USAF-NOAA records directly because they could not match each newly emerged BMR in the former to the corresponding sunspot group in the latter. The Wang & Sheeley (1989) database includes spot-free BMRs, which are not represented in the USAF-NOAA record and should have been excluded. Wang & Sheeley (1989) had also determined Φ_BMR approximately from photographic magnetograms, introducing further uncertainty to the LCC analysis. SH and LCC were hampered by limitations in the data available to them. As stated in the introduction, it would be straightforward to derive the relationship between Φ_BMR and A_BMR,S from an extensive and reliable database of the magnetic and physical properties of newly emerged BMRs, but no such record currently exists. The Wang & Sheeley (1989) database is clearly not ideal for this purpose. The BMR databases compiled more recently by Yeates et al. (2007) and by Yeates (2020), although of excellent quality, are not suitable either as they tabulate the physical properties of each BMR at the time it crosses the central meridian, not when it newly emerges.

The YSK model and the analysis based on it presented in this study have their own assumptions and limitations, discussed in Sects. 2.1 and 2.2. Nonetheless, it is worth noting that the YSK model allowed us to circumvent the lack of any proper database of newly emerged BMRs, an issue that severely limited the SH and LCC analyses, and constrains the relationship between Φ_BMR and A_BMR,S from ΣΦ_FN and ΣA_S (i.e. disc-integrated observations) instead.

3. Application to a surface flux transport model

In this section we look at the result of applying the empirical relationship between Φ_BMR and A_BMR,S presented in this study (Eq. (8)) and the constraints on Φ_BMR reported by SH (Eq. (10)) and by LCC (Eq. (13)) to the SFTM by CJSS. We generated three surface flux transport simulation runs, denoted CJSS-SH, CJSS-LCC, and CJSS-YEA (Table 1). The CJSS-SH run is generated with the CJSS model as it is and set up following the ‘reference case’ presented by Cameron et al. (2010). As we explain below, the CJSS model (incidentally and not by design) already complies with the constraint on Φ_BMR reported by SH. The CJSS-LCC and CJSS-YEA runs are generated with the CJSS model modified to comply with the constraints on Φ_BMR reported by LCC and derived here instead. In the following, we describe the CJSS model and the simulation set-up (Sect. 3.1) before discussing the results (Sect. 3.2).

3.1. Simulation set-up

The governing equation of the CJSS model is

$$\begin{array}{cc}\hfill \frac{\partial B}{\partial t}=& -\mathrm{\Omega }\left(\lambda \right)\frac{\partial B}{\partial \varphi }-\frac{1}{R_{\odot }cos\lambda }\frac{\partial }{\partial \lambda }[\upsilon \left(\lambda \right)Bcos\lambda ]\hfill \\ \hfill & +{\eta }_{\mathrm{S}}[\frac{1}{R_{\odot }^{2}cos\lambda }\frac{\partial }{\partial \lambda }(cos\lambda \frac{\partial B}{\partial \lambda })+\frac{1}{R_{\odot }^2{cos}^2\lambda }\frac{{\partial }^{2}B}{\partial {\varphi }^2}]\hfill \\ \hfill & +D\left({\eta }_{\mathrm{R}}\right)+S(\lambda ,\varphi ,t),\hfill \end{array}$$(14)

where R_⊙, λ, and ϕ denote the radius of the Sun, latitude, and longitude, respectively. This describes the passive transport of the radial component of the photospheric magnetic field, B, under the influence of differential rotation, Ω; meridional flow, υ; and diffusion, which is characterised by the surface and radial diffusivity, η_S and η_R (DeVore et al. 1984; Wang et al. 1989; Mackay et al. 2000; Schrijver et al. 2002; Baumann et al. 2004). Radial diffusion is introduced through a linear operator, D (Baumann et al. 2006), and BMR emergence through a source term, S. As noted in the previous paragraph, we set up the simulation runs following the reference case presented by Cameron et al. (2010). We refer the reader to this earlier paper for details on the definition of the initial condition and the Ω, υ, η_S, and η_R parameters. For the CJSS-SH run, we set up the source term following Cameron et al. (2010). For the CJSS-LCC and CJSS-YEA runs, we modified the source term in such a way as to conform the model to the constraints on Φ_BMR reported by LCC (Eq. (13)) and derived here (Eq. (8)). We describe the definition of the source term and these modifications next.

The source term (i.e. details of BMR emergence) is inferred from the composite of the RGO and USAF-NOAA sunspot area records by Hathaway & Choudhary (2008). The location and area of each sunspot group when the latter is at its maximum is taken to denote the location and A_BMR,S of a newly emerged BMR. Following van Ballegooijen et al. (1998), each BMR is described as two circular magnetic patches of opposite polarity. Let B₊(λ,ϕ) represent the magnetic flux density of the positive polarity patch at (λ,ϕ), (λ₊,ϕ₊) the centre of the patch, and β₊(λ,ϕ) the heliocentric angle between (λ,ϕ) and (λ₊,ϕ₊), and let B₋(λ,ϕ), (λ₋,ϕ₋) and β₋(λ,ϕ) represent the same for the negative polarity patch. The magnetic flux density profile is given by

$$\begin{array}{c}\hfill B_{\pm }(\lambda ,\varphi )=B_{\max }{\left(\frac{0.4\mathrm{\Delta }\beta }{\delta }\right)}^{2}exp(-\frac{2[1-cos{\beta }_{\pm }(\lambda ,\varphi )]}{{\delta }^2}),\end{array}$$(15)

where Δβ = β₊(λ₋,ϕ₋) = β₋(λ₊,ϕ₊) denotes the heliocentric angle between the two patches and δ = 4° the size of each patch. The constant factor B_max modulates the absolute scale of the model output and is fixed at 374 G by matching the output average unsigned magnetic flux density to the observed values from the Mount Wilson Observatory (MWO) and Wilcox Solar Observatory (WSO) (Arge et al. 2002), illustrated in Fig. 3a. We note that BMR magnetic flux density, B_±, and therefore the amount of magnetic flux in each BMR, Φ_BMR, are assumed to scale linearly with the square of the angular separation, (Δβ)². The angular separation, Δβ, is given by

$$\begin{array}{c}\hfill \mathrm{\Delta }\beta =6.46\times {10}^{-2}\sqrt{A_{\mathrm{BMR}}},\end{array}$$(16)

Fig. 3. Surface flux transport simulation runs. (a) Average unsigned photospheric magnetic flux density in the CJSS-SH (blue), CJSS-LCC (red), and CJSS-YEA (black) simulation runs (Table 1). The CJSS-SH and CJSS-LCC time series are largely hidden by the CJSS-YEA time series due to the close similarity. The circles and squares denote the observed values from MWO and WSO (Arge et al. 2002). (b) Open flux corresponding to each run, derived by extrapolating the model output photospheric magnetic field into the heliosphere, and the reconstruction of this quantity from the aa index by Lockwood et al. (2009) (green solid line) and by Lockwood et al. (2014) (green dashed line). In (a) and (b) the CJSS-LCC and CJSS-YEA time series are dashed and dotted, respectively, to aid visibility. (c) Mean magnetic field over the north (solid lines) and the south polar caps (dashed lines) in each run. The green crosses indicate the observed timing of the polar field reversals reported by Makarov et al. (2004) (see Sects. 3.1 and 3.2 for details).

where Δβ and A_BMR are in units of degrees and μHem, respectively, and A_BMR is calculated from A_BMR,S using the CCD relationship (Eq. (11)). Cameron et al. (2010) established this empirical relationship with the Δβ and sunspot umbral area measurements from the Kodaikanal Solar Observatory (KSO) and MWO. The tilt angle of each BMR is calculated from the observed cycle-average tilt angle, also from KSO and MWO, through the procedure detailed in Cameron et al. (2010). The source term, and therefore the simulation runs, extend from 1913 to 1986 (solar cycle 15 to 21), the period over which sunspot group location, area, and tilt angle measurements are available. The simulation runs and the analysis based on them (see the following sections) are not sensitive to which sunspot area composite we take into the CJSS model as input. The Hathaway & Choudhary (2008), Balmaceda et al. (2009), and Mandal et al. (2020) composites are essentially different efforts to extend the RGO record, which ends in 1976, to more recent times with sunspot area measurements from other observatories. The various composites are therefore exactly identical to one another up to this time (i.e. over most of the simulation period of 1913 to 1986).

The CJSS model relies on the combination of the assumption that B_± scales linearly with (Δβ)² (Eq. (15)) and the empirical relationship between Δβ and A_BMR (Eq. (16)) to define B_±. Substituting Eq. (16) into Eq. (15) we obtain

$$\begin{array}{c}\hfill B_{\pm }(\lambda ,\varphi )={\stackrel{\sim }{B}}_{\max }{A}_{\mathrm{BMR}}exp(-\frac{2[1-cos{\beta }_{\pm }(\lambda ,\varphi )]}{{\delta }^2}),\end{array}$$(17)

where ${\stackrel{\sim }{B}}_{\text{max}}=B_{\text{max}}0.4^2{(6.46\times {10}^{-2})}^2{\delta }^{-2}=0.0156$ G μHem⁻¹ is a constant factor that serves the same purpose as B_max. The consequence of this approach is that B_±, and therefore also Φ_BMR, scales linearly with A_BMR in the model. Incidentally, this is similar to the constraint on Φ_BMR reported by SH, where this quantity also scales linearly with A_BMR (Eq. (10)).

Let us generalise Eq. (17) by rewriting it as

$$\begin{array}{c}\hfill B_{\pm }(\lambda ,\varphi )={\stackrel{\sim }{B}}_{\max }fexp(-\frac{2[1-cos{\beta }_{\pm }(\lambda ,\varphi )]}{{\delta }^2}),\end{array}$$(18)

where f = A_BMR determines the A_BMR,S-dependence of B_± and Φ_BMR. In the CJSS-LCC run we modify the CJSS model to comply with the constraint on Φ_BMR reported by LCC (Eq. (13)) by redefining f as $A_{\text{BMR},\text{S}}^{0.71}$. Similarly, in the CJSS-YEA run we apply the constraint on Φ_BMR derived in this study (Eq. (8)) to the CJSS model by redefining f as ${\stackrel{\sim }{A}}_{\text{BMR}}$. By this redefinition of f, we replace the assumption in the CJSS model that B_± scales linearly with (Δβ)² with direct observational constraints on the relationship between Φ_BMR and A_BMR,S. As in the CJSS-SH run, ${\stackrel{\sim }{B}}_{\text{max}}$ is set in the CJSS-LCC and CJSS-YEA runs at the value that brings the modelled average unsigned magnetic flux density to the absolute scale of the observed values (Fig. 3a). The value of ${\stackrel{\sim }{B}}_{\text{max}}$ and the definition of f in each simulation run is summarised in Table 1.

3.2. Simulation results

3.2.1. Total photospheric magnetic flux

In terms of the total photospheric magnetic flux, represented in Fig. 3a by the average unsigned magnetic flux density, the three simulation runs (black, blue, and red curves) are very similar. The agreement in absolute scale is moot since ${\stackrel{\sim }{B}}_{\text{max}}$ is set in each run at the level that brings the model output to the absolute scale of the observed values (circles and boxes). What is notable is that the time variability is also very similar even with the different A_BMR,S-dependence of Φ_BMR imposed for each run.

3.2.2. Solar open flux and polar field reversals

Cameron et al. (2010) and Jiang et al. (2011b) had verified their model by demonstrating its ability to reproduce the evolution of the solar open flux and the timing of the polar field reversals. The authors derived the open flux corresponding to the model output photospheric magnetic field by extrapolating it into the heliosphere using the Current Sheet Source Surface (CSSS) model by Zhao & Hoeksema (1995a,b). This was matched to the independent reconstruction of the open flux from the aa geomagnetic activity index by Lockwood et al. (2009). They compared the timing of the polar field reversals in the model output to the timings reported by Makarov et al. (2004) based on an examination of the polar filaments in Hα synoptic charts.

Following CJSS, we applied the CSSS model to each simulation run to calculate the corresponding solar open flux. In Fig. 3b we compare the results (black, blue, and red) to the Lockwood et al. (2009) reconstruction (green solid line) and the update by Lockwood et al. (2014) (green dashed line). The various runs do not differ from one another more than they do from the Lockwood et al. (2009, 2014) time series.

In Fig. 3c we chart the evolution of the polar field in the simulation runs as given by the mean magnetic flux density over the north (λ > 75°, solid lines) and the south polar caps (λ < −75°, dashed lines). Polar field reversal is indicated by the mean polar field switching polarity. In the same plot we indicate the observed timing of the polar field reversals reported by Makarov et al. (2004) (green crosses). Analogous to the solar open flux, the various runs do not differ more from one another than they do from the Makarov et al. (2004) observations. Applying the constraints on Φ_BMR presented by LCC (Eq. (13)) and derived in this study (Eq. (8)) to the CJSS model retained its ability to reproduce the evolution of the solar open flux and the timing of the polar field reversals.

3.2.3. Correlation between the polar field and solar cycle strength

A key finding of the CJSS studies is that the amount of polar flux at the end of each solar cycle in the model appears to be highly correlated to the observed strength of the succeeding cycle. This is consistent with the existence of a Babcock–Leighton-type solar dynamo (Babcock 1961; Leighton 1964, 1969), where the polar field, which peaks around the end of each cycle, is the main source of the toroidal field in the convection zone (see discussion in Cameron & Schüssler 2015). We examine if this feature of the CJSS model is sensitive to the A_BMR,S-dependence of Φ_BMR assumed. Taking each simulation run, we compare the peak of the mean unsigned magnetic flux density over both polar caps at the end of each cycle to the maximum total sunspot area over the following cycle, representing modelled polar flux and measured cycle strength (Fig. 4). The Pearson’s correlation coefficient, R, between the two quantities is 0.933 for the CJSS-SH run, 0.925 for the CJSS-LCC run, and 0.937 for the CJSS-YEA run. The different A_BMR,S-dependence of Φ_BMR imposed on the various runs made no appreciable difference to the correlation between the model output polar field at the end of each cycle and the observed strength of the succeeding cycle.

Fig. 4. For each simulation run, the peak of the mean unsigned magnetic flux density over the polar caps at the end of solar cycle n vs the maximum total surface coverage by sunspots over cycle n + 1, where n ranges from 15 to 21. We excluded solar cycle 14 as the end of this cycle occurred at the start of the simulation, and is therefore sensitive to the arbitrary initial condition. The crosses and lines denote the individual cycle values and the corresponding linear fit. The Pearson’s correlation coefficient, R, is indicated (see Sect. 3.2.3 for details).

Whether in terms of the time variation of the total photospheric magnetic flux (Fig. 3a) and the solar open flux (Fig. 3b), the timing of the polar field reversals (Fig. 3c), or the correlation between the polar field at the end of each cycle and the observed strength of the succeeding cycle (Fig. 4), the various simulation runs are very similar. It appears that the CJSS model is, at least to the extent tested here, not particularly sensitive to the A_BMR,S-dependence of Φ_BMR assumed. This is likely, at least in part, due to the fact that the magnetic flux in each newly emerged BMR can persist for weeks or months, such that the state of a simulation run at a given time depends not only on the BMRs injected then but also on all the BMRs that were introduced in the weeks to months prior. The amount of magnetic flux in individual BMRs depends strongly on the A_BMR,S-dependence of Φ_BMR assumed. However, the aggregate magnetic flux introduced over an extended period by all the BMRs that had emerged then would not be as sensitive to the A_BMR,S-dependence of Φ_BMR assumed due to the summation. Nonetheless, it is noteworthy that the application of the empirical relationship between Φ_BMR and A_BMR,S presented here to what is an established SFTM retained its key features, including its ability to replicate various independent datasets.

4. Balance between bright faculae and network and dark sunspots at various solar activity levels

As discussed in the Introduction, the balance between bright faculae and network and dark sunspots on the solar disc and how that changes with solar activity level is of interest for the relevance to the understanding of the brightness variability of the Sun and other cool stars. Making use of the CJSS-YEA run (Sect. 3), let us examine how the amount of magnetic flux on the solar disc in faculae and network, ΣΦ_FN, and in sunspots, ΣΦ_S, relate to one another and how the ratio $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ changes with solar activity level. The quantity ΣΦ_S is calculated from the Balmaceda et al. (2009)ΣA_S time series (Fig. 1a) using Eq. (7). We obtained the total magnetic flux on the solar disc from the CJSS-YEA run and subtracted ΣΦ_S from it to yield ΣΦ_FN. In the following discussion we examine the annual mean of each quantity.

In Fig. 5a we chart ΣΦ_FN estimated from the CJSS-YEA run as a function of ΣΦ_S (black crosses). As a check, we compare this to the YSK model reconstruction of ΣΦ_FN from ΣΦ_S (red crosses). Before we discuss the trend in ΣΦ_FN with ΣΦ_S, let us address the following. Firstly, both the CJSS and YSK models describe the evolution of the photospheric magnetic field due to sunspot-bearing BMRs alone, ignoring ERs. Secondly, the two estimates of ΣΦ_FN differ in absolute scale and we have to rescale the values from the YSK model by a factor of 1.46 to bring them to the absolute scale of the values from the CJSS-YEA run. Thirdly, plotted against ΣΦ_S, the values from the CJSS-YEA run are markedly more scattered than the values from the YSK model.

Fig. 5. Disc-integrated facular and network magnetic flux, ΣΦFN, as indicated by the CJSS-YEA run (black crosses) and by the YSK model (red crosses), as a function of the disc-integrated sunspot magnetic flux, ΣΦS. The first solar cycle in the CJSS-YEA run (cycle 15, 1913 to 1923) has been excluded as the simulation is sensitive to the arbitrary initial condition here. The annual mean over the period of 1924 to 1986, where there is data from both models, is shown. The black line and the grey region represents the quadratic polynomial fit to the CJSS-YEA run values and the corresponding 95% confidence interval, while the red line represents the fit to the YSK model values. See Sect. 4 for details.

Even though the CJSS model neglects ERs, considering the ability of the CJSS-YEA run to reproduce various independent datasets (Sect. 3.2), the effect of this limitation on the corresponding estimate of ΣΦ_FN is likely minor. A similar argument applies to the YSK model estimate of ΣΦ_FN. The absolute scale of the CJSS model is set to that of the average unsigned magnetic flux density values reported by Arge et al. (2002), which is based on MWO and WSO magnetograms (circles and squares, Fig. 4a). The absolute scale of the YSK model is set to that of the ΣΦ_FN time series reported by Yeo et al. (2014b), which is based on Kitt Peak, MDI, and HMI magnetograms (black, Fig. 1b). We attribute the different absolute scale of the CJSS and YSK model outputs to systematic uncertainties in how Arge et al. (2002) and Yeo et al. (2014b) derived their respective datasets and in the absolute magnetogram signal scale of the various solar telescopes (Riley et al. 2014). In the YSK model, the variation in ΣΦ_FN due to sunspot-bearing BMRs is given by the convolution of $\mathrm{\Sigma }{A}_{\text{S}}^{h_1}$ with a FIR filter, H (Eqs. (1) and (2)). The implicit assumption is that all BMRs exhibit the same time evolution as delineated by H. This ignores the fact that even BMRs containing the same amount of magnetic flux can evolve rather differently due to factors such as latitude, tilt angle, and the magnetic field in their surroundings. These factors are taken into account in the CJSS model, producing the greater scatter in the corresponding ΣΦ_FN versus ΣΦ_S plot.

To more clearly show the trend in ΣΦ_FN with ΣΦ_S, we fit a quadratic polynomial in ΣΦ_S to each of the two estimates of ΣΦ_FN. The fit to the CJSS-YEA run values (black line, Fig. 5) is given by

$$\begin{array}{cc}\hfill \mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}=& (0.52\pm 0.11)+(3.33\pm 0.53)\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}\hfill \\ \hfill & -(0.73\pm 0.46)\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}^2\hfill \end{array}$$(19)

and the fit to the YSK model values (red line) by

$$\begin{array}{cc}\hfill \mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}=& (0.334\pm 0.019)+(3.367\pm 0.086)\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}\hfill \\ \hfill & -(0.425\pm 0.074)\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}^2.\hfill \end{array}$$(20)

This revealed that ΣΦ_FN increases with increasing ΣΦ_S at a diminishing rate such that it appears to gradually saturate, analogous to what earlier studies noted comparing facular area and various disc-integrated quantities that are sensitive to the amount of faculae and network present to sunspot area and number (Foukal 1993, 1996, 1998; Chapman et al. 1997; Solanki & Fligge 1998; Shapiro et al. 2014; Yeo et al. 2020), discussed in the Introduction. The fit to the CJSS-YEA run values (black line) and the fit to the YSK model values (red line) do differ slightly, though the latter lies within the 95% confidence interval of the former (shaded region). This means that, taking the scatter in the CJSS-YEA run values (black crosses) into account, that there is little separating the two estimates of ΣΦ_FN in terms of the ΣΦ_S-dependence. The CJSS-YEA run and the YSK model are not entirely independent of one another as we had constrained the former using the empirical relationship between Φ_BMR and A_BMR,S derived from the latter. Nonetheless, it is significant that two models that adopt vastly different approaches to modelling the evolution of the photospheric magnetic field returned similar results.

Next, we examine what the apparent trend in ΣΦ_FN with ΣΦ_S implies for the balance between the amount of faculae and network and the amount of sunspots on the solar disc and how that changes with solar activity level. Sunspot indices are often employed as proxies of solar activity, while stellar activity is often characterised by the quantity $log{R}_{\text{HK}}^{\prime }$ defined in Noyes et al. (1984). Taking this into consideration, let us examine $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$, as estimated from the CJSS-YEA run, as a function of ΣA_S and of $log{R}_{\text{HK}}^{\prime }$ (crosses, Fig. 6). We make use of the model reconstruction of $log{R}_{\text{HK}}^{\prime }$ from ΣA_S presented by YSK.

Fig. 6. Variation in the balance between bright faculae and network and dark sunspots with solar activity. (a) The ratio $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$, where ΣΦFN is obtained from the CJSS-YEA run, as a function of ΣAS (crosses). As in Fig. 5, the annual mean over the period of 1924 to 1986 is depicted. The solid line follows the power-law fit (Eq. (22)) and the dotted line gives the position of the horizontal asymptote as given by the constant term in this fit. (b) The same, but for $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ as a function of $log{R}_{\text{HK}}^{\prime }$. The vertical dashed line indicates $log{R}_{\text{HK}}^{\prime }=-4.9$ and the horizontal dashed lines gives the value of $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ at $log{R}_{\text{HK}}^{\prime }=-4.9$ and $log{R}_{\text{HK}}^{\prime }=-4.7$, as given by the power-law fit (Eq. (23)) (see Sect. 4 for details).

The gradual saturation in ΣΦ_FN as we go from low to high ΣΦ_S (Fig. 5) translates into $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ decreasing with increasing solar activity level at a diminishing rate such that it appears to approach a horizontal asymptote (Fig. 6). Let us model $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ as conforming to a power-law relationship with activity level such that

$$\begin{array}{c}\hfill \frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}}=g_1{[\alpha +g_2]}^{g_3}+g_4,\end{array}$$(21)

where α represents activity as given by either ΣA_S or $log{R}_{\text{HK}}^{\prime }$ and g₁ to g₄ are fit parameters. The g₁[α + g₂]^g₃ term describes the variation in $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ with α, while g₄ denotes the position of the horizontal asymptote. We fit Eq. (21) to the $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ versus ΣA_S plot (Fig. 6a) and the $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ versus $log{R}_{\text{HK}}^{\prime }$ plot (Fig. 6b) in turn. For the former, g₂ is fixed at null. At convergence we arrive at

$$\begin{array}{c}\hfill \frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}}=(690\pm 110)\mathrm{\Sigma }{A}_{\mathrm{S}}^{(-0.876\pm 0.043)}+(2.63\pm 0.29)\end{array}$$(22)

and

$$\begin{array}{cc}\hfill \frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\mathrm{S}}}=& (0.038\pm 0.060){[logR{\prime }_{\mathrm{HK}}+(4.9405\pm 0.0046)]}^{(-1.35\pm 0.45)}\hfill \\ \hfill & +(2.39\pm 0.72),\hfill \end{array}$$(23)

represented in Fig. 6 by the solid lines. In both instances the g₄ term or the position of the horizontal asymptote (dotted lines) is above unity, implying that the amount of magnetic flux on the solar disc in sunspots never exceeds that in faculae and network. This is consistent with Parnell et al. (2009) and Thornton & Parnell (2011), who came to a similar conclusion by examining the distribution of magnetic features with the amount of flux in each feature.

The observation here that most of the magnetic flux on the solar disc is in faculae and network does not contradict our earlier assertion that most of the magnetic flux in newly emerged BMRs is in the enclosed sunspots (Sect. 2.2). The magnetic flux in BMRs persists long after the comparatively short-lived sunspots they enclose have dissipated, and manifests as faculae and network (van Driel-Gesztelyi & Green 2015). In combination with the fact that the photospheric magnetic field at a given time is comprised of not just the BMRs that had emerged then, but also what remains from all the BMRs that had emerged in the weeks to months prior, the two assertions we are making do not exclude one another.

In the context of the brightness variability of cool stars, the transition between the faculae-dominated and spot-dominated regimes is estimated to be in the range of $log{R}_{\text{HK}}^{\prime }$ of −4.9 to −4.7 (Lockwood et al. 2007; Hall et al. 2009; Shapiro et al. 2014; Radick et al. 2018; Reinhold et al. 2019). Looking at Fig. 6b, the $log{R}_{\text{HK}}^{\prime }$ of the Sun (crosses) is in the range of −4.93 to −4.84, around the lower end of the range of activity where this transition occurs. The extrapolation of the power-law fit (solid line) indicates that in the range of $log{R}_{\text{HK}}^{\prime }$ of −4.9 to −4.7, $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ varies slowly with activity, such that it stays within a relatively narrow range of about 3 to 5 (dashed lines). Taken together, these two observations imply that the $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ level of the Sun might not be far above the level where facular and network brightening and sunspot darkening would cancel one another out. This concurs with the findings of Shapiro et al. (2016) and Radick et al. (2018), who noted that while the Sun is faculae-dominated, it is only marginally so as facular and network brightening and sunspot darkening appear to be closely balanced.

5. Summary

The relationship between sunspot-bearing bipolar magnetic regions (BMRs) and the sunspots they enclose, connected to how solar magnetism and activity relate to sunspot prevalence, is important. It is particularly relevant to surface flux transport models (SFTMs) that derive the details of BMR emergence they require as input from sunspot observations. However, attempts to constrain how the amount of magnetic flux in each newly emerged BMR, Φ_BMR, relates to its area (Schrijver & Harvey 1994) or the area of the enclosed sunspots, A_BMR,S (Lemerle et al. 2015), were limited by the fact that there is no extensive and reliable database of the magnetic and physical properties of newly emerged BMRs currently available. In the absence of proper observational constraints, existing SFTMs adopt various simplifications to determine the details of BMR emergence from sunspot observations. In this study we aimed to establish the relationship between Φ_BMR and A_BMR,S.

To this end, we made use of the empirical model of the relationship between the disc-integrated facular and network magnetic flux, ΣΦ_FN, and the total surface coverage by sunspots, ΣA_S, presented by Yeo et al. (2020). The structure of the Yeo et al. (2020) model enabled us to circumvent the lack of any proper BMR database, and instead derive an empirical relationship between Φ_BMR and A_BMR,S (Eq. (8)) from ΣA_FN and ΣA_S (i.e. disc-integrated quantities). The empirical relationship between Φ_BMR and A_BMR,S indicates that for newly emerged BMRs, the enclosed sunspots occupy a small proportion of their area (14% to 18%), but contain most of their magnetic flux (70% to 90%) (Fig. 2b). We found evidence that the earlier results by Schrijver & Harvey (1994) (Eq. (10)) and Lemerle et al. (2015) (Eq. (13)) might underestimate the Φ_BMR of larger BMRs (Fig. 2a).

To verify the empirical relationship between Φ_BMR and A_BMR,S presented here, we examined the result of applying this and the earlier results by Schrijver & Harvey (1994) and Lemerle et al. (2015) to a SFTM. We made use of the SFTM reported by Cameron et al. (2010) and Jiang et al. (2011a,b). This model, incidentally and not by design, already complies with the Schrijver & Harvey (1994) result. We generated three simulation runs, one with the model as it is and the other two with the model modified to comply with the Lemerle et al. (2015) results and the current results. We compared the evolution of the solar open flux to the independent reconstruction of this quantity from the aa geomagnetic activity index by Lockwood et al. (2009, 2014) (Fig. 3b) and the timing of the polar field reversals to the observed timing reported by Makarov et al. (2004) (Fig. 3c). The various runs reproduce the Lockwood et al. (2009, 2014) and Makarov et al. (2004) results relatively well and almost equally. An important feature of the Cameron et al. (2010) and Jiang et al. (2011a,b) model is that the polar field at the end of each cycle is highly correlated to the observed strength of the following cycle. This is significant as it is consistent with the existence of a Babcock–Leighton-type solar dynamo. In terms of this correlation, the various runs are again very similar (Fig. 4). Surface flux transport simulations do not appear to be, to the extent tested here, particularly sensitive to the A_BMR,S-dependence of Φ_BMR assumed. Nonetheless, it is worth noting that applying the current result to what is an established SFTM retained its key features.

Making use of the surface flux transport simulation based on the relationship between Φ_BMR and A_BMR,S presented here, we examined how the amount of magnetic flux on the solar disc in faculae and network, ΣΦ_FN, and in sunspots, ΣΦ_S, relate to one another and the relationship between the ratio $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ and the solar activity level. We found ΣΦ_FN to increase with increasing ΣΦ_S at a diminishing rate such that it appears to gradually saturate (Fig. 5). This complies with the analogous results from the earlier studies that compared facular area and various disc-integrated quantities sensitive to the amount of faculae and network present to sunspot area and number (Foukal 1993, 1996, 1998; Chapman et al. 1997; Solanki & Fligge 1998; Shapiro et al. 2014; Yeo et al. 2020). Taking ΣA_S and log $log{R}_{\text{HK}}^{\prime }$ as proxies of solar activity, we found $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ to decrease with increasing activity at a diminishing rate such that it appears to approach a horizontal asymptote (Fig. 6). We modelled $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ as conforming to a power-law relationship with activity (Eqs. (22) and (23)). The power-law model indicates that, at least within the range of activity examined, the proportion of the magnetic flux on the solar disc that is in sunspots remains smaller than that in faculae and network. This does not contradict our claim that most of the magnetic flux in newly emerged BMRs is in the enclosed sunspots (see Sect. 4). The extrapolation of the power-law model to higher activity levels indicates that the $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ level of the Sun might be close to the level where facular and network brightening and sunspot darkening would cancel one another out, consistent with Shapiro et al. (2016) and Radick et al. (2018).

Acknowledgments

The authors are grateful to R. H. Cameron of the Max-Planck Institut für Sonnensystemforschung for the useful discussions and for granting access to the surface flux transport model featured in this work. KLY and NAK received funding from the German Federal Ministry of Education and Research (project 01LG1909C), KLY and SKS from the European Research Council through the Horizon 2020 research and innovation program of the European Union (grant 695075) and JJ from the National Science Foundation of China (grant 11873023 and 11522325) and the Fundamental Research Funds for the Central Universities of China.

References

All Tables

All Figures

Fig. 1. Observed ΣAS and ΣΦFN, YSK model (Eq. (1)) reconstruction of ΣΦFN from ΣAS and H term in this model (Eq. (2)). (a) Composite time series of the total surface coverage by sunspots, ΣAS, by Balmaceda et al. (2009). (b) Composite time series of the disc-integrated facular and network magnetic flux, ΣΦFN, by Yeo et al. (2014b) (black), the reconstruction of this quantity from ΣAS with the YSK model (red), and the 365-day running mean of the difference between measurement and model (blue). The Balmaceda et al. (2009) time series goes back to 1874 but the range of the plot is confined to 1974 to 2017, the period where it overlaps with the Yeo et al. (2014b) time series. The gaps in the Yeo et al. (2014b) time series correspond to periods where there no suitable magnetogram observations with which to determine ΣΦFN were available. (c) The H term, which describes the response of ΣΦFN to the appearance of a sunspot per unit sunspot area. The dotted line follows the exponential envelope of H and t⊙ denotes the synodic rotation period of the Sun (see Sect. 2.1 for details).

In the text

Fig. 2. Relationship between BMR magnetic flux and sunspot area. (a) Amount of magnetic flux in newly emerged BMRs, ΦBMR, as a function of the area of the enclosed sunspots, ABMR,S, as indicated by the empirical relationship between ΦBMR and ABMR,S presented in this study (black solid line, Eq. (8)) and the constraints on ΦBMR reported by SH (blue, Eq. (10)) and LCC (red, Eq. (13)). For the estimate of ΦBMR based on the SH result, the lower and upper bounds, corresponding to k = 1 and k = 2, are given. The dashed line follows the amount of magnetic flux in newly emerged BMRs contained in the embedded sunspots, ΦBMR,S (Eq. (7)). (b) Corresponding plot of $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$. To an approximation, $\frac{{\mathrm{\Phi }}_{\text{BMR,S}}}{{\mathrm{\Phi }}_{\text{BMR}}}$ is directly proportional to the fraction of the area of each BMR occupied by their sunspots, $\frac{A_{\text{BMR,S}}}{A_{\text{BMR}}}$, indicated on the right axis (Eq. (9)) (see Sects. 2.2 and 2.3 for details).

In the text

Fig. 3. Surface flux transport simulation runs. (a) Average unsigned photospheric magnetic flux density in the CJSS-SH (blue), CJSS-LCC (red), and CJSS-YEA (black) simulation runs (Table 1). The CJSS-SH and CJSS-LCC time series are largely hidden by the CJSS-YEA time series due to the close similarity. The circles and squares denote the observed values from MWO and WSO (Arge et al. 2002). (b) Open flux corresponding to each run, derived by extrapolating the model output photospheric magnetic field into the heliosphere, and the reconstruction of this quantity from the aa index by Lockwood et al. (2009) (green solid line) and by Lockwood et al. (2014) (green dashed line). In (a) and (b) the CJSS-LCC and CJSS-YEA time series are dashed and dotted, respectively, to aid visibility. (c) Mean magnetic field over the north (solid lines) and the south polar caps (dashed lines) in each run. The green crosses indicate the observed timing of the polar field reversals reported by Makarov et al. (2004) (see Sects. 3.1 and 3.2 for details).

In the text

Fig. 4. For each simulation run, the peak of the mean unsigned magnetic flux density over the polar caps at the end of solar cycle n vs the maximum total surface coverage by sunspots over cycle n + 1, where n ranges from 15 to 21. We excluded solar cycle 14 as the end of this cycle occurred at the start of the simulation, and is therefore sensitive to the arbitrary initial condition. The crosses and lines denote the individual cycle values and the corresponding linear fit. The Pearson’s correlation coefficient, R, is indicated (see Sect. 3.2.3 for details).

In the text

Fig. 5. Disc-integrated facular and network magnetic flux, ΣΦFN, as indicated by the CJSS-YEA run (black crosses) and by the YSK model (red crosses), as a function of the disc-integrated sunspot magnetic flux, ΣΦS. The first solar cycle in the CJSS-YEA run (cycle 15, 1913 to 1923) has been excluded as the simulation is sensitive to the arbitrary initial condition here. The annual mean over the period of 1924 to 1986, where there is data from both models, is shown. The black line and the grey region represents the quadratic polynomial fit to the CJSS-YEA run values and the corresponding 95% confidence interval, while the red line represents the fit to the YSK model values. See Sect. 4 for details.

In the text

Fig. 6. Variation in the balance between bright faculae and network and dark sunspots with solar activity. (a) The ratio $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$, where ΣΦFN is obtained from the CJSS-YEA run, as a function of ΣAS (crosses). As in Fig. 5, the annual mean over the period of 1924 to 1986 is depicted. The solid line follows the power-law fit (Eq. (22)) and the dotted line gives the position of the horizontal asymptote as given by the constant term in this fit. (b) The same, but for $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ as a function of $log{R}_{\text{HK}}^{\prime }$. The vertical dashed line indicates $log{R}_{\text{HK}}^{\prime }=-4.9$ and the horizontal dashed lines gives the value of $\frac{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{FN}}}{\mathrm{\Sigma }{\mathrm{\Phi }}_{\text{S}}}$ at $log{R}_{\text{HK}}^{\prime }=-4.9$ and $log{R}_{\text{HK}}^{\prime }=-4.7$, as given by the power-law fit (Eq. (23)) (see Sect. 4 for details).

In the text