How Good Is the Bipolar Approximation of Active Regions for Surface Flux Transport?

Abstract

We investigate how representing active regions with bipolar magnetic regions (BMRs) affects the end-of-cycle polar field predicted by the surface flux transport model. Our study is based on a new database of BMRs derived from the SDO/HMI active region patch data between 2010 and 2020. An automated code is developed for fitting each active region patch with a BMR, matching both the magnetic flux and axial dipole moment of the region and removing repeat observations of the same region. By comparing the predicted evolution of each of the 1090 BMRs with the predicted evolution of their original active region patches, we show that the bipolar approximation leads to a 24% overestimate of the net axial dipole moment, given the same flow parameters. This is caused by neglecting the more complex multipolar and/or asymmetric magnetic structures of many of the real active regions, and may explain why previous flux transport models had to reduce BMR tilt angles to obtain realistic polar fields. Our BMR database and the Python code to extract it are freely available.

1 Introduction

Despite its simplicity, the surface flux transport (SFT) model introduced by Leighton (1964) has proved remarkably effective at mimicking the evolution of the large-scale magnetic field on the solar surface (see the reviews by Sheeley, 2005; Jiang et al., 2014; Wang, 2017). An important recent application is solar cycle prediction. Since the polar field at Solar Minimum is the best predictor for the following solar cycle amplitude (Schatten et al., 1978; Muñoz-Jaramillo et al., 2013; Pesnell, 2016; Petrovay, 2020), the SFT model offers the possibility of extending predictions earlier since it – in turn – predicts this polar field. Several predictions of Cycle 25 have been made with this methodology (Iijima et al., 2017; Jiang et al., 2018; Upton and Hathaway, 2018), as well as by coupling SFT to interior dynamo models (Bhowmik and Nandy, 2018; Labonville, Charbonneau, and Lemerle, 2019).

A major component of randomness in the solar cycle is now believed to arise from active region emergence. There are significant fluctuations both in the number of emerging active regions and in their properties such as emergence latitude, magnetic flux, or tilt angle (see the review by van Driel-Gesztelyi and Green, 2015). General cycle-to-cycle trends remain unclear, with some studies suggesting active regions in stronger solar cycles to have lower tilt angles (e.g., McClintock and Norton, 2013) and others finding no significant variation (Tlatova et al., 2018). Such a trend would act to stabilise the solar cycle by reducing the polar field strength at the end of strong cycles. Another stabilising effect is the tendency for active regions in stronger cycles to emerge at higher latitudes, reducing the cross-equatorial flux transport needed to build the polar field (Jiang et al., 2014). On the other hand, it is now believed that individual “rogue” active regions with statistically extreme properties can have a significant effect on the polar field (Cameron et al., 2013; Jiang, Cameron, and Schüssler, 2014; Nagy et al., 2017; Whitbread, Yeates, and Muñoz-Jaramillo, 2018).

To account for these fluctuations within particular solar cycles, SFT models need to be driven with input data tailored to real observations of the Sun, rather than generic statistical distributions. The most direct method is to assimilate observed magnetograms directly from observed regions of the Sun (e.g., Worden and Harvey, 2000; Upton and Hathaway, 2014; Hickmann et al., 2015). However, for many applications it is useful to isolate individual active region sources. Reasons to do this might include: efficiency and speed of calculation, particularly when running ensembles for prediction purposes (Petrovay, Nagy, and Yeates, 2020); coupling to three-dimensional models for either the solar corona (Mackay and Yeates, 2012; Yeates, 2014) or interior (Bhowmik and Nandy, 2018; Labonville, Charbonneau, and Lemerle, 2019); determining the contribution of individual active regions (Yeates, Baker, and van Driel-Gesztelyi, 2015; Whitbread, Yeates, and Muñoz-Jaramillo, 2018); and extending simulations to previous epochs where sunspot observations are available but magnetograms are not (Virtanen et al., 2019a). Whilst some authors have adopted a compromise approach of assimilating magnetograms for individual regions (Yeates, Baker, and van Driel-Gesztelyi, 2015; Whitbread et al., 2017; Virtanen et al., 2019a), the majority of models fit each source region with a simple bipolar magnetic region (BMR). Public catalogues of magnetic BMRs determined from NSO data are available for Cycle 21 from Wang and Sheeley (1989) and Cycles 23–24 from Yeates (2016).

This paper concerns the consequences of fitting active regions with (symmetric) BMR shapes, when these BMRs are tailored to individual observed active regions. A difficulty with this approach has been found when the magnetic flux, size and tilt angle of each BMR are all matched to the corresponding observed active region. Using the observed values for all of these properties tends to lead to an over-strong polar field after the ensuing SFT evolution (e.g., Yeates, 2014; Cameron et al., 2010; Jiang, Cameron, and Schüssler, 2014; Bhowmik and Nandy, 2018). To correct the discrepancy, all of these authors reduce the tilt angles of the inserted BMRs compared to their observed values, with a typical reduction of 20% to 30%. This reduction is argued to compensate for the lack of inflows towards active regions in the SFT model, which could have the effect of reducing the tilt angle during subsequent evolution of the real region (Cameron et al., 2010). The BMR-driven SFT model for Cycle 24 described in this paper also overestimates polar field. However, by comparing with the same model driven by non-BMR sources, we will show that this overestimate results simply from the assumption of bipolar shape, rather than any systematic error in the transport per se.

A hint of the limitations of bipolar shape comes from the recent work of Iijima, Hotta, and Imada (2019) on the asymmetry of bipolar regions. These authors explore the consequences of the fact that, in real active regions, leading polarity flux tends to be more spatially concentrated than following polarity flux (Fisher et al., 2000). Iijima, Hotta, and Imada (2019) used an SFT model to show that if the following polarity is more dispersed than the leading polarity, then sufficient following flux can escape into the opposite hemisphere to reduce the predicted polar field. Here, we test what is lost in the BMR approximation more generally, using a new automated BMR database for Cycle 24 (Section 2). The effect of fitting BMRs on the SFT evolution is discussed in Section 3 and concluding remarks are made in Section 4.

2 Bipole Database

Our database is derived from the Spaceweather HMI Active Region Patch (SHARP) data from the Helioseismic and Magnetic Imager on Solar Dynamics Observatory (Bobra et al., 2014). Each SHARP represents an individual active region tracked over a single disk passage. Using SHARP data, we avoid the need to repeat the identification of individual active regions (except for regions with more than one disk passage), and can take advantage of the existing metadata to query the vast HMI database. We use the hmi.sharp_cea_720s series, which comprise definitive data remapped to a cylindrical equal-area projection. Here, “definitive” means that the regions are identified and defined only after their full-disk passage. This series is chosen because it includes vector magnetic field data (and derived quantities) that may be included in our database in future—for example, to estimate the helicity of BMRs. However, for the initial database described in this paper, we use only the (remapped) line-of-sight magnetic field.

Our open-source Python code for BMR extraction is available online at https://github.com/antyeates1983/sharps-bmrs, and the resulting database generated for this paper is available on the Harvard Dataverse (Yeates, 2020). We stress that the philosophy in developing this new BMR database has been to minimise subjectivity through development of a fully-automated BMR fitting code.

2.1 Magnetogram Extraction

For each SHARP, we select a single representative magnetogram for our database. Each SHARP has multiple associated magnetograms, observed at high cadence over the region’s disk passage. Because we are using line-of-sight magnetic data, we select the observation with flux-weighted centroid closest to Central Meridian. This information is conveniently available in the SHARP metadata. One such magnetogram is illustrated in Figure 1(a).

Figure 1

An example SHARP region, number 7147 which passed Central Meridian on 30 September 2017. Panel (a) shows the HMI line-of-sight magnetogram with the provided SHARP mask indicated in black. Panel (b) shows the result of the Gaussian smoothing. Panel (c) shows the interpolated magnetogram on the computational grid, and panel (d) shows the approximating BMR. All panels are saturated at $\pm 100~\mathrm{G}$ (red positive, blue negative). This region is in the southern hemisphere, so is an “anti-Hale” region for Cycle 24. In (c) and (d), the overall centroid $(s_{0},\phi _{0})$ is shown by the black circle and the polarity centroids $(s_{+},\phi _{+})$, $(s_{-},\phi _{-})$ by the black crosses.

Having downloaded a single line-of-sight magnetogram for the given SHARP, we set pixels outside the masked region to zero. Then, we apply a Gaussian smoothing and interpolate the data to a lower resolution grid that is more appropriate for global simulations. The results in this paper, and the published database, use a uniform grid of $180\times 360$ cells in sine-latitude and longitude, although this resolution can be modified in the accompanying code, as can the parameters of the Gaussian filter. We took $\sigma =4$ pixels for the standard deviation of the Gaussian kernel, and used cubic interpolation. The results of these two stages are illustrated in Figure 1(b) and (c).

Once interpolated to the computational grid, we compute and record the flux imbalance, defined as the ratio of net flux to absolute flux over all pixels $i,j$,

$$ \Delta \Phi = \frac{\sum _{i,j}B^{i,j}}{\sum _{i,j}|B^{i,j}|}. $$

(1)

Note that all cells have equal area on this grid. The original SHARP regions are not flux balanced, so, once interpolated, regions will lie somewhere between $\Delta \Phi =0$ (perfectly flux balanced) and $\Delta \Phi =1$ (unipolar). Regions with large imbalance will be discarded from the resulting BMR database, as discussed in Section 2.3. For a region that is retained, the magnetogram is corrected for flux balance by applying different multiplicative scaling factors to the positive and negative pixels, in such a way that both the positive and the negative fluxes are scaled to the original mean of the two, and the overall unsigned flux is unchanged. This ensures that the polarity inversion lines do not change position.

2.2 Bipolar Approximation

To compute the approximating bipolar magnetic region (BMR) for a given SHARP, we first compute the centroids $(s_{+}, \phi _{+})$ and $(s_{-},\phi _{-})$ of positive and negative $B$ on the computational grid. Here $s$ denotes sine-latitude and $\phi$ denotes (Carrington) longitude. Based on these polarity centroids, we compute:

• i) the overall centroids

  $$ s_{0} = \frac{1}{2}(s_{+} + s_{-}),\qquad \phi _{0} = \frac{1}{2}(\phi _{+} + \phi _{-}), $$

  (2)

• ii) the polarity separation, which is the heliographic angle

  $$ \rho = \arccos \left [s_{+}s_{-} + \sqrt{1-s_{+}^{2}}\sqrt{1 - s_{-}^{2}} \cos (\phi _{+}-\phi _{-}) \right ], $$

  (3)

  and

• iii) the tilt angle with respect to the equator, given by

  $$ \gamma = \arctan \left [ \frac{\arcsin (s_{+}) - \arcsin (s_{-})}{\sqrt{1-s_{0}^{2}}(\phi _{-} - \phi _{+})} \right ]. $$

  (4)

Together with the unsigned flux, $|\Phi |$, these parameters define the BMR for our chosen functional form. For an untilted BMR centered at $s=\phi =0$, this functional form is defined as

$$ B(s,\phi ) = F(s,\phi ) = -B_{0}\frac{\phi }{\rho }\exp \left [- \frac{\phi ^{2} + 2\arcsin ^{2}(s)}{(a\rho )^{2}}\right ], $$

(5)

where the amplitude $B_{0}$ is scaled to match the corrected flux of the observed region on the computational grid. To account for the location $(s_{0},\phi _{0})$ and tilt $\gamma$ of a general region, we set $B(s,\phi ) = F(s',\phi ')$, where $(s',\phi ')$ are spherical coordinates in a frame where the region is centered at $s'=\phi '=0$ and untilted (explicit expressions are given in Appendix A). Figure 1(d) shows an example BMR.

The parameter $a$ in Equation 5 controls the size of the BMR relative to the separation, $\rho$, of the original polarity centroids. For given values of $\lambda _{0}$, $\gamma$, and $\rho$, and $B_{0}$ chosen to give the required magnetic flux, the parameter $a$ may be chosen to control the axial dipole moment of the BMR. In Section 2.5, we will see that a good match to the axial dipole moment of the original SHARP is obtained with $a=0.56$, and the same value works for every region.

2.3 Filtering

To ensure that only meaningful BMRs are included in the database, we filter out those SHARPs that do not meet the following criteria:

• i) The flux imbalance (before correction) should be less than a given threshold – we choose $\Delta \Phi \leq 0.5$. This effectively discounts unipolar SHARPs, which cannot be approximated by BMRs.

• ii) The polarity separation should be resolved on the computational grid, i.e., $\rho \ge \Delta \phi$ where $\Delta \phi$ is the longitudinal grid spacing.

The third and fourth columns of Table 1 show the number of SHARPs per year which were rejected at each of these two filtering steps. The figures in parentheses indicate total unsigned flux in Mx. Overall, about 13% of the SHARP flux is rejected as being unipolar (step i), and only a further 0.3% because of insufficient polarity separation (step ii). Following these two steps, we further identify and remove regions that correspond to repeat observations. These are shown in the fifth column of Table 1, and our procedure for identifying them is described next.

Table 1 Number of identified and filtered regions (parentheses show unsigned flux in Mx).

2.4 Removal of Repeat Observations

To avoid double-counting BMRs in our database, we have implemented a procedure to remove SHARPs that correspond to a repeat observation of a decaying region already in the database from a previous disk passage. This is done by comparing each SHARP with every existing SHARP that passed Central Meridian between 20 and 34 days earlier. To illustrate the procedure, Figure 2 shows a chain of three SHARPs where both successive pairs are candidate repeat regions.

Figure 2

Two successive pairs of candidate repeat regions: SHARPs 3563/3686 and 3686/3793. The left column shows the regions in context overlaid on HMI synoptic magnetograms (with the SHARP masks in black). The right column shows each SHARP on the computational grid. Dashed green contours in (b) and (d) show the derotated boundaries of the later SHARPs, 3686 and 3793, respectively. In all cases radial magnetic field $B$ is shown with red positive and blue negative, saturated at $\pm 200~\mathrm{G}$. The fluxes $\Phi _{1}$, $\Phi _{2}$, $\Phi _{3}$ give the total unsigned flux for each SHARP.

To determine automatically whether a SHARP $B^{i,j}_{n}$ is a repeat observation of an earlier region $B^{i,j}_{m}$, we first derotate $B^{i,j}_{n}$ to remove the effect of differential rotation, producing a SHARP $\overline{B}^{i,j}_{n}$ that may be directly compared with $B^{i,j}_{m}$. The derotation uses the (Carrington frame) angular velocity profile

$$ \Omega (s) = \Omega _{A} + \Omega _{B}s^{2} + \Omega _{C}s^{4}, $$

(6)

where $\Omega _{A}=0.18^{\circ }~\text{day}^{-1}$, $\Omega _{B}=-2.396^{\circ }~\text{day}^{-1}$, $\Omega _{C}=-1.787^{\circ }~\text{day}^{-1}$. We then classify $B^{i,j}_{n}$ as a repeat observation of $B^{i,j}_{m}$ if $B^{i,j}_{m}$ had higher unsigned flux than $\overline{B}^{i,j}_{n}$ in the latter’s envelope region $\overline{D}_{n}$. Specifically, if $R_{nm}>1$ where

$$ R_{nm} = \frac{\sum _{\overline{D}_{n}}|B^{i,j}_{m}|}{\sum _{\overline{D}_{n}} |\overline{B}^{i,j}_{n}|}. $$

(7)

For the two successive pairs in Figure 2, we find $R_{12}=0.94$ and $R_{23}=2.02$. Although most of the flux of SHARP 3563 lies within the derotated envelope of SHARP 3686 (dashed green line in Figure 2b), the ratio $R_{12}$ is less than unity because there is still rather strong flux present in SHARP 3686, and there is a sufficient mismatch in its (derotated) position compared to the earlier region. Thus SHARP 3686 is classified as a new region. On the other hand, $R_{23}$ is substantially greater than unity, so SHARP 3793 is classified as a repeat observation of SHARP 3686. Both decisions concur with a manual determination from the original HMI magnetograms in the left column.

The numbers of repeat regions identified in each year are shown in Table 1. There are 143 repeat regions identified in total, with about 12% of the original SHARP flux removed in this final stage of the filtering procedure. This leaves 1090 BMRs in the database, with a total unsigned flux of $1.0\times 10^{25}~\text{Mx}$.

2.5 Summary of Observed BMR Properties

Figure 3 summarizes the main properties of the 1090 BMRs in the database. In all three panels, the blue/red colour indicates the sign of the leading (westward) polarity, so that Hale’s law is evident in Figure 3. For comparison with previous studies, Figures 3(b) and (c) show flux versus polarity separation, and tilt versus latitude, respectively. Note that the tilt angle shown here has been adjusted to lie within the range $\pm 90^{\circ }$ by disregarding the magnetic polarity. Thus we plot the unsigned tilt angle

$$ \overline{\gamma } = \left \{ \textstyle\begin{array}{l@{\quad }l} \gamma + 180^{\circ }& \text{for}\ \gamma < -90^{\circ }, \\ \gamma & \text{for}\ -90^{\circ }\leq \gamma \leq 90^{\circ }, \\ \gamma - 180^{\circ }& \text{for}\ \gamma > 90^{\circ }. \end{array}\displaystyle \right . $$

(8)

Figure 3

Properties of the 1090 BMRs in the database, including (a) sine-latitude versus time; (b) unsigned flux versus polarity separation; and (c) tilt angle against latitude. In all three plots, points are coloured by unsigned flux, with the sign indicating the leading polarity. Panel (c) shows the unsigned tilt angle $\overline{\gamma }$. Black curves in (b) and (c) show least-squares fits as indicated (over all BMRs), while dashed curves show similar fits to the BMR databases of Wang and Sheeley (1989) for Cycle 21 and Yeates (2016) for Cycles 23 and 24. In (c) the fit coefficients are 0.52 for WS21 and 0.45 for Y23+24.

The least-squares fit in Figure 3(b) shows that the unsigned flux scales as $\rho ^{1.8}$. Here $|\Phi |$ is the total absolute flux over both polarities. This is consistent with flux being roughly proportional to the area of an active region. The dashed lines show fits for two other BMR catalogues: (i) for Cycle 21, derived from Kitt Peak Vacuum Telescope (KPVT) full-disk magnetograms (Wang and Sheeley, 1989; Sheeley and Wang, 2016), and (ii) for Cycles 23 and 24, derived from KPVT and later SOLIS/VSM synoptic magnetograms (Yeates, 2016). The steeper power law in our new database seems to arise from the scatter in flux at the small-$\rho$ end, which may be affected by our resolution and the thresholds applied to define SHARPs. As shown by Figure 3(c), our least-squares fit for tilt angle (Joy’s law) is close to that in the Yeates (2016) database. The slope found for the Wang and Sheeley (1989) BMRs is a little higher, but given the spread in tilt angles, it is difficult to be sure of the statistical significance of this difference. The fits are also broadly in line with Stenflo and Kosovichev (2012), who found $\overline{\gamma } = 32.1^{\circ }\sin \lambda _{0}$ for Cycle 23 using MDI data.

For our purposes in Section 3, an important BMR parameter is the axial dipole moment,

$$ b_{1,0} = \frac{3}{4\pi }\int _{0}^{2\pi }\int _{-1}^{1} sB(s,\phi )\, \mathrm{d}s\,\mathrm{d}\phi . $$

(9)

For each region, we record both $b_{1,0}^{\mathrm{Bi}}$—the axial dipole moment of the fitted BMR—and $b_{1,0}^{\mathrm{Si}}$—the axial dipole moment of the original SHARP region on the computational grid (as in Figure 1c). Figure 4(a) shows that there is a tight linear correlation between $b_{1,0}^{\mathrm{Bi}}$ and $b_{1,0}^{\mathrm{Si}}$, indicating that the fitted BMRs are well able to reproduce both the magnetic flux and the axial dipole moment of each individual SHARP. If the $a$ parameter in Equation 5 is reduced from its optimum value of 0.56, the linear relation remains but the slope reduces, as shown by the fainter points in Figure 4(a). Furthermore, as pointed out by Wang and Sheeley (1991), $b_{1,0}^{\mathrm{Bi}}$ depends linearly on the flux, cosine-latitude, and latitudinal spread of the BMR (cf. Jiang, Cameron, and Schüssler, 2014). Figure 4(b) shows that a similar relationship holds for our BMRs.

Figure 4

Scatter plots of BMR (initial) axial dipole moment against (a) the (initial) axial dipole moment of the original SHARP region, and (b) the parameter combination $|\Phi |\rho \sin \gamma \cos \lambda _{0}~R_{\odot }^{-2}$. Least-squares linear fits are indicated. In (a), the fainter points show the effect of reducing the BMR size in Equation 5 from $a=0.56$ to $a=0.4$.

Finally, summing $b_{1,0}^{\mathrm{Si}}$ over all regions gives a net axial dipole input of $3.32~\mathrm{G}$, with total positive/negative contributions of $2.15~\mathrm{G}$/$-0.64~\mathrm{G}$ in the northern hemisphere and $2.40~\mathrm{G}$/$-0.59~\mathrm{G}$ in the southern hemisphere. Thus we estimate a slightly higher total dipole input than Virtanen et al. (2019b), who found a net input of $2.91~\mathrm{G}$ for Cycle 24 when they extracted active regions from a combination of NSO/SOLIS and HMI synoptic maps. We conjecture that the slightly higher total may relate to more flux being included in our extraction technique using SHARP data. On the other hand, we reproduce the trends found by Virtanen et al. (2019b) whereby, for Cycle 24, the southern hemisphere has a stronger positive dipole input but a weaker negative dipole input than the northern hemisphere. In Section 3, we consider how these initial dipole moments would evolve over the solar cycle as the active regions spread out and decay.

3 Predicted Asymptotic Contributions

The contribution of an active region to the polar field at Solar Minimum is determined not by its axial dipole moment at emergence time, but by its eventual contribution to the axial dipole at Solar Minimum (e.g., Wang and Sheeley, 1991; Mackay, Priest, and Lockwood, 2002). In this section, we use surface flux transport (SFT) modelling to predict this net contribution for each region in our database. In particular, we consider how the bipolar approximation affects the accuracy of this prediction.

3.1 Surface Flux Transport Model

We evolve $B(s,\phi ,t)$ forward in time using the SFT model, which in our coordinates $s\equiv \sin \lambda$ and $\phi$ reads

$$ \frac{\partial B}{\partial t} = \frac{D}{R_{\odot }^{2}}\left [ \frac{\partial }{\partial s}\left ((1-s^{2}) \frac{\partial B}{\partial s}\right ) + \frac{1}{1-s^{2}} \frac{\partial ^{2}B}{\partial \phi ^{2}}\right ] - \frac{\partial }{\partial s}\left [\frac{v_{s}(s)}{R_{\odot }}\sqrt{1-s^{2}} B\right ]- \Omega (s)\frac{\partial B}{\partial \phi }. $$

(10)

Here $D$ is the supergranular diffusivity, $v_{s}(s)$ the meridional flow speed, and $\Omega (s)$ the angular velocity of differential rotation. In fact, it is not necessary to solve the two-dimensional equation, since $b_{1,0}$ depends only on the longitude-averaged field (DeVore, Boris, and Sheeley, 1984; Cameron and Schüssler, 2007; Iijima et al., 2017; Petrovay and Talafha, 2019). Specifically, writing $\overline{B}(s,t)=(2\pi )^{-1}\int _{0}^{2\pi }B(s,\phi ,t)\, \mathrm{d}\phi$, we have

$$ b_{1,0}(t) = \frac{3}{2}\int _{-1}^{1} s\overline{B}(s,t)\,\mathrm{d}s. $$

(11)

Taking the longitude-average of Equation 10 leaves the evolution equation

$$ \frac{\partial \overline{B}}{\partial t} = \frac{\partial }{\partial s}\left [\frac{D}{R_{\odot }^{2}}(1-s^{2}) \frac{\partial \overline{B}}{\partial s} - \frac{v_{s}(s)}{R_{\odot }} \sqrt{1-s^{2}}\overline{B}\right ], $$

(12)

which we solve numerically using a finite-volume method on our computational grid. We use the meridional flow profile

$$ v_{s}(s) = D_{u}s(1-s^{2})^{p/2}, $$

(13)

as used in the optimisation exercise of Whitbread, Yeates, and Muñoz-Jaramillo (2018). We adopt their meridional flow profile shape of $p=2.33$, but not their values of $D=466.8~\text{km}^{2}\,\text{s}^{-1}$ or $D_{u} = 0.025~\text{km}\,\text{s}^{-1}$. This is because these values were fitted to Cycles 21–23, and we find them to produce too strong a polar field with our Cycle 24 input data. We found a better match of the observed evolution—particularly of $b_{1,0}$—with $D$ reduced to $350~\text{km}^{2}\,\text{s}^{-1}$ and $D_{u}$ increased to $0.041~\text{km}\,\text{s}^{-1}$. The latter value corresponds to a peak poleward flow speed of $15~\text{m}\,\text{s}^{-1}$ (see Appendix B). We did not include an exponential decay time as this did not seem to be necessary to obtain a reasonable match to the observed evolution.

3.2 Complete Simulation

First, we illustrate the SFT model with all regions included. Figures 5(a) and (b) show the resulting time–latitude distributions of $\overline{B}$ in the SFT model where the sources are, respectively, the original SHARPs (mapped to the computational grid and with flux corrected, as in Figure 1c) or the fitted BMRs. For comparison, Figure 5(c) shows an observed supersynoptic map determined from HMI data. Specifically, we use the radial component, pole-filled synoptic maps in the hmi.synoptic_mr_polfil_720s series (Sun, 2018). A smoothing filter of the form $\mathrm{e}^{-b_{0}l(l+1)}$ is applied to the spherical harmonic coefficients (with $b_{0}=10^{-4}$) before mapping to the computational grid and correcting the flux balance multiplicatively. It is evident that both models do a reasonable job of matching the observed evolution of the large-scale field outside of active regions. Indeed, it is quite difficult to distinguish the two simulations visually (but compare, for example, the year 2014).

Figure 5

Full SFT simulations. Panels (a) and (b) show $\overline{B}$ as a function of time and latitude, capped at $\pm 10~{\text{G}}$, for the simulations with (a) SHARP regions, and (b) approximating BMRs. Panel (c) shows an observed supersynoptic map constructed from HMI data (see text), while panel (d) shows the axial dipole moments for the simulations and the HMI data (computed from panel c, with the early part being in good agreement with Sun et al., 2015). Dashed lines show simulations with increased $D_{u}$ (peak flow speed $20~\text{km}\,\text{s}^{-1}$).

If we plot the axial dipole moment, $b_{1,0}$, as in Figure 5(d), then the difference between the two simulations becomes clearer. The solid lines show simulations with our standard parameters, chosen so that the SHARPs simulation best reproduces the observed evolution of $b_{1,0}$. Although $b_{1,0}$ from the BMR simulation remains close until the time of dipole reversal in 2014, it later increases, and overestimates the axial dipole in 2020 substantially. Specifically, the net change in $b_{1,0}$ over the whole simulation is $+2.51~{\text{G}}$ for the SHARPs but $+3.11~{\text{G}}$ for the BMRs. This amounts to a $24\%$ overestimate by the BMR simulation. Conversely, if the meridional flow speed is chosen so that the BMR simulation matches the observed axial dipole – as shown by the dashed lines in Figure 5(d) – then the dipole reversal is delayed by several months and the evolution of $b_{1,0}$ is a poorer match for the observations between 2012 and 2017. To investigate why the BMR simulation overestimates the dipole moment, we next consider the evolution of the individual regions.

3.3 Individual Regions

Since the SFT equation 12 is linear, we can determine the dipole contribution of an individual active region by simulating that region alone (Yeates, Baker, and van Driel-Gesztelyi, 2015). When it is initialised with a single BMR and there is no further emergence, the $\overline{B}$ distribution will evolve toward an asymptotic steady state satisfying

$$ D(1-s^{2})\frac{\partial \overline{B}}{\partial s} - v_{s}(s)\sqrt{1-s^{2}} \overline{B} = \mathrm{constant}. $$

(14)

There is an exact solution of this equation giving the latitudinal profile (see Appendix B), but here we are interested in the magnitude of this asymptotic state, which we must compute numerically for each of the 1090 regions by evolving $\overline{B}$ forward in time for 10 years through Equation 12. This is ample time to reach the asymptotic profile, which typically takes only 2–3 years for these SFT parameters. From this asymptotic $\overline{B}$ profile, we compute the “final” axial dipole moment $b_{1,0}^{Bf}$ for each region numerically. As an example of this evolution, the rightmost column of Figure 6 shows the time evolution of $b_{1,0}$ for the five largest regions (by flux). For each region, we simulate both the evolution of the SHARP (first column) and of the fitted BMR (second column). For SHARPs 4698 and 1879, the dipole evolves in a similar way for both the SHARP and the approximating BMR. But for the other three regions in Figure 6, there are significant discrepancies.

Figure 6

The 5 regions with largest $|\Phi |$, showing their SHARPs (first column), approximating BMRs (second column), longitudinal averages (third column) and predicted evolution of their axial dipole moment with the SFT model (fourth column). In the third and fourth columns the red line shows the SHARP and the grey line the BMR.

Figure 7 shows a scatter plot of the predicted axial dipole moments using the BMRs versus the original SHARPs. Overall there is a strong correlation, but with some scatter. The slight bias for $b_{1,0}^{\mathrm{Bf}} > b_{1,0}^{\mathrm{Sf}}$ that leads to the overestimate in the complete simulation is present but difficult to discern by eye from this figure.

Figure 7

Scatterplot of predicted axial dipole moments after 10-year evolution of the BMRs versus the corresponding SHARPs.

In the SFT model, the dipole “amplification factor” for a BMR – in other words, the ratio $b_{1,0}^{\mathrm{Bf}}/b_{1,0}^{\mathrm{Bi}}$ – is known to be a Gaussian function of emergence latitude. This was first noted by Jiang, Cameron, and Schüssler (2014) and explained mathematically by Petrovay, Nagy, and Yeates (2020). Figure 8(a) shows that our BMR-driven simulations do indeed follow this pattern. On the other hand, when we plot the same ratio for the SHARP regions, as in Figure 8(b), there is still a Gaussian pattern but now considerable scatter. Similar scatter was found by Whitbread, Yeates, and Muñoz-Jaramillo (2018) who considered active regions from NSO data without making a BMR approximation. For clarity, outliers are not shown in Figure 8, and when these are included the average ratio $b_{1,0}^{\mathrm{Sf}}/b_{1,0}^{\mathrm{Si}}$ at each latitude falls below the Gaussian fit from the BMR simulation. Again, this reflects the tendency of the BMR-driven simulation to overestimate the asymptotic dipole moment.

Figure 8

Dipole amplification by the 10-year evolution for (a) the BMR approximation, and (b) the SHARPs, shown as a function of absolute latitude of the initial region. A least-squares fit to (a) is shown in both panels. In (b) there are a small number of extreme outliers not shown: the actual minimum value is −343.8 and the maximum is 75.3.

To understand why BMRs tend, on average, to overestimate the asymptotic dipole moment despite matching the initial dipole moment, it is instructive to look at the four regions with the greatest disparity, shown in Figure 9. From the first column, it is evident that the shapes of these regions diverge significantly from symmetric BMRs. In fact, it is only the longitude-average, $\overline{B}$, that matters, since this entirely determines the evolution of $b_{1,0}$. Accordingly, the third columns of Figures 6 and 9 show the initial $\overline{B}$ profiles for both the SHARP (in red) and the BMR (in grey). For example, SHARP 1879 in Figure 6 is well predicted by the BMR since the initial $\overline{B}$ is close to that of the SHARP, even though the original SHARP consists of (at least) two separate bipolar regions. On the other hand, all four regions in Figure 9 differ significantly from a symmetric BMR, even from the viewpoint of $\overline{B}$, leading to the disparity in predicted $b_{1,0}$ evolution.

Figure 9

The four regions with largest disparity $|b_{1,0}^{\mathrm{Bf}}-b_{1,0}^{\mathrm{Sf}}|$ in the predicted asymptotic dipole moment, showing their SHARPs (first column), approximating BMRs (second column), longitudinal averages (third column) and predicted evolution of their axial dipole moment (fourth column). In the third and fourth columns the red line shows the SHARP and the grey line the BMR.

As discussed by Iijima, Hotta, and Imada (2019), one way this can happen is if one polarity is more diffuse than the other; an example is SHARP 3376. Those authors argued that the effect will weaken the asymptotic dipole moment compared to the BMR, but here the effect is strong enough to reverse the sign. Overall, we find that the SHARP and BMR models predict final dipole moments of differing sign in 53 of our 1090 regions. More generally, the evolution can differ if $\overline{B}$ for the SHARP contains significant variation that cannot be captured by a BMR, true for all of the regions in Figure 9 (and several in Figure 6 too). Jiang et al. (2019) presented a case study of just such a $\delta$-spot region where the predicted dipole is weaker than expected, again changing sign during the evolution because of the more complex initial shape. In other cases, the SHARP predicts essentially zero dipole moment while the BMR predicts a non-zero value, or vice versa as in SHARP 1722.

Finally, Figure 10 represents an attempt to quantify the idea that regions with the greatest disparity between $b_{1,0}^{\mathrm{Bf}}$ and $b_{1,0}^{\mathrm{Sf}}$ tend to have either complex structure or asymmetry in the original SHARP. Specifically, Figure 10(a) shows—as a function of the disparity—the number of connected subregions with (absolute) field strength above $100~\mathrm{G}$ in the SHARP, which is a measure of complexity. Similarly, Figure 10(c) shows the difference between positive and negative peaks in the SHARP’s $\overline{B}$ profile, which is a measure of asymmetry. In both cases there are significant correlations, although stronger for the complexity measure than for the asymmetry. In addition, Figure 10(b) shows significant clusters of complex regions during the Solar Maximum 2014–2015, which is the period when the axial dipole moments from the SHARP and BMR simulations are seen to diverge in Figure 5(d).

Figure 10

Structure of the SHARPS, as measured by either the number of connected subregions with $|B^{i,j}|>100~\mathrm{G}$ within each SHARP (top), or the asymmetry $|\max |\overline{B}^{\,i}| - \min |\overline{B}^{\,i}||$ (bottom). The left column (a and c) shows each measure plotted against the disparity in predicted axial dipole, with $R$ showing the Spearman rank correlation coefficient (with negligible $p$ values). The right column (b and d) shows the SHARPs coloured by each measure as a function of latitude and time.

4 Conclusions

In summary, using a new BMR database derived from HMI/SHARP data, we find that the SFT simulation driven by BMRs overestimates the axial dipole at the end of the solar cycle by 24%. By contrast, using the SHARP maps themselves to drive the SFT model allows it to match the observed axial dipole evolution with the same SFT parameters.

Since the initial BMRs and SHARPs had exactly the same dipole moment, our results suggest that it is likely not the omission of active region inflows that is responsible for the overestimation of the axial dipole by BMR-driven models. In fact, one could argue that the effect of such inflows is already indirectly included in our SFT simulations, through our optimisation of the meridional flow speed. The observed meridional flow is found to be inversely proportional to solar activity (cf. Hathaway and Upton, 2014), consistent with the effect of active region inflows and likely explaining why we needed a faster flow speed for this Cycle 24 simulation compared to the Whitbread, Yeates, and Muñoz-Jaramillo (2018) model.

We saw in Figure 5(d) that changing the meridional flow speed can give the correct axial dipole at the end of the BMR-driven simulation, but at the expense of unrealistic intermediate evolution (the dashed grey curve in Figure 5d). In particular, this delays the axial dipole reversal by a few months. In their recent parameter study for the SFT model, where the combined source term is represented by a pair of flux rings, Petrovay and Talafha (2019) found that an exponential decay term was required in order to avoid such a delay in the dipole reversal. Here we find that such a term is not required in order to reverse the dipole at the correct time, provided the original SHARPs are used rather than BMRs.

Could the accuracy of the SFT model be improved within the context of BMR sources? To mimic past simulations (e.g., Yeates, 2014), we also tried a run with all of the BMR tilts reduced from $\tilde{\gamma }$ to $g\tilde{\gamma }$. With $g=0.82$, the resulting simulation reproduced the correct final dipole moment, but again at the expense of the intermediate evolution – in fact, the resulting curve for $b_{1,0}$ was close to the dashed grey curve in Figure 5(d). It is possible that improved results could be obtained by selectively reducing the tilt angles according to the complexity of individual regions, but this remains for future investigation. In fact, it may be more effective to look at ways of fitting highly asymmetric or multipolar SHARPs (such as those shown in Figure 9) with asymmetric BMRs, or indeed with multiple BMRs. The manual databases of Wang and Sheeley (1989) and Yeates (2016) were able to fit multiple BMRs to the same active region where appropriate, but an objective, automated procedure for this remains to be developed.

Appendices

Appendix A: Coordinate Rotation

Here we derive the coordinate transformation from the frame $(s,\phi )$ where the BMR is centred at $s=\phi =0$ to the frame $(s',\phi ')$ where it is centred at $(s_{0},\phi _{0})$ with tilt $\gamma$. This amounts to a rotation, which is easiest to express in Cartesian coordinates

$$ x = \cos \phi \sqrt{1-s^{2}}, \quad y=\sin \phi \sqrt{1-s^{2}},\quad z=s. $$

(15)

Multiplying by the rotation matrices for the sequence of rotations indicated in Figure 11 shows that Cartesian coordinates in the rotated frame are

$$\begin{aligned} \left[ \textstyle\begin{array}{c} x' \\ y' \\ z' \end{array}\displaystyle \right ] =& \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} 1\; & 0\; & 0 \\ 0\; & \cos \gamma \; & -\sin \gamma \\ 0\; & \sin \gamma \; & \cos \gamma \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \cos \lambda _{0}\; & 0\; & \sin \lambda _{0} \\ 0\; & 1\; & 0 \\ -\sin \lambda _{0}\; & 0\; & \cos \lambda _{0} \end{array}\displaystyle \right ] \cdot \end{aligned}$$

(16)

$$\begin{aligned} &\cdot \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \cos \phi _{0}\; & \sin \phi _{0}\; & 0 \\ -\sin \phi _{0}\; & \cos \phi _{0}\; & 0 \\ 0 & 0 & 1 \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} x \\ y \\ z \end{array}\displaystyle \right ], \end{aligned}$$

(17)

where $s_{0}=\sin \lambda _{0}$. From these we determine $\phi '=\arctan (y'/x')$ and $s' = z'$.

Figure 11

The sequence of rotations to transform to coordinates where the BMR is at $s'=\phi '=0$ with $\gamma '=0$: (i) rotate angle $\phi _{0}$ around the $z$-axis; (ii) rotate angle $\lambda _{0}=\arcsin (s_{0})$ around the (new) $y$-axis; and (iii) rotate angle $\gamma$ around the (new) $x$-axis.

Appendix B: Analytical Expressions

Here we note some analytical expressions for a meridional flow of the form in Equation 13, which the reader may find useful (see also DeVore, Boris, and Sheeley, 1984). Firstly, one can show that the flow peaks at latitudes $s=\pm (1+p)^{-1/2}$ with speed

$$ v_{\mathrm{max}} = \pm D_{u} p^{p/2}(1+p)^{-(1+p)/2}. $$

(18)

Secondly, one can solve Equation 14 exactly to find the shape of the asymptotic profile, which is

$$ \overline{B}(s) = \overline{B}(1) \frac{\exp [-C(1-s^{2})^{(1+p)/2}] - \mathrm{e}^{-C}}{1 - \mathrm{e}^{-C}}, \qquad C=\frac{D_{u}R_{\odot }}{(p+1)D}. $$

(19)

Notice that the diffusivity $D$ enters only through the dimensionless parameter $C$. The actual amplitude of the profile, written here in terms of the value $\overline{B}(1)$ at the north pole, must be determined by following the evolution (but for an approximate estimate of the asymptotic flux of a single active region, see Petrovay, Nagy, and Yeates, 2020). Once the asymptotic value of $\overline{B}(1)$ is known, the axial dipole moment is found from Equation 19 to be

$$ b_{1,0} = \frac{3\overline{B}_{1}}{2(1-\mathrm{e}^{-C})}\left [q \frac{ \Gamma (q) - \Gamma (q,C)}{C^{q}} - \mathrm{e}^{-C}\right ], $$

(20)

where $q=2(1+p)^{-1}$ and $\Gamma (q,C)$ is the incomplete Gamma function. In this paper, $b_{1,0}$ was simply computed numerically from the asymptotic solution $\overline{B}(s)$.