How good is the bipolar approximation of active regions for surface flux transport?

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.


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;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 . 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 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 BMRs 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, 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 in Section 4.

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.

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).
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 × 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 σ = 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, 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 ∆Φ = 0 (perfectly flux balanced) and ∆Φ = 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 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.

Bipolar approximation
To compute the approximating bipolar magnetic region (BMR) for a given SHARP, we first compute the centroids (s + , φ + ) and (s − , φ − ) of positive and negative B on the computational grid. Here s denotes sine-latitude and φ denotes (Carrington) longitude. Based on these polarity centroids, we compute (i) the overall centroids (ii) the polarity separation, which is the heliographic angle and (iii) the tilt angle with respect to the equator, given by .  (d) shows the approximating BMR. All panels are saturated at ±100 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 , φ 0 ) is shown by the black circle and the polarity centroids (s + , φ + ), (s − , φ − ) by the black crosses.
Together with the unsigned flux, |Φ|, these parameters define the BMR for our chosen functional form. For an untilted BMR centered at s = φ = 0, this functional form is defined as 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 , φ 0 ) and tilt γ of a general region, we set coordinates in a frame where the region is centered at s ′ = φ ′ = 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, ρ, of the original polarity centroids. For given values of λ 0 , γ, and ρ, 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.

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 ∆Φ ≤ 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., ρ ≥ ∆φ where ∆φ 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.

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  (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 ±200 G. The fluxes Φ 1 , Φ 2 , Φ 3 give the total unsigned flux for each SHARP.
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.
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 B i,j n that may be directly compared with B i,j m . The derotation uses the (Carrington frame) angular velocity profile where 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 × 10 25 Mx. 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 ±90 • by disregarding the magnetic polarity. Thus we plot the unsigned tilt angle

Summary of observed BMR properties
The least-squares fit in Figure 3(b) shows that the unsigned flux scales as ρ 1.8 . Here |Φ| 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-ρ 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 γ = 32.1 • sin λ 0 for Cycle 23 using MDI data.
For our purposes in Section 3, an important BMR parameter is the axial dipole moment, For each region, we record both b Bi 1,0 -the axial dipole moment of the fitted BMR -and b Si 1,0 -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 Bi 1,0 and b Si 1,0 , 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 Bi 1,0 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.
Finally, summing b Si 1,0 over all regions gives a net axial dipole input of 3.32 G, with total positive/negative contributions of 2.15 G/−0.64 G in the Northern hemisphere and 2.40 G/−0.59 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 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.

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.

Surface flux transport model
We evolve B(s, φ, t) forward in time using the SFT model, which in our coordinates s ≡ sin λ and φ reads (10) Here D is the supergranular diffusivity, v s (s) the meridional flow speed, and Ω(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 longitudeaveraged field (DeVore, Boris, and Sheeley, 1984;Cameron and Schüssler, 2007;Iijima et al., 2017;Petrovay and Talafha, 2019). Specifically, writing B(s, t) = Taking the longitude-average of Equation (10) leaves the evolution equation which we solve numerically using a finite-volume method on our computational grid. We use the meridional flow profile 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 km 2 s −1 or D u = 0.025 km 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 km 2 s −1 and D u increased to 0.041 km s −1 . The latter value corresponds to a peak poleward flow speed of 15 m 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.

Complete simulation
First, we illustrate the SFT model with all regions included. Figures 5(a) and (b) show the resulting time-latitude distributions of 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 e −b0l(l+1) is applied to the spherical harmonic coefficients (with b 0 = 10 −4 ) before mapping to the our 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). 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 G for the SHARPs but +3.11 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.

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 B distribution will evolve toward an asymptotic steady state satisfying 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 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 B profile, we compute the "final" axial dipole moment b Bf 1,0 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 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 Bf 1,0 > b Sf 1,0 that leads to the overestimate in the complete simulation is present but difficult to discern by eye from this figure.
In the SFT model, the dipole "amplification factor" for a BMR -in other words, the ratio b Bf 1,0 /b Bi 1,0 -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 Sf 1,0 /b Si 1,0 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.
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. first column, it is evident that the shapes of these regions diverge significantly from symmetric BMRs. In fact, it is only the longitude-average, 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 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 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 B, leading to the disparity in predicted b 1,0 evolution. 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 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 δ-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 Bf 1,0 and b Sf 1,0 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 G in the SHARP, which is a measure of complexity. Similarly, Figure 10 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).

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 allow 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 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γ to gγ. 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 gray curve in Figure 5(d). It is possible that improved results could be obtained by selectively reducing the tilt angles according to the complexity Figure 11. The sequence of rotations to transform to coordinates where the BMR is at s ′ = φ ′ = 0 with γ ′ = 0: (i) rotate angle φ 0 around the z-axis; (ii) rotate angle λ 0 = arcsin(s 0 ) around the (new) y-axis; and (iii) rotate angle γ around the (new) x-axis.
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. one can show that the flow peaks at latitudes s = ±(1 + p) −1/2 with speed v max = ±D u p p/2 (1 + p) −(1+p)/2 .
Secondly, one can solve the ODE (14) exactly to find the shape of the asymptotic profile, which is 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 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 B(1) is known, the axial dipole moment is found from (19) to be where q = 2(1 + p) −1 and Γ(q, C) is the incomplete Gamma function. In this paper, b 1,0 was simply computed numerically from the asymptotic solution B(s).