Evolution of the Sun's non-axisymmetric toroidal field

TL;DR

This study models the Sun's non-axisymmetric azimuthal magnetic flux using a kinematic flux transport model, revealing that non-axisymmetric components are minor and have limited impact on flux emergence locations.

Contribution

It introduces a new kinematic flux transport model based on the Babcock-Leighton framework to infer subsurface non-axisymmetric magnetic flux from surface observations.

Findings

Non-axisymmetric azimuthal flux constitutes about 3% of total flux.

Non-axisymmetric structures have limited influence on flux emergence.

Emergence probability depends on the ratio of active region flux to underlying azimuthal flux.

Abstract

Aims. We aim to infer the sub-surface distribution of the Sun's non-axisymmetric 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 Babcock-Leighton 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 constrain 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 non-axisymmetric 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.

Figures

Equations

$$B_j=⟨B_j⟩+B_j^{\prime },$$

$$⟨B_j⟩(\theta ,\varphi )=\frac{{\int }_{\theta -\delta }^{\theta +\delta }{\int }_{\varphi -\delta }^{\varphi +\delta }{B}_j\sin {\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }\mathrm{d}{\theta }^{\prime }}{{\int }_{\theta -\delta }^{\theta +\delta }{\int }_{\varphi -\delta }^{\varphi +\delta }\sin {\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }\mathrm{d}{\theta }^{\prime }},$$

$$\frac{\partial B_r}{\partial t}=$$

$$+\frac{{\eta }_H}{R_{\odot }^2}\left[\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial B_r}{\partial \theta }\right)+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2{B}_r}{\partial {\varphi }^2}\right]$$

$$+S(\theta ,\varphi ,t),$$

$${\Omega }_{R_{\odot }}(\theta )=14.437-1.48{\cos }^2\theta -2.99{\cos }^4\theta {[}^{\circ }/day].$$

$$u_M=-11\cos (2\theta )[\mathrm{m}{\mathrm{s}}^{-1}].$$

$$b_{\varphi }={\int }_{R_{\mathrm{b}}}^{R_{\odot }}{B}_{\varphi }r\mathrm{d}r.$$

$$\frac{\partial b_{\varphi }}{\partial t}$$

$$+b_{\theta }\sin \theta {\frac{\mathrm{d}\Omega }{\mathrm{d}\theta }|}_{R_{\mathrm{NSSL}}}+\frac{\partial (b_{\theta }\sin \theta )}{\partial \theta }{\Omega }_{R_{\mathrm{NSSL}}}$$

$$+\frac{{\eta }_0}{R_{\odot }^2}[\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial b_{\varphi }}{\partial \theta }\right)+2\frac{\cos \theta }{{\sin }^2\theta }\frac{\partial b_{\theta }}{\partial \varphi }-\frac{b_{\varphi }}{{\sin }^2\theta }$$

$$+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2{b}_{\varphi }}{\partial {\varphi }^2}]-\frac{1}{R_{\odot }}\frac{\partial (\bar{u}{b}_{\varphi })}{\partial \theta }+S_{\varphi }(\theta ,\varphi ,t),$$

$$b_{\theta }={\int }_{R_{\mathrm{T}}}^{R_{\mathrm{NSSL}}}{B}_{\theta }r\mathrm{d}r.$$

$${\Omega }_{R_{\mathrm{NSSL}}}(\theta )={\Omega }_{R_{\odot }}(\theta )+0.53{[}^{\circ }/day].$$

$$\eta (r)={\eta }_0\frac{r^2}{R_{\odot }^2},$$

$$\bar{u}=u_0\cos (2\theta )[\mathrm{m}{\mathrm{s}}^{-1}],$$

$$\Delta B_r(\theta ,\varphi ,t^n)=B_r^n-{\tilde{B}}_r^n,$$

$$\Delta B_r(\theta ,\varphi ,t^n)={\int }_{t^{n-1}}^{t^n}S(\theta ,\varphi ,t)\mathrm{d}t.$$

$$\Delta b_{\varphi }(\theta ,\varphi ,t^n)={\int }_{t^{n-1}}^{t^n}{S}_{\varphi }(\theta ,\varphi ,t)\mathrm{d}t.$$

$$B_r(\theta ,\varphi ,t^n)=\sum _m{a}_m(\theta ,t^n)e^{im\varphi }.$$

$${\frac{\partial a_0(\theta ,t^n)}{\partial \theta }|}_{\theta =0,\pi }=0;$$

$$a_{m\ne 0}(\theta =0,\pi ;t^n)=0.$$

$$a_m({\theta }^{\ast },t^n)=a_m({\theta }^{-},t^n)+\frac{a_m({\theta }^{+},t^n)-a_m({\theta }^{-},t^n)}{{\theta }^{+}-{\theta }^{-}}({\theta }^{\ast }-{\theta }^{-}),$$

$$a_0({\theta }^{\ast },t^n)=a_0({\theta }_0^{+},t^n);$$

$$a_{m\ne 0}({\theta }^{\ast },t^n)=\frac{a_{m\ne 0}({\theta }_0^{+},t^n)}{{\theta }_0^{+}}{\theta }^{\ast }.$$

$$a_0({\theta }^{\ast },t^n)=a_0({\theta }_0^{-},t^n);$$

$$a_{m\ne 0}({\theta }^{\ast },t^n)=a_{m\ne 0}({\theta }_0^{-},t^n)-\frac{a_{m\ne 0}({\theta }_0^{-},t^n)}{\pi -{\theta }_0^{-}}({\theta }^{\ast }-{\theta }_0^{-}).$$

$$\xi =\frac{{\sum }_c{N}_c|b_{\varphi ,c}|}{{\sum }_c{N}_c/|b_{\varphi ,c}|},$$

$${\bar{b}}_{\varphi ,i}^{\prime }(⟨b_{\varphi ,i}⟩<0)$$

$${\bar{b}}_{\varphi ,i}^{\prime }(⟨b_{\varphi ,i}⟩>0)$$

$$r_i=\frac{{\Phi }_i}{2{\int }_{1\mathrm{deg}}{b}_{\varphi ,i}\mathrm{d}\theta }.$$

$$b_{\varphi }={\int }_{R_{\mathrm{b}}}^{R_{\odot }}{B}_{\varphi }r\mathrm{d}r,$$

$$r\frac{\partial B_{\varphi }}{\partial t}$$

$$-\frac{1}{\sin \theta }\frac{\partial \eta }{\partial r}\frac{\partial B_r}{\partial \varphi }+\frac{\partial }{\partial r}\left[\eta \frac{\partial }{\partial r}(r{B}_{\varphi })\right]$$

$$+\frac{\eta }{r}[\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial B_{\varphi }}{\partial \theta }\right)+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2{B}_{\varphi }}{\partial {\varphi }^2}$$

$$+\frac{2}{\sin \theta }\frac{\partial B_r}{\partial \varphi }+\frac{2\cos \theta }{{\sin }^2\theta }\frac{\partial B_{\theta }}{\partial \varphi }-\frac{B_{\varphi }}{{\sin }^2\theta }],$$

Full text

*Methods. *We have built a kinematic flux transport model of the solar dynamo based on the Babcock-Leighton 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 constrain 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, $\sim 3\%$ of the total azimuthal flux. Owing to its relative weakness, the non-axisymmetric 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.

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 Babcock-Leighton models.

In Babcock-Leighton 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 cross-equatorial transport of preceding polarity flux, and leads to the cancellation of the polar fields and the build-up of a new, reversed axial dipole, which is the source of azimuthal field for the new activity cycle.

Babcock-Leighton models have gained substantial support in recent years. Dasi-Espuig et al. (2010) found a strong correlation between the strength of a Babcock-Leighton 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 build-up 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 Babcock-Leighton 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 sub-surface 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 Babcock-Leighton 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,

[TABLE]

where $j=\{r,\theta,\phi\}$,

[TABLE]

and $\delta$ is the scale over which the average is performed. Using $\delta\sim 2-3^{\circ}$ 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 $j$-th component of the averaged magnetic field by $B_{j}$.

Our model consists of two two-dimensional domains in the $(\phi,\theta)$ 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):

[TABLE]

where $\Omega_{R_{\odot}}$ is the differential rotation, $u_{M}$ is the meridional flow and $\eta_{H}$ is the surface diffusivity associated to the convective flows. The emergence of new flux is described by the source term $S(\theta,\phi,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):

[TABLE]

The second term in Eq. (3) corresponds to the surface meridional flow. Following van Ballegooijen et al. (1998), we model the meridional flow as:

[TABLE]

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; Martin-Belda & Cameron 2016). We use a surface diffusivity of $\eta_{H}=250\,\mathrm{km^{2}s^{-1}}$, as indicated by observations (Schrijver & Martin 1990; Jafarzadeh et al. 2014).

The source term $S(\theta,\phi,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:

[TABLE]

In the above expression, $R_{\rm b}$ refers to the bottom of the convection zone and $R_{\odot}$ 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:

The magnetic field is purely radial in the near-surface shear layer (NSSL), owing to strong downwards turbulent pumping. 2. 2.

The magnetic field does not penetrate the radiative interior. 3. 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. 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 (Christensen-Dalsgaard & Schou 1988).

To derive an evolution equation for $b_{\phi}$, we integrated the azimuthal component of the induction equation in radius. The resulting equation reads (see Appendix A):

[TABLE]

where $\Omega_{R_{\rm NSSL}}$ is the differential rotation at the bottom of the NSSL, $\bar{u}$ is an effective return meridional flow, $\eta_{0}$ is the effective diffusivity of the azimuthal field, $S_{\phi}(\theta,\phi,t)$ is a source term associated to flux emergence, and

[TABLE]

In the above expression, $R_{\rm T}$ refers to the top of the tachocline and $R_{\rm NSSL}$ refers to the bottom of the NSSL. The quantity $b_{\theta}$ is calculated in terms of $b_{\phi}$ 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 $\theta$-component of the magnetic field, $b_{\theta}$, can be obtained in terms of $b_{\phi}$ and $B_{r}$ from the solenoidality condition $\nabla\cdot\mathbf{B}=\mathbf{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

[TABLE]

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):

[TABLE]

where $\eta_{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 $\eta_{0}\sim 150-450\,\mathrm{km^{2}s^{-1}}$.

The fifth term corresponds to the advection of the azimuthal flux by an effective equatorward flow, which we modeled as

[TABLE]

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 sub-surface 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_{\phi}(\theta,\phi,t)$, we extrapolated the surface sources downward via a potential field solution (see Sect. 2.2).

Our model, therefore, consists of: (a) a two-dimensional domain representing the surface of the Sun, in which the (radial) surface field evolves according to Eq. (3); (b) a two-dimensional 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 $B_{r}^{n}$ 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

[TABLE]

where $\tilde{B}_{r}^{n}$ 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(\theta,\phi,t)$ through

[TABLE]

In practice, we add $\Delta 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 $\nabla\cdot\mathbf{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 $\Delta B_{r}(\theta,\phi,t^{n})$. Let $\Delta B_{\phi}(\theta,\phi,t^{n})$ be the azimuthal component of the extrapolated field, and $\Delta b_{\phi}=\int_{R_{b}}^{R_{\odot}}\Delta B_{\phi}r\,\mathrm{d}r$. The source term in Eq. (7), $S_{\phi}(\theta,\phi,t^{n})$, is related to $\Delta b_{\phi}(\theta,\phi,t^{n})$ through

[TABLE]

The extrapolated $\Delta b_{\phi}$ is added to the simulated $b_{\phi}$ every rotation, at the same time $\Delta B_{r}$ is added to the surface field.

The potential field extrapolation of $\Delta B_{r}(\theta,\phi,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 $\theta$, rather than in $\cos\theta$. 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

[TABLE]

At the poles, regularity of $B_{r}$ translates into the following boundary conditions for $a_{m}$:

[TABLE]

The value of the Fourier coefficient at point of colatitude $\theta^{*}$ of the new grid is given by:

[TABLE]

where $\theta^{-}$ and $\theta^{+}$ are the points of the old grid adjacent to $\theta^{*}$ in the $\theta$ direction. For the points of the new grid that are located between the north pole and the northernmost grid point of the old grid, $\theta_{0}^{+}$, Eqs. (16) and (17) lead to:

[TABLE]

For grid points between the southernmost point of the old grid, $\theta_{0}^{-}$, and the south pole, we have:

[TABLE]

Since the maximum angular order of the spherical harmonic analysis is limited by an anti-aliasing condition, the decomposition of $B_{r}$ in spherical harmonics has the effect of slightly smoothing the magnetograms (see the top row of Fig. 1).

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_{\phi}$ and $b_{\theta}$. The sign of $b_{\phi}$ indicates the direction of the integrated azimuthal field. In our chosen coordinate system, a positive value of $b_{\phi}$ corresponds to azimuthal field pointing in the sense of rotation of the Sun. The sign of $b_{\theta}$ is predominantly negative, which reflects the sign of the axial dipole at this cycle minimum. The irregularities of $b_{\theta}$ 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\,\mathrm{m\,s^{-1}}$, and the diffusivity $\eta_{0}$ was varied between $25$ and $600\,\mathrm{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_{\phi,\,i}=b_{\phi}(\lambda_{i},\phi_{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 non-axisymmetric 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 $\eta_{0}$ give rise to different global magnitudes of the azimuthal flux density so, for easier comparison, we normalized $b_{\phi,\,i}$ to its maximum value in each simulation. The emergences were grouped in bins of width $4\cdot 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\,\mathrm{m\,s^{-1}}$ and $\eta_{0}=100\,\mathrm{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 two-peak 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 well-separated peaks, we considered the quantity

[TABLE]

where $c$ runs through the $b_{\phi,\,i}$ bins and $b_{\phi,\,c}$ is the mid-point of each bin. The value of $\xi$ will be bigger for the simulations where the emergences occur farther away from the places where $b_{\phi}\sim 0$. Figure 3 shows the value of $\xi$ for all test runs carried out. We find that the combination $u=3\,\mathrm{m\,s^{-1}}$ and $\eta_{0}=100\,\mathrm{km^{2}s^{-1}}$ maximizes $\xi$. These parameters are close to the range found by Cameron & Schüssler (2017) for the operation of the solar dynamo in an updated Leighton 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.

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 sub-surface 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 low-resolution 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 $\sim 15^{\circ}$ 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 winding-up of the reversed poloidal field.

The bottom row of Fig. 4 shows the non-axisymmetric part of the integrated azimuthal field, calculated as $b_{\phi}^{\prime}=b_{\phi}-\langle b_{\phi}\rangle$, where $\langle b_{\phi}\rangle$ is the azimuthal mean of $b_{\phi}$. The magnitude of $b_{\phi}^{\prime}$ represents, on average, $\sim 3\%$ of the total azimuthal flux density. This non-axisymmetric structure arises from the emergence process (which is the only non-axisymmetric 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 non-axisymmetric structure on the emergence process

In order to investigate whether the non-axisymmetric 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, $\langle b_{\phi,\,i}\rangle$. 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 $\langle b_{\phi,\,i}\rangle<0$ are mainly located in the north hemisphere, while those with $\langle b_{\phi,\,i}\rangle>0$ correspond to the south hemisphere.

An influence of the non-axisymmetric structure of the azimuthal flux on the triggering of the emergence process would lead to a non-zero average value of $b_{\phi,\,i}^{\prime}$ 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 $b_{\phi,\,i}^{\prime}$ over all the emergences for which $\langle b_{\phi,\,i}\rangle>0$ would yield a positive value. In the other hemisphere, the average $b_{\phi,\,i}^{\prime}$ would be smaller than zero. Computing these averages yields

[TABLE]

The average value of $b_{\phi,\,i}^{\prime}$ 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 sub-surface azimuthal flux.

Next, we consider the possible influence of the non-axisymmetric structure on the active region areas. Figure 6 shows the excess of azimuthal flux density above the azimuthal average beneath each active region, $b_{\phi,\,i}^{\prime}/{\rm sign}(\langle b_{\phi,\,i}\rangle)$, versus the active region area, $A_{i}$. Again, there is no significant deviation from zero, which suggests that the inferred non-axisymmetric structure of the sub-surface 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, $\Phi_{i}$, versus the unsigned underlying azimuthal flux density, $|b_{\phi,\,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 cross-calibration factors in Table 2 of Muñoz-Jaramillo et al. (2015). The resulting relationship is $\Phi_{i}\,\mathrm{[Mx]}=1.44\cdot 10^{19}A_{i}\,\mathrm{[\mu 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 $\Omega$-shaped magnetic structure that crosses the solar surface twice.

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,

[TABLE]

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_{\phi,\,i}$ (with bins of size $2.8\cdot 10^{21}\,{\rm Mx\,deg^{-1}}$). We thus obtain the number of emergences in each $(r_{i},b_{\phi,\,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_{\phi,\,i}$ bin by the cycle-averaged 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 non-axisymmetric model of the magnetic flux transport in the Sun, based on the Babcock-Leighton 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\,\mathrm{ms^{-1}}$ and an effective diffusivity for the azimuthal field of $\eta_{0}=100\,\mathrm{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 non-axisymmetric 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.

Aknowledgements

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.

Appendix A The radial integral of the $\phi$ component of the induction equation

We derive the evolution equation for the azimuthal flux per unit colatitude,

[TABLE]

under the assumptions specified in Sect. 2.1.

The $\phi$ component of the induction equation in spherical coordinates, $(r,\theta,\phi)$, can be written as follows:

[TABLE]

where we have assumed $\eta=\eta(r)$.

The first row of Eq. (26) includes the advection and shear terms. Integrating the shear terms radially, we obtain

[TABLE]

In the above equation, $B_{r}|_{R_{\rm b}}$ vanishes since we assume that the magnetic field does not penetrate the radiative interior. With $B_{\theta}=0$ in the NSSL, we can change the upper limit of integration of the integral in Eq. (27) to $R_{\rm NSSL}$. The lower limit can be changed to the top of the tachocline, $R_{\rm T}$, since there is no poloidal field in the convection zone part of the tachocline. In the new integration domain, $\Omega$ depends only on $\theta$, so we can move it outside of the integral. This yields

[TABLE]

where we have defined

[TABLE]

The radial integral of the radial advection term in Eq. (26), $\int_{R_{\rm b}}^{R_{\odot}}\partial(u_{r}B_{\phi})/\partial r\,\mathrm{d}r$, vanishes since both $u_{r}$ and $B_{\phi}$ vanish at $R_{\rm b}$ and $R_{\odot}$. Integrating the latitudinal advection term yields

[TABLE]

where $\bar{u}$ is a weighted average of the meridional flow over the convection zone,

[TABLE]

The diffusion term in Eq. (26) reads:

[TABLE]

The integral of the radial part, $\int_{R_{\rm b}}^{R_{\odot}}\frac{\partial}{\partial r}\left[\eta\frac{\partial}{\partial r}(rB_{\phi})\right]\,\mathrm{d}r$, 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

[TABLE]

Substituting the above expression in Eq. (32) and integrating the rest of the diffusion terms leads to:

[TABLE]

Combining Eqs. (28), (30), and (34), and introducing the source term that is necessary to ensure connectivity with the surface sources, we obtain the final form of the equation:

[TABLE]

To simplify the notation and keep consistency with Eq. (3), we write $B_{r}$ instead of $B_{r}|_{R_{\odot}}$ to refer to the radial field at the surface.

Appendix B Calculation of $b_{\theta}$

The quantity $b_{\theta}$ can be calculated in terms of $b_{\phi}$ and $B_{r}|_{R_{\odot}}$ from the solenoidality condition, $\nabla\cdot\mathbf{B}=0$. Writing the divergence operator in spherical coordinates leads to

[TABLE]

Integrating over the convection zone and using the definitions (25) and (29) we obtain

[TABLE]

where we have used $B_{r}|_{R_{\rm b}}=0$ and $B_{\theta}|_{R_{\rm b}<r<R_{\rm T}}=B_{\theta}|_{R_{\rm NSSL}<r<R_{\odot}}=0$. Integrating now in $\theta$ yields

[TABLE]

where, again, $B_{r}$ denotes now the surface field.