Evolution of the Sun’s nonaxisymmetric toroidal field
^{1} Institut für Astrophysik, GeorgAugustUniversität Göttingen, 37077 Göttingen, Germany
email: belda@mps.mpg.de
^{2} MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
Received: 27 January 2017
Accepted: 28 March 2017
Aims. We aim to infer the subsurface distribution of the Sun’s nonaxisymmetric azimuthal magnetic flux from observable quantities, such as the surface magnetic field and the large scale plasma flows.
Methods. We have built a kinematic flux transport model of the solar dynamo based on the BabcockLeighton framework. We constructed the source term for the poloidal field using SOLIS magnetograms spanning three solar cycles. Based on this source we calculated the azimuthal flux below the surface. The flux transport model has two free parameters which we constrained using sunspot observations from cycle 22. We compared the model results with observations from cycle 23.
Results. The structure of the azimuthal field is mainly axisymmetric. The departures from axisymmetry represent, on average, ~3% of the total azimuthal flux. Owing to its relative weakness, the nonaxisymmetric structure of the azimuthal field does not have a significant impact on the location in which the emergences appear or on the amount of flux contained in them. We find that the probability of emergence is a function of the ratio between the flux content of an active region and the underlying azimuthal flux.
Key words: Sun: magnetic fields / Sun: activity
© ESO, 2017
1. Introduction
The magnetic activity of the Sun and other stars is a manifestation of their internal magnetic field, which is thought to be sustained by a hydromagnetic dynamo. In the case of the Sun, it is generally thought that the differential rotation in the convection zone generates the toroidal magnetic field out of the poloidal field, but where exactly this field is amplified and stored is still an open question (see, e.g. Charbonneau 2010). The mechanism for the regeneration of the poloidal field from the azimuthal component is less agreed upon, with modelling approaches falling mainly into two categories: the turbulent dynamo models and the BabcockLeighton models.
In BabcockLeighton models (Babcock 1961; Leighton 1969), the poloidal field is regenerated by the surface transport of the magnetic flux of decaying active regions. Newly emerged bipolar magnetic regions (BMRs) show a systematic tilt with respect to the E−W direction, with the preceding polarity (in the Sun’s sense of rotation) appearing closer to the equator than the trailing polarity (Joy’s law). In addition, the preceding polarity of a BMR emerging in a given hemisphere tends to be of the same sign as the polar field in that hemisphere at the beginning of the ongoing activity cycle (Hale’s law). This facilitates the crossequatorial transport of preceding polarity flux, and leads to the cancellation of the polar fields and the buildup of a new, reversed axial dipole, which is the source of azimuthal field for the new activity cycle.
BabcockLeighton models have gained substantial support in recent years. DasiEspuig et al. (2010) found a strong correlation between the strength of a BabcockLeighton type source term in a given cycle, calculated from the observed tilt angle of active regions, and the strength of the next cycle. Kitchatinov & Olemskoy (2011) found that the aggregate contribution of active regions to the poloidal field during one cycle and the strength of the global dipole at the end of the same cycle (as inferred from the AA index) also correlate closely. Wang et al. (2009) showed that the buildup of the polar fields during cycles 20 to 23 is consistent with the passive transport of magnetic flux by the observed surface flows. On the theoretical side, Cameron & Schüssler (2015) showed that the main source of net azimuthal flux in each hemisphere is the winding up of poloidal flux that is connected to the polar fields at the surface.
One of the key features of dynamo models is the formulation of the poloidal source term as a function of the azimuthal field. One often considered possibility in BabcockLeighton flux transport models is that magnetic flux tubes are stored in an overshoot region at the base of the convection zone. These develop a magnetic buoyancy instability and rise through the convection zone to emerge at the surface in the form of BMRs. For a review on this topic, see, for example, Fan (2009).
Alternatively, 3D numerical calculations indicate that persistent, coherent azimuthal magnetic structures can arise in a turbulent convection zone, owing to turbulent intermittency (Brown et al. 2010; Nelson et al. 2013, 2014). Moreover, the transport of magnetic flux to the surface can be achieved by means of convective upflows, which might be at least as relevant as magnetic buoyancy.
In this work we consider the evolution of the radially integrated azimuthal flux density as a function of longitude and latitude. Our aim is to infer the distribution and evolution of the subsurface azimuthal flux from observable quantities in order to gain insight on its relation with the observed properties of active regions. To do so, we have constructed a model of the magnetic flux transport in the Sun, based on the BabcockLeighton framework. The poloidal field source term is determined by observational data (synoptic magnetograms). The paper is structured as follows: in Sect. 2 we introduce and calibrate our model; in Sect. 3 we present and discuss our results; and in Sect. 4 we briefly summarize our conclusions. In the appendices we derive the equation for the evolution of the azimuthal flux density and other supplementary quantities in our model.
2. Methods
2.1. Model
We have considered a mean field approach based on horizontal averages, (1)where j = { r,θ,φ }, (2)and δ is the scale over which the average is performed. Using δ ~ 2–3° is enough to ensure some scale separation with respect to the larger turbulence correlation lengths at the surface (those of supergranulation). In the remainder of the paper we drop the angle brackets for clarity, and refer to the jth component of the averaged magnetic field by B_{j}.
Our model consists of two twodimensional domains in the (φ,θ) plane, representing the surface of the Sun and the convection zone, respectively. The evolution of the surface magnetic field, assumed to be radial, is governed by the surface flux transport equation (DeVore et al. 1984): (3)where Ω_{R⊙} is the differential rotation, u_{M} is the meridional flow and η_{H} is the surface diffusivity associated to the convective flows. The emergence of new flux is described by the source term S(θ,φ,t).
The first term in Eq. (3)describes the transport of the surface field by the solar differential rotation. We use the differential rotation profile inferred by correlation tracking of magnetic features by Hathaway & Rightmire (2011): (4)The second term in Eq. (3)corresponds to the surface meridional flow. Following van Ballegooijen et al. (1998), we model the meridional flow as: (5)This expression captures the main characteristics of the observed meridional flow (see Hathaway & Rightmire 2011).
The third term in Eq. (3)describes the dispersal of magnetic flux on the surface by means of random convective flows modelled as a diffusion process (see Leighton 1964; MartinBelda & Cameron 2016). We use a surface diffusivity of η_{H} = 250 km^{2} s^{1}, as indicated by observations (Schrijver & Martin 1990; Jafarzadeh et al. 2014).
The source term S(θ,φ,t) represents new emergences, and is built from synoptic magnetograms (see Sect. 2.2).
The azimuthal field in the convection zone is represented in our model by the azimuthal flux per unit colatitude, that is: (6)In the above expression, R_{b} refers to the bottom of the convection zone and R_{⊙} is the solar radius.
Following Cameron & Schüssler (2017), we made the following assumptions regarding the plasma flows and the structure of the internal magnetic field of the Sun:

1.
The magnetic field is purely radial in the nearsurface shearlayer (NSSL), owing to strong downwards turbulent pumping.

2.
The magnetic field does not penetrate the radiative interior.

3.
The poloidal field does not penetrate the tachocline. This assumption of the model is partly justified by Spruit (2011), who noted that the tachocline cannot support large shear stresses, which would be present if the poloidal field did penetrate the tachocline.

4.
The radial shear is negligible in the region between the tachocline and the NSSL. This is based on helioseismic inference of the rotation rate in the deep interior and the NSSL (ChristensenDalsgaard & Schou 1988).
To derive an evolution equation for b_{φ}, we integrated the azimuthal component of the induction equation in radius. The resulting equation reads (see Appendix A): (7)where Ω_{RNSSL} is the differential rotation at the bottom of the NSSL, is an effective return meridional flow, η_{0} is the effective diffusivity of the azimuthal field, S_{φ}(θ,φ,t) is a source term associated to flux emergence, and (8)In the above expression, R_{T} refers to the top of the tachocline and R_{NSSL} refers to the bottom of the NSSL. The quantity b_{θ} is calculated in terms of b_{φ} and B_{r} in Appendix B.
The first three terms in Eq. (7)describe the generation of azimuthal flux by differential rotation and the azimuthal flux transport. Here, B_{r} is the radial field at the surface, which is constrained by observation and evolves according to Eq. (3). The radially integrated θcomponent of the magnetic field, b_{θ}, can be obtained in terms of b_{φ} and B_{r} from the solenoidality condition ∇·B = 0 (see Appendix B). The differential rotation profile is evaluated at the bottom of the NSSL. The analysis of helioseismic data by Barekat et al. (2014) suggests that the radial shear in this layer is independent of latitude. Following these authors, we adopted (9)The fourth term of Eq. (7)describes the turbulent diffusion of azimuthal flux. Following Cameron & Schüssler (2016), we assumed the following form for the diffusivity in the derivation of Eq. (7): (10)where η_{0} is a free parameter of our model. Cameron & Schüssler (2016) used the properties of the decay phase of the sunspot cycles to estimate η_{0} ~ 150−450 km^{2} s^{1}.
The fifth term corresponds to the advection of the azimuthal flux by an effective equatorward flow, which we modeled as (11)where u_{0} is a free parameter of the model. This flow can correspond to a return meridional flow (Wang et al. 1991; Durney 1995) or equatorward pumping (Guerrero & de Gouveia Dal Pino 2008).
The solenoidality condition requires that the surface magnetic field connects to the subsurface field. Hence, a source term in Eq. (7)is needed to ensure the connectivity of the surface field sources to the field in the convection zone. In order to calculate S_{φ}(θ,φ,t), we extrapolated the surface sources downward via a potential field solution (see Sect. 2.2).
Our model, therefore, consists of: (a) a twodimensional domain representing the surface of the Sun, in which the (radial) surface field evolves according to Eq. (3); (b) a twodimensional domain representing the convection zone, in which the radial integral of the azimuthal magnetic field evolves according to Eq. (7); and (c) the coupling of both domains through the emergences, represented by the source terms, and the solenoidality condition.
2.2. Treatment of the source terms
The source term in Eq. (3), which represents the emergence of flux on the solar surface, was calculated using SOLIS synoptic magnetograms. Let be the radial magnetic field corresponding to Carrington rotation (CR) n, as given by the corresponding synoptic magnetogram. The field associated with emergences during CR n was computed as (12)where is the magnetic field from the previous magnetogram evolved for one rotation using Eq. (3). This expression is related to the surface source term S(θ,φ,t) through (13)In practice, we add ΔB_{r} to the simulated surface field every carrington rotation. The synoptic magnetograms are corrected by multiplying the positive values by a factor such that the resulting net magnetic flux on the surface is zero. Apart from describing flux emergences, the surface source term also corrects for the errors of the SFT model and errors in individual synoptic magnetograms.
The surface field B_{r} must connect with the field in the convection zone so that ∇·B = 0 is maintained. The connectivity of the surface sources with the subsurface field is achieved through the source term in Eq. (7), which we calculated by performing a downwards potential field extrapolation of ΔB_{r}(θ,φ,t^{n}). Let ΔB_{φ}(θ,φ,t^{n}) be the azimuthal component of the extrapolated field, and . The source term in Eq. (7), S_{φ}(θ,φ,t^{n}), is related to Δb_{φ}(θ,φ,t^{n}) through (14)The extrapolated Δb_{φ} is added to the simulated b_{φ} every rotation, at the same time ΔB_{r} is added to the surface field.
The potential field extrapolation of ΔB_{r}(θ,φ,t^{n}) is more easily done using spherical harmonics, for which it is convenient to remap the synoptic magnetograms onto a grid equally spaced in θ, rather than in cosθ. The remapping prevents large errors near the poles, which arise from the poorer spatial resolution of the synoptic magnetograms at high latitudes. A discussion of this problem can be found in Tóth et al. (2011).
The value of the magnetic field at each point of the new equiangular grid was interpolated linearly from the magnetogram (old grid). The magnetic field at the poles is not known, which makes it appropriate to perform the interpolation in Fourier space by expanding (15)At the poles, regularity of B_{r} translates into the following boundary conditions for a_{m}: The value of the Fourier coefficient at point of colatitude θ^{∗} of the new grid is given by: (18)where θ^{−} and θ^{+} are the points of the old grid adjacent to θ^{∗} in the θ direction. For the points of the new grid that are located between the north pole and the northernmost grid point of the old grid, , Eqs. (16)and (17)lead to: For grid points between the southernmost point of the old grid, , and the south pole, we have: Since the maximum angular order of the spherical harmonic analysis is limited by an antialiasing condition, the decomposition of B_{r} in spherical harmonics has the effect of slightly smoothing the magnetograms (see the top row of Fig. 1).
Fig. 1 Initial condition of the simulations. Top left: synoptic magnetogram corresponding to CR 1625. Top right: synoptic magnetogram corresponding to CR 1625, remapped to an equiangular grid and resampled to the highest angular degree order, l, used to compute the potential field extrapolation. Bottom left: radial integral of the azimuthal field across the convection zone, b_{φ}, extrapolated from the surface field. Bottom right: radial integral of the θ component of the magnetic field across the convection zone, b_{θ}, extrapolated from the surface field. Red and blue indicate opposite polarities in all the maps. 

Open with DEXTER 
2.3. Setup and calibration
The initial condition of the simulations was computed as a potential field extrapolation of the first magnetogram of the series. This was taken at CR 1625, which corresponds to the end of cycle 20 in March, 1975. Figure 1 shows the raw magnetogram, the remapped version used to build the surface field source, and the extrapolated b_{φ} and b_{θ}. The sign of b_{φ} indicates the direction of the integrated azimuthal field. In our chosen coordinate system, a positive value of b_{φ} corresponds to azimuthal field pointing in the sense of rotation of the Sun. The sign of b_{θ} is predominantly negative, which reflects the sign of the axial dipole at this cycle minimum. The irregularities of b_{θ} near the north pole are probably related to the noise of the magnetogram, and diffuse very quickly once the simulation starts.
To calibrate the model, we required that the simulated azimuthal flux lie radially below the active regions observed on the Sun during cycle 22. We ran simulations with different values of the free parameters for cycles 21, 22 and 23. The strength of the return flow, u, was varied between 1 and 6 m s^{1}, and the diffusivity η_{0} was varied between 25 and 600 km^{2} s^{1}. The simulations were let run for cycle 21 for initialization. The active region data (latitude, longitude and area) was extracted from the USAF/NOAA sunspot group database. We evaluate the azimuthal flux density underlying every emergence location, b_{φ,i} = b_{φ}(λ_{i},φ_{i}), as close in time as possible, but always prior, to the time of observed maximum area of the active region (here the index i runs over the emergences). In the case of backside emergences, there can be a significant delay between these two times (as large as half a rotation). We note, however, that the change in azimuthal flux density on timescales shorter than one rotation is, in most cases, small as the nonaxisymmetric component represents less than 3% of the total azimuthal flux density on average.
Figure 2 shows a few examples of the distribution of emergences according to their underlying azimuthal flux density for various combinations of parameters. Different values of η_{0} give rise to different global magnitudes of the azimuthal flux density so, for easier comparison, we normalized b_{φ,i} to its maximum value in each simulation. The emergences were grouped in bins of width 4 × 10^{2} to reduce noise. The two peaks in the distribution of the active regions reflect the equatorial antisymmetry of the azimuthal flux. In cycle 22, the underlying azimuthal flux is mainly positive in the northern hemisphere and negative in the southern hemisphere. The case with u_{0} = 3 m s^{1} and η_{0} = 100 km s^{2} (black line) corresponds best to the requirement that the simulated azimuthal flux be preferentially located underneath the emergences. In the other three cases, the twopeak structure is not so conspicuous, and there are more emergences where there is little or no simulated azimuthal flux.
To find the parameter combinations that yield two wellseparated peaks, we considered the quantity (23)where c runs through the b_{φ,i} bins and b_{φ,c} is the midpoint of each bin. The value of ξ will be bigger for the simulations where the emergences occur farther away from the places where b_{φ} ~ 0. Figure 3 shows the value of ξ for all test runs carried out. We find that the combination u = 3 m s^{1} and η_{0} = 100 km^{2} s^{1} maximizes ξ. These parameters are close to the range found by Cameron & Schüssler (2017) for the operation of the solar dynamo in an updated version of the BabcockLeighton model. We proceed to the analysis of the data from cycle 23 using the azimuthal flux density maps generated in the simulation using the above parameter values.
Fig. 2 Number of emergences (N) as a function of the normalized underlying azimuthal flux density (). The different curves correspond to different choices of the free parameters. Continuous black line: u_{0} = 3 m s^{1} and η_{0} = 100 km^{2} s^{1}. Dashed blue line: u_{0} = 5 m s^{1} and η_{0} = 400 km^{2} s^{1}. Dashed red line: u_{0} = 1 m s^{1} and η_{0} = 25 km^{2} s^{1}. Continuous green line: u_{0} = 6 m s^{1} and η_{0} = 600 km^{2} s^{1}. 

Open with DEXTER 
3. Results for cycle 23
3.1. Angular distribution and evolution of azimuthal flux
Our analysis, which integrates the azimuthal field in the radial direction, allows us to infer the latitudinal and longitudinal structure of the subsurface field from the observed surface field and large scale flows.
The top row of Fig. 4 shows the observed surface field near the activity maximum of Cycle 23 (CR 1987) and towards the end of that cycle (CR 2060). Squares indicate the emergence sites from the USAF/NOAA sunspot record. Some of the emergences do not seem to correspond to strong concentrations of magnetic field in the magnetograms, and some features in the magnetograms do not have a counterpart in the active region record. A possible cause for the mismatch in the first case could be the loss of information in the lowresolution magnetograms. In the second case, one possibility is that small sunspot groups that emerged on the far side of the Sun lacked spots when the region rotated onto the visible side. In this case, the flux content of the active regions is still present in the synoptic magnetogram, and therefore included in the source term.
The middle row of Fig. 4 shows the inferred maps of azimuthal flux density. The magnetic activity sits mainly on top of the azimuthal flux system. The azimuthal flux corresponding to CR 1987 presents a structure that is strongly axisymmetric and antisymmetric about the equator. The strongest concentration of azimuthal flux occurs at ~ 15° of latitude in both hemispheres. At CR 2060, most of the azimuthal flux has diffused and cancelled across the equator, and a new azimuthal flux system of opposite polarity, corresponding to the new cycle, has begun to develop at higher latitudes from the windingup of the reversed poloidal field.
The bottom row of Fig. 4 shows the nonaxisymmetric part of the integrated azimuthal field, calculated as , where ⟨ b_{φ} ⟩ is the azimuthal mean of b_{φ}. The magnitude of represents, on average, ~ 3% of the total azimuthal flux density. This nonaxisymmetric structure arises from the emergence process (which is the only nonaxisymmetric ingredient of our model), and tends to diffuse away towards the end of the cycle, when the number of emergences is smaller.
3.2. Impact of the nonaxisymmetric structure on the emergence process
In order to investigate whether the nonaxisymmetric structure of the azimuthal flux influences the emergence process, we consider the deviation of the azimuthal flux density underlying each active region from the azimuthally averaged azimuthal flux density at the latitude of emergence, ⟨ b_{φ,i} ⟩. The result is shown in Fig. 5. As seen in Sect. 2.3, the bipolar distribution of events reflects the strong antisymmetry of the azimuthal field about the equator. The active regions for which ⟨ b_{φ,i} ⟩ < 0 are mainly located in the north hemisphere, while those with ⟨ b_{φ,i} ⟩ > 0 correspond to the south hemisphere.
Fig. 3 Distribution of the quantity as a function of the parameters u_{0} and η_{0}. Each circle represents a run. The size of each circle is proportional to the quantity ξ, which measures the adjustment of the run to our requirement that the simulated azimuthal flux lies underneath the observed active regions. For better visualization, we also encode the value of ξ in the colour of the circles. 

Open with DEXTER 
Fig. 4 Surface field (B_{r}), azimuthal flux density (b_{φ}) and its nonaxisymmetric component (), for CR 1987 (close to the middle of cycle 23) and CR 2060 (near the end of that cycle). Red and blue represent opposite polarities. The squares on the top and middle rows represent observed emergence sites, extracted from the USAF/NOAA sunspot group database. 

Open with DEXTER 
An influence of the nonaxisymmetric structure of the azimuthal flux on the triggering of the emergence process would lead to a nonzero average value of in each hemisphere. For example, if emergences at a given latitude tended to occur at longitudes where the azimuthal flux density is above the azimuthal mean, averaging over all the emergences for which ⟨ b_{φ,i} ⟩ > 0 would yield a positive value. In the other hemisphere, the average would be smaller than zero. Computing these averages yields The average value of in each hemisphere results very close to zero, in relative terms. Therefore, we do not find a significant correlation between the location of the emergence events and the departures from axisymmetry of the subsurface azimuthal flux.
Next, we consider the possible influence of the nonaxisymmetric structure on the active region areas. Figure 6 shows the excess of azimuthal flux density above the azimuthal average beneath each active region, , versus the active region area, A_{i}. Again, there is no significant deviation from zero, which suggests that the inferred nonaxisymmetric structure of the subsurface azimuthal flux is unrelated to the area of the emerged active region.
3.3. Relationship between azimuthal and emerged flux
Here we study the relationship between the flux contained in an active region and the azimuthal flux density underneath the emergence site. Figure 7 shows the flux of each active region in the sunspot group record at the time of maximum development, Φ_{i}, versus the unsigned underlying azimuthal flux density,  b_{φ,i} . Since we used SOLIS synoptic magnetograms to feed our simulations, we want to compare with fluxes comparable to those from SOLIS. The flux contained in each active region was calculated from its sunspot group area (obtained from the USAF/NOAA sunspot database) by using the crosscalibration factors in Table 2 of MuñozJaramillo et al. (2015). The resulting relationship is Φ_{i} [Mx] = 1.44 × 10^{19}A_{i} [μHem]. A factor 1/2 is introduced to account for the fact that the two polarities of the active region are part of a single Ωshaped magnetic structure that crosses the solar surface twice.
Fig. 5 Nonaxisymmetric component of the azimuthal flux underlying each emergence in cycle 23 () versus its axisymmetric component (⟨ b_{φ,i} ⟩). The emergences are represented by green points. The black line represents the average of over all emergences inside bins of 2 × 10^{21} Mx deg^{1} width. The error bars denote the standard error of the mean. 

Open with DEXTER 
Fig. 6 Excess azimuthal flux density above the azimuthal mean underlying each emergence in cycle 23 () versus the area of the emergence (A_{i}). The emergences are represented by green points. The black line represents the average of over all emergences inside area bins of width 100 μHem. The error bars denote the standard error of the mean. 

Open with DEXTER 
Fig. 7 Emerged magnetic flux (Φ_{i}) versus underlying azimuthal flux density for the active regions recorded over cycle 23. Each green point represents an emergence. The stripes in the background indicate azimuthal flux density ranges. The area of the emergences is represented on the righthandside vertical axis. 

Open with DEXTER 
Fig. 8 Probability of emergence per unit time, per unit area and per flux ratio bin (P) as a function of the ratio between the emerged flux and the azimuthal flux underlying the emergence site (r). The different colours correspond to the ranges of azimuthal flux density represented in Fig. 7. The numbers in the legend refer to the mid points of the azimuthal flux ranges. The stripes in the background indicate flux ratio bins of width 0.1. 

Open with DEXTER 
Fig. 9 Ratio of the emerged flux to the underlying azimuthal flux (r_{i}) versus the azimuthal flux density beneath the emergence location for the active regions recorded during cycle 23. Each green point represents an emergence. The stripes in the background indicate azimuthal flux ranges. The two continuous curves separate ephemeral regions (below the lower curve), mediumsized regions (between the two curves) and large active regions (above the upper curve) by their flux content, according to Table 5.1 of Schrijver & Zwaan (2008). These fluxes have been converted from fluxes as given by Kitt Peak magnetograms to fluxes as given by SOLIS magnetograms by using the crosscalibration constants in MuñozJaramillo et al. (2015). 

Open with DEXTER 
Using the emergences shown in Fig. 7 we estimate the probability of emergence as a function of the ratio between the flux content of the emerged active region and the azimuthal flux available within one degree colatitude directly beneath it, (24)To do so, we bin the data according to the ratio r_{i} (with bins of size 0.1) and the underlying azimuthal flux density, b_{φ,i} (with bins of size 2.8 × 10^{21} Mx deg^{1}). We thus obtain the number of emergences in each (r_{i},b_{φ,i}) bin. These are converted to a probability of emergence per unit area and unit time as a function of r by dividing the number of emergences in each b_{φ,i} bin by the cycleaveraged area of the subsurface domain covered by the corresponding azimuthal flux density range, and the duration of cycle 23. The resulting probability distributions are plotted in Fig. 8. Each coloured line corresponds to a different azimuthal flux density range. The distributions decrease rapidly for flux ratios greater than 0.4, suggesting that emergences whose flux comprises more than 40% of the azimuthal flux available underneath the emergence site are rare events.
The probability distributions shown in Fig. 8 seem to converge as we consider stronger azimuthal flux ranges. For the upper end of azimuthal flux ranges the probability of emergence is very similar. The lower probabilities obtained for emergences with smaller underlying azimuthal fluxes are due to a detection bias. To illustrate this, we plot the flux ratio of the emergences as a function of the underlying azimuthal flux density (Fig. 9). The two curves separate ephemeral, medium and large active regions. Active regions lying closer to the lower curve have a lifetime of days, while the lifetime of those closer to the upper curve approaches weeks. Many smaller active regions will not appear in the USAF/NOAA sunspot catalogue, either because they emerge and decay on the backside of the Sun or because they do not have enough flux to form spots or pores. Thus, the probability distributions corresponding to lower azimuthal density fluxes (indicated by dashed lines in Fig. 8) are substantially affected by this detection bias. The fact that the less affected distributions (corresponding to larger amounts of underlying azimuthal flux) seem to converge suggests that the probability of emergence is a function of the ratio of the emerged flux and the azimuthal flux underlying the emergence site.
4. Summary and conclusion
We have provided a nonaxisymmetric model of the magnetic flux transport in the Sun, based on the BabcockLeighton dynamo framework. Using synoptic magnetograms as an input, we inferred the latitudinal and longitudinal distribution of azimuthal flux (per unit colatitude) and its evolution over three cycles.
We calibrated our model by requiring that the azimuthal flux in Cycle 22 in our simulations lied mainly radially underneath the activity belts. This led to a return meridional flow (and/or latitudinal pumping) having an amplitude of u_{0} = 3 m s^{1} and an effective diffusivity for the azimuthal field of η_{0} = 100 km^{2} s^{1}. These values are in the range found by Cameron & Schüssler (2017) for the operation of the solar dynamo.
The azimuthal flux system is highly axisymmetric and antisymmetric about the equator. The departures from axisymmetry represent, on average, approximately 3% of the azimuthal flux at a given location. We found that the nonaxisymmetric structure does not have a significant impact on the location of the emergences or their observed properties. We also found that the probability of emergence is a function of the ratio of the flux content of the emerged active region and the underlying azimuthal flux.
Acknowledgments
D.M.B. acknowledges postgraduate fellowship of the International Max Planck Research School on Physical Processes in the Solar System and Beyond. This work utilizes SOLIS data obtained by the NSO Integrated Synoptic Program (NISP), managed by the National Solar Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under a cooperative agreement with the National Science Foundation. This work was carried out in the context of Deutsche Forschungsgemeinschaft SFB 963 “Astrophysical Flow Instabilities and Turbulence” (Project A16). We thank Manfred Schüssler for his valuable suggestions and his thorough revision of this manuscript.
References
 Babcock, H. W. 1961, ApJ, 133, 572 [NASA ADS] [CrossRef] [Google Scholar]
 Barekat, A., Schou, J., & Gizon, L. 2014, A&A, 570, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., & Schüssler, M. 2015, Science, 347, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., & Schüssler, M. 2016, A&A, 591, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cameron, R., & Schüssler, M. 2017, A&A, 599, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charbonneau, P. 2010, Liv. Rev. Sol. Phys., 7, 91 [Google Scholar]
 ChristensenDalsgaard, J., & Schou, J. 1988, in Seismology of the Sun and SunLike Stars, ed. E. J. Rolfe, ESA SP, 286 [Google Scholar]
 DasiEspuig, M., Solanki, S. K., Krivova, N. A., Cameron, R., & Peñuela, T. 2010, A&A, 518, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DeVore, C. R., Boris, J. P., & Sheeley, Jr., N. R. 1984, Sol. Phys., 92, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Durney, B. R. 1995, Sol. Phys., 160, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2009, Liv. Rev. Sol. Phys., 6, 4 [Google Scholar]
 Guerrero, G., & de Gouveia Dal Pino, E. M. 2008, A&A, 485, 267 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hathaway, D. H., & Rightmire, L. 2011, ApJ, 729, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Jafarzadeh, S., Cameron, R. H., Solanki, S. K., et al. 2014, A&A, 563, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2011, Astron. Lett., 37, 656 [NASA ADS] [CrossRef] [Google Scholar]
 Leighton, R. B. 1964, ApJ, 140, 1547 [NASA ADS] [CrossRef] [Google Scholar]
 Leighton, R. B. 1969, ApJ, 156, 1 [NASA ADS] [CrossRef] [Google Scholar]
 MartinBelda, D., & Cameron, R. H. 2016, A&A, 586, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MuñozJaramillo, A., Senkpeil, R. R., Windmueller, J. C., et al. 2015, ApJ, 800, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, N. J., Brown, B. P., Sacha Brun, A., Miesch, M. S., & Toomre, J. 2014, Sol. Phys., 289, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Martin, S. F. 1990, Sol. Phys., 129, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Zwaan, C. 2008, Solar and Stellar Magnetic Activity (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Spruit, H. C. 2011, Theories of the Solar Cycle: A Critical View, eds. M. P. Miralles, & J. Sánchez Almeida, 39 [Google Scholar]
 Tóth, G., van der Holst, B., & Huang, Z. 2011, ApJ, 732, 102 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A. A., Cartledge, N. P., & Priest, E. R. 1998, ApJ, 501, 866 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M., Sheeley, Jr., N. R., & Nash, A. G. 1991, ApJ, 383, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y.M., Robbrecht, E., & Sheeley, Jr., N. R. 2009, ApJ, 707, 1372 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: The radial integral of the φ component of the induction equation
We derive the evolution equation for the azimuthal flux per unit colatitude, (A.1)under the assumptions specified in Sect. 2.1.
The φ component of the induction equation in spherical coordinates, (r,θ,φ), can be written as follows: (A.2)where we have assumed η = η(r).
The first row of Eq. (A.2)includes the advection and shear terms. Integrating the shear terms radially, we obtain (A.3)In the above equation, B_{r}  _{Rb} vanishes since we assume that the magnetic field does not penetrate the radiative interior. With B_{θ} = 0 in the NSSL, we can change the upper limit of integration of the integral in Eq. (A.3)to R_{NSSL}. The lower limit can be changed to the top of the tachocline, R_{T}, since there is no poloidal field in the convection zone part of the tachocline. In the new integration domain, Ω depends only on θ, so we can move it outside the integral. This yields (A.4)where we have defined (A.5)The radial integral of the radial advection term in Eq. (A.2), , vanishes since both u_{r} and B_{φ} vanish at R_{b} and R_{⊙}. Integrating the latitudinal advection term yields (A.6)where is a weighted average of the meridional flow over the convection zone, (A.7)The diffusion term in Eq. (A.2)reads: (A.8)The integral of the radial part, , vanishes since there is no diffusive flux transport across the boundaries. Following results of Cameron & Schüssler (2016), we assume an effective diffusivity of azimuthal flux (A.9)Substituting the above expression in Eq. (A.8)and integrating the rest of the diffusion terms leads to: (A.10)Combining Eqs. (A.4), (A.6), and (A.10), and introducing the source term that is necessary to ensure connectivity with the surface sources, we obtain the final form of the equation: (A.11)To simplify the notation and keep consistency with Eq. (3), we write B_{r} instead of B_{r}  _{R⊙} to refer to the radial field at the surface.
Appendix B: Calculation of b_{θ}
The quantity b_{θ} can be calculated in terms of b_{φ} and B_{r}  _{R⊙} from the solenoidality condition, ∇·B = 0. Writing the divergence operator in spherical coordinates leads to (B.1)Integrating over the convection zone and using the definitions (A.1)and (A.5)we obtain (B.2)where we have used B_{r}  _{Rb} = 0 and B_{θ}  _{Rb<r<RT} = B_{θ}  _{RNSSL<r<R⊙} = 0. Integrating now in θ yields (B.3)where, again, B_{r} denotes now the surface field.
All Figures
Fig. 1 Initial condition of the simulations. Top left: synoptic magnetogram corresponding to CR 1625. Top right: synoptic magnetogram corresponding to CR 1625, remapped to an equiangular grid and resampled to the highest angular degree order, l, used to compute the potential field extrapolation. Bottom left: radial integral of the azimuthal field across the convection zone, b_{φ}, extrapolated from the surface field. Bottom right: radial integral of the θ component of the magnetic field across the convection zone, b_{θ}, extrapolated from the surface field. Red and blue indicate opposite polarities in all the maps. 

Open with DEXTER  
In the text 
Fig. 2 Number of emergences (N) as a function of the normalized underlying azimuthal flux density (). The different curves correspond to different choices of the free parameters. Continuous black line: u_{0} = 3 m s^{1} and η_{0} = 100 km^{2} s^{1}. Dashed blue line: u_{0} = 5 m s^{1} and η_{0} = 400 km^{2} s^{1}. Dashed red line: u_{0} = 1 m s^{1} and η_{0} = 25 km^{2} s^{1}. Continuous green line: u_{0} = 6 m s^{1} and η_{0} = 600 km^{2} s^{1}. 

Open with DEXTER  
In the text 
Fig. 3 Distribution of the quantity as a function of the parameters u_{0} and η_{0}. Each circle represents a run. The size of each circle is proportional to the quantity ξ, which measures the adjustment of the run to our requirement that the simulated azimuthal flux lies underneath the observed active regions. For better visualization, we also encode the value of ξ in the colour of the circles. 

Open with DEXTER  
In the text 
Fig. 4 Surface field (B_{r}), azimuthal flux density (b_{φ}) and its nonaxisymmetric component (), for CR 1987 (close to the middle of cycle 23) and CR 2060 (near the end of that cycle). Red and blue represent opposite polarities. The squares on the top and middle rows represent observed emergence sites, extracted from the USAF/NOAA sunspot group database. 

Open with DEXTER  
In the text 
Fig. 5 Nonaxisymmetric component of the azimuthal flux underlying each emergence in cycle 23 () versus its axisymmetric component (⟨ b_{φ,i} ⟩). The emergences are represented by green points. The black line represents the average of over all emergences inside bins of 2 × 10^{21} Mx deg^{1} width. The error bars denote the standard error of the mean. 

Open with DEXTER  
In the text 
Fig. 6 Excess azimuthal flux density above the azimuthal mean underlying each emergence in cycle 23 () versus the area of the emergence (A_{i}). The emergences are represented by green points. The black line represents the average of over all emergences inside area bins of width 100 μHem. The error bars denote the standard error of the mean. 

Open with DEXTER  
In the text 
Fig. 7 Emerged magnetic flux (Φ_{i}) versus underlying azimuthal flux density for the active regions recorded over cycle 23. Each green point represents an emergence. The stripes in the background indicate azimuthal flux density ranges. The area of the emergences is represented on the righthandside vertical axis. 

Open with DEXTER  
In the text 
Fig. 8 Probability of emergence per unit time, per unit area and per flux ratio bin (P) as a function of the ratio between the emerged flux and the azimuthal flux underlying the emergence site (r). The different colours correspond to the ranges of azimuthal flux density represented in Fig. 7. The numbers in the legend refer to the mid points of the azimuthal flux ranges. The stripes in the background indicate flux ratio bins of width 0.1. 

Open with DEXTER  
In the text 
Fig. 9 Ratio of the emerged flux to the underlying azimuthal flux (r_{i}) versus the azimuthal flux density beneath the emergence location for the active regions recorded during cycle 23. Each green point represents an emergence. The stripes in the background indicate azimuthal flux ranges. The two continuous curves separate ephemeral regions (below the lower curve), mediumsized regions (between the two curves) and large active regions (above the upper curve) by their flux content, according to Table 5.1 of Schrijver & Zwaan (2008). These fluxes have been converted from fluxes as given by Kitt Peak magnetograms to fluxes as given by SOLIS magnetograms by using the crosscalibration constants in MuñozJaramillo et al. (2015). 

Open with DEXTER  
In the text 