1 Introduction

Magnetic fields play a key role in the existence and variability of a wide variety of phenomena on the Sun. These range from relatively stable, slowly evolving features such as sunspots (Solanki, 2003), coronal loops (Reale, 2010), and solar prominences (Labrosse et al., 2010; Mackay et al., 2010), to highly dynamic phenomena such as solar flares (Benz, 2008) and coronal mass ejections (Forbes et al., 2006; Chen, 2011). Solar magnetic fields may directly or indirectly affect the Earth through the Sun’s open magnetic flux (Balogh et al., 1995), solar wind (Hollweg, 2008), and total solar irradiance variations (Fröhlich and Lean, 1998).

Our present day understanding of solar magnetic fields dates back to 1908 when G.E. Hale made the first magnetic field observations of sunspots (Hale, 1908). However, it was not until the systematic mapping of the Sun’s magnetic field carried out by the Babcocks (Babcock and Babcock, 1952; Babcock and Livingston, 1958; Babcock, 1959) that the true nature of solar magnetic activity became apparent. While significant advances in observations have been made over the last 50 years, only the strength and distribution of the line-of-sight component at the level of the photosphere has been regularly measured over solar cycle time-scales. However, over the last 5 years, vector magnetic field measurements at the photospheric level have been systematically carried out by the satellite missions of Hinode since 2006 (Kosugi et al., 2007) and SDO (Solar Dynamics Observatory) since 2010 (Pesnell et al., 2012). In addition, systematic ground based measurements of vector magnetic fields have been made by SOLIS (Synoptic Optical Long-term Investigations of the Sun) since 2009 (Keller et al., 2003). While Hinode only observes vector magnetic fields over localised areas, the global capabilities of SDO and SOLIS provide us with a unique opportunity to observe the large-scale spatial distribution of vector magnetic fields across the solar surface. Over time scales of years, SDO and SOLIS will significantly enhance our understanding of the emergence and transport of magnetic fields at the level of the photosphere and their subsequent effect on the global corona. While magnetic fields may be observed directly at the photospheric level, due to low densities the same is not true for the corona. As many important phenomena occur in the solar corona, a key component in our understanding of solar magnetic fields - and the build up of energy within them - is the use of theoretical models to construct (or extrapolate) coronal magnetic fields from photospheric data.

While solar magnetic fields have been observed in detail over long periods of time, the same is not true for other stars. In recent years however, the technique of Zeeman Doppler Imaging (Semel, 1989) has lead to a significant advance in our understanding of magnetic fields on other stars. Results show a wide range of magnetic distributions across stars of varying mass and spectral type (Donati et al., 2009). With the accurate measurement of stellar magnetic fields, techniques developed to model solar magnetic fields, both at photospheric and coronal levels, are now widely applied in the stellar context (Cohen et al., 2010; Jardine et al., 2011).

In this review, we primarily focus on our present day understanding of global solar magnetic fields from both observations and theoretical models. The reader should note that we focus on quasi-static or steady-state coronal models and will not consider fully dynamic eruptive models. Within the review we also restrict ourselves to global aspects of the Sun’s magnetic field so will consider neither the Sun’s small-scale field nor limited field-of-view models. The review will focus in particular on advances made in the last 15 years, although we have not aimed at completeness in material or references. Additional topics will be added in future revisions. Where appropriate, we expand this discussion into stellar magnetic fields to summarise new results or to describe the application of models developed for the Sun in the stellar context. The review is split into three distinct parts, where each part is largely self-contained. Thus the reader may focus on each part separately if desired. The review is split in the following way:

  • Global Photospheric Magnetic Fields (Section 2)

    Observations of global solar and stellar photospheric magnetic fields are briefly discussed (Section 2.1). Following from this, magnetic flux transport models used to simulate the spatial and temporal evolution of photospheric magnetic fields are described, along with the variety of extensions and applications of these models (Section 2.2). For a historical discussion of magnetic flux transport models, the reader is referred to the article by Sheeley Jr (2005).

  • Global Coronal Models (Section 3)

    This section of the review surveys the wide variety of techniques used to model global coronal magnetic fields. This includes the various approximations that are applied from static extrapolation techniques, to time dependent quasi-static models and finally recent advances in global non-eruptive MHD models.

  • Application of Global Models (Section 4)

    The final part of the review considers the combined application of both magnetic flux transport models and coronal modeling techniques to model a variety of phenomena found on the Sun. These include the Sun’s open magnetic flux (Section 4.1), coronal holes (Section 4.2), and the hemispheric pattern of solar filaments (Section 4.3). In addition, recent advances in MHD models for modeling the plasma emission of the corona are also considered (Section 4.4).

Finally, in Section 5 a brief summary is given and some outstanding problems or areas of advancement are outlined.

2 Photospheric Magnetic Fields

2.1 Observations

Presently, three solar cycles of continuous data have been collected by a variety of groundand space-based observatories (Mount Wilson Observatory, Wilcox Solar Observatory, Kitt Peak, SoHO/MDI, SOLIS) mapping the distribution and evolution of the Sun’s normal magnetic field component at the level of the photosphere. An illustration of this can be seen in Figure 1a (from Hathaway, 2010). The image is known as the solar magnetic “butterfly diagram” and illustrates the longitude-averaged radial magnetic field as a function of time (horizontal axis) versus sine-latitude (vertical axis). Yellow represents positive flux and blue negative flux, where the field saturates at ± 10 G. The main features in the long-term evolution of the global magnetic field are:

  • — At the start of each solar cycle (ca. 1975, 1986, 1996) the majority of magnetic flux sits in the polar regions which are of opposite polarity. New magnetic flux then emerges in the form of sunspots or large magnetic bipoles in two latitude bands between ± 30° (Maunder, 1913). As the cycle progresses these bands approach the equator.

  • — As new flux emerges in the form of magnetic bipoles, the polarities lie mainly in an eastwest orientation. In each hemisphere the leading polarity (in the direction of rotation) has the same sign as the polar field of the hemisphere in which it lies, the following polarity has the opposite sign (Hale’s Polarity Law, Hale and Nicholson, 1925). In addition, the majority of bipoles emerge subject to Joy’s Law, where the leading polarity lies equatorward of the following (Hale et al., 1919; Howard, 1991). The effect of Joy’s Law is clearly seen in Figure 1a as a latitudinal separation of magnetic fluxes: following polarity dominates at latitudes poleward of ± 30°; leading polarity dominates at latitudes closer to the equator.

  • — This latitudinal separation combined with the effects of meridional flow (Duvall Jr, 1979; Hathaway and Rightmire, 2010) and the dispersal of magnetic flux out from complexes of activity, leads to the preferential transport poleward of the following polarity. In contrast, the leading polarity in each hemisphere, which lies at lower latitudes, partially escapes the effect of meridional flow to disperse and cancel across the equator. As a result, in each hemisphere more following flux is transported poleward, where it first cancels the existing polar field and then builds up a new polar field of following polarity. This transport of flux poleward occurs in quasi-episodic poleward streams which extend out of the active latitudes. Such poleward streams are anchored in closely packed clusters of sunspots. The clusters may persist in the same location for several solar rotations as new bipoles emerge to refresh the cluster. These sites are referred to as complexes of activity (Gaizauskas et al., 1983) or as activity nests (Schrijver and Zwaan, 2000). Harvey and Zwaan (1993) found that at least half of the active regions larger than 3.5 deg2 emerge within nests. Nests themselves are sometimes grouped together as interacting units; the grouping format determines how flux streams poleward as the groups decay (Gaizauskas, 2008).

  • — The reversal in sign of the polar field typically occurs 1 – 2 years after cycle maximum. In recently observed cycles this process occurs with every 11-year activity cycle, with a further 11-year activity cycle required before the polar fields in each hemisphere return to their initial sign. A 22-year magnetic cycle therefore overlies the 11-year activity cycle.

Figure 1
figure 1

(a) The Solar Butterfly Diagram (reproduced from Hathaway, 2010). Yellow represents positive flux and blue negative flux where the field saturates at ± 10 G. (b) Example of a typical radial magnetic field distribution for AB Dor taken through Zeeman Doppler Imaging (ZDI, data from Donati et al., 2003) where the image saturates at ± 300 G. White/black denotes positive/negative flux. Due to the tilt angle of AB Dor, measurements can only be made in its northern hemisphere. Image reproduced by permission from Mackay et al. (2004), copyright by RAS.

While continuous global measurements of magnetic activity only exist from around 1975, observations of the numbers of sunspots may be used to provide a long running data set of solar activity back to 1611 (Hoyt and Schatten, 1998; Vaquero et al., 2011). These show that on top of the approximate 11 or 22-year activity cycle there are strong modulations in the number of sunspots (or magnetic flux emergence rate) over periods of centuries. It is also possible for large-scale magnetic activity to disappear. Such an event occurred between 1645 and 1715 where it is known as the Maunder minimum (Eddy, 1976; Ribes and Nesme-Ribes, 1993). Before 1611 indicators of magnetic activity on the Sun may be found through the use of proxies such as 14C (Stuiver and Quay, 1980) and 10Be isotopes (Beer et al., 1990). Through this, reconstructions of the level of magnetic activity over the past 10,000 years may be made (Usoskin, 2008) and show that many such “grand minima” have occurred over the last 10,000 years.

While our present day knowledge of solar magnetic fields is vast, the majority of this knowledge comes from observing the line-of-sight component at the level of the photosphere. To gain a much fuller understanding of the Sun’s magnetic field, vector field measurements are required (Lites et al., 1994). However, these measurements are complicated to make, with problems including low signal-to-noise ratios and resolving the 180° ambiguity. In addition, in the past such measurements have not been regularly made until the launch of Hinode (ca. 2006) which can make vector magnetic field measurements over localised areas of the Sun. However, with the new space mission of SDO (launched in 2010) and the availability of ground based vector magnetic field measurements (SOLIS), such measurements should now be made regularly over the full solar disk in strong field regions. This - combined with modeling techniques - should significantly enhance our understanding of solar magnetic fields in years to come. While current vector magnetic field measurements should increase our knowledge, to fully understand the Sun’s magnetic field vector, measurements are also required in weak field regions over the entire Sun. This poses a significant technical challenge, but is something that future instrument designers must consider.

While our understanding of magnetic fields on stars other than the Sun is at an early stage, significant progress has been made over the last 10 years. For young, rapidly rotating solar-like stars, very different magnetic field distributions may be found compared to the Sun. An example of this can be seen in Figure 1b, where a typical radial magnetic field distribution for AB Dor taken through Zeeman Doppler Imaging (ZDI, Semel, 1989) is shown. AB Dor has a rotation period of around 1/2 day, which is significantly shorter than that of the Sun (27 days). Compared to the Sun, key differences include kilogauss polar fields covering a large area of the pole and the mixing of both positive and negative polarities at the poles. While this is an illustration of a single star at a single time, many such observations have been made across a wide range of spectral classes. In Figure 3 of Donati et al. (2009) the varying form of morphology and strength of the magnetic fields for a number of stars ranging in spectral class from early F to late M can be seen compared to that of the Sun. The plot covers stars with a rotation period ranging from 0.4 to 30 days and masses from (0.09 to 2M). While ZDI magnetic field data sets are generally too short to show cyclic variations, recent observations of the planet-hosting star, τ Bootis, have shown that it may have a magnetic cycle with period of only 2 years (Fares et al., 2009). Indirect evidence for cyclic magnetic field variations on other stars can also be seen from the Mt. Wilson Ca II H+K observations, which use chromospheric observations as a proxy for photospheric magnetic activity (Baliunas et al., 1995). These show that magnetic activity on stars of spectral types G2 to K5V has three main forms of variation. These are (i) moderate activity and regular oscillations similar to the Sun, (ii) high activity and irregular variations (mainly seen on young stars), and finally (iii) stars with flat levels of activity. The final set are assumed to be in a Maunder like state. In the next section magnetic flux transport models used to simulate the evolution of the radial magnetic field at the level of the photosphere on the Sun and other stars are discussed.

2.2 Magnetic flux transport simulations

On large spatial scales, once new magnetic flux has emerged on the Sun, it evolves through the advection processes of differential rotation (Snodgrass, 1983) and meridional flow (Duvall Jr, 1979; Hathaway and Rightmire, 2010). In addition, small convective cells such as super-granulation lead to a random walk of magnetic elements across the solar surface. On spatial scales much larger than super-granules this random walk may be modeled as a diffusive process (Leighton, 1964). Magnetic flux transport simulations (Sheeley Jr, 2005) apply these effects to model the large-scale, long-time evolution of the radial magnetic field B r (θ, φ, t) across the solar surface. In Section 2.2.1 the basic formulation of these models is described. In Section 2.2.2 extensions to the standard model are discussed and, finally, in Sections 2.2.32.2.6 applications of magnetic flux transport models are considered.

2.2.1 Standard model

The standard equation of magnetic flux transport arises from the radial component of the magnetic induction equation under the assumptions that v r = 0 and /∂r = 0.Footnote 1 These assumptions constrain the radial field component to evolve on a spherical shell of fixed radius, where the time evolution of the radial field component is decoupled from the horizontal field components. Under these assumptions, the evolution of the radial magnetic field, B r , at the solar surface (R = 1) is governed by

$$\frac{{\partial B_r }} {{\partial t}} = \frac{1} {{\sin \theta }}\frac{\partial } {{\partial \theta }}\left( {\sin \theta \left( { - u(\theta )B_r + D\frac{{\partial B_r }} {{\partial \theta }}} \right)} \right) - \Omega (0)\frac{{\partial B_r }} {{\partial \varphi }} + \frac{D} {{\sin ^2 \theta }}\frac{{\partial ^2 B_r }} {{\partial \varphi ^2 }} + S(\theta ,\varphi ,t), $$
((1))

where Ω(θ) and u(θ) represent the surface flows of differential rotation and meridional flow, respectively, which passively advect the field, D is the isotropic diffusion coefficient representing superganular diffusion, and finally S(θ, φ, t) is an additional source term added to represent the emergence of new magnetic flux. Figure 2 illustrates two numerical solutions to the flux transport equation when S = 0. Both are initialised with a single bipole in the northern hemisphere, which is then evolved forward in time for 30 solar rotations. The simulations differ only in the tilt of the initial bipole: in the left column, the bipole satisfies Joy's law, while in the right column both polarities lie at the same latitude. The effect of Joy's law has a significant impact on the strength and distribution of B r across the surface of the Sun and, in particular, in the polar regions.

Figure 2
figure 2

Evolution of the radial component of the magnetic field (B r ) at the solar surface for a single bipole in the northern hemisphere with (a) – (c) initial tilt angle of γ = 20° and (d) – (f) γ = 0°. The surface distributions are shown for (a) and (d) the initial distribution, (b) and (e) after 15 rotations, and (c) and (f) after 30 rotations. White represents positive flux and black negative flux and the thin solid line is the Polarity Inversion Line. The saturation levels for the field are set to 100 G, 10 G, and 5 G after 0, 15, and 30 rotations, respectively.

A wide range of studies have been carried out to determine the best fit profiles for the advection and diffusion processes (see DeVore et al., 1985b,a; Wang et al., 1989b; van Ballegooijen et al., 1998). The parameter study of Baumann et al. (2004) demonstrates the effect of varying many of the model parameters. The most commonly accepted values are the following:

  1. 1.

    Differential Rotation: The form that best agrees with the evolution of magnetic flux seen on the Sun (DeVore et al., 1985b) is that of Snodgrass (1983). The profile (Figure 3a) is given by

    $$\Omega (\theta ) = 13.38 - 2.30\cos ^2 \theta - 1.62\cos ^4 \theta \deg /day,$$
    ((2))

    and was determined by cross-correlation of magnetic features seen on daily Mt. Wilson magnetogram observations. The key effect of differential rotation is to shear magnetic fields in an east-west direction, where the strongest shear occurs at mid-latitudes (see Figure 3a). This produces bands of alternating positive and negative polarity as you move poleward (visible after 15 rotations in Figure 2). These bands result in steep meridional gradients which accelerate the decay of the non-axisymmetricFootnote 2 field during periods of low magnetic activity (DeVore, 1987). On the Sun the timescale for differential rotation to act is τdr = 2π/(Ω(0) – Ω(90)) ∼ 1/4 year. Within magnetic flux transport simulations the profile (2) is either applied directly, simulating the actual rotation of the Sun, or with the Carrington rotation rate (13.2 deg/day) subtracted (van Ballegooijen et al., 1998; McCloughan and Durrant, 2002).

  2. 2.

    Supergranular Diffusion: This represents the effect on the Sun's large-scale magnetic field of the small-scale convective motions of supergranules. The term was first introduced by Leighton (1964) to represent the effect of the non-stationary pattern of supergranular cells dispersing magnetic flux across the solar surface. In addition it describes the cancellation of magnetic flux (see Figure 2) when positive and negative magnetic elements encounter oneanother. Initial estimates based on obtaining the correct reversal time of the polar fields (without including meridional flow) were of a diffusion coefficient of D ∼ 770 – 1540 km2 s-1. However, once meridional flow was included to aid the transport of magnetic flux poleward, the value was lowered to around 200 – 600 km2 s-1 which better agrees with observational estimates and is commonly used in models today (DeVore et al., 1985b; Wang et al., 1989b). Globally this gives a time-scale τmf = R 2 /D ∼ 34 – 80 yr, however when considering non-global length scales such as that of individual active regions the timescale is much shorter. Some models have introduced a discrete random walk process as an alternative to the diffusion term which will be discussed in Section 2.2.2 (e.g., Wang and Sheeley Jr, 1994; Schrijver, 2001).

  3. 3.

    Meridional Flow: This effect was the last to be added to what is now known as the standard magnetic flux transport model. It represents an observed weak flow that pushes magnetic flux from the equator to the poles in each hemisphere. The exact rate and profile applied varies from author to author, but peak values of 10 – 20 m s-1 are commonly used. Figure 3b shows a number of meridional flow profiles that have been used by different authors. The typical time-scale for meridional flow is τmf = R/u(θ) ∼ 1 – 2 yr. As current measurements of meridional flow are at the limits of detection, for both the flow rate and profile, this has lead to some variations in the exact profile and values used which will be discussed in more detail in Section 2.2.4. The first systematic study of the effects of meridional flow on the photospheric field was carried out by DeVore et al. (1984). This showed that to obtain a realistic distribution and strength for the polar field, the meridional flow profile must peak at low to mid-latitudes and rapidly decrease to zero at high latitudes. The inclusion of meridional flow was a key development which meant that much lower rates of the diffusion coefficient could be applied, while still allowing the polar fields to reverse at the correct time. It also aids in producing the observed “topknot” latitudinal profile of the polar field (more concentrated than a dipole, Sheeley Jr et al., 1989), and in reproducing the strong poleward surges observed in the butterfly diagram (Wang et al., 1989a). In recent years more significant changes to the profile of meridional flow have been suggested. This will be discussed in Section 2.2.4 and Section 2.2.5.

  4. 4.

    Magnetic Flux Emergence: The final term in Equation (1) is a time-dependent source term which represents a contribution to the radial magnetic field from the emergence of new magnetic bipoles. Most flux transport simulations carry out emergence in a semi-empirical way where the emergence is carried out instantaneously, so that growth of the new bipole is not considered. Instead, only its decay under the action of the advection and diffusion process described above are followed. Inclusion of this term is critical in simulations extending over more than one rotation, to ensure that accuracy of the photospheric field is maintained. Within the literature the source term has been specified in a number of ways and reproduces the main properties of the butterfly diagram, along with varying levels of magnetic activity through single cycles and from one cycle to the next. Different ways in which the source term has been specified include:

    1. (a)

      Observationally determining the properties of new bipoles from daily or synoptic magnetograms so that actual magnetic field configurations found on the Sun may be reproduced (Sheeley Jr et al., 1985; Yeates et al., 2007). Statistical variations of these properties have been applied to model multiple solar cycles of varying activity.

    2. (b)

      Producing synthetic data sets from power law distributions (Harvey and Zwaan, 1993; Schrijver and Harvey, 1994) where the power laws specify the number of bipoles emerging at a given time with a specific area or flux (van Ballegooijen et al., 1998; Schrijver, 2001).

    3. (c)

      Using observations of sunspot group numbers to specify the number of bipoles emerging (Schrijver et al., 2002; Baumann et al., 2006; Jiang et al., 2010a) where the flux within the bipoles can be specified through empirical sunspot area-flux relationships (Baumann et al., 2006; Jiang et al., 2010a). Recently Jiang et al. (2011a) extended this technique back to 1700 using both Group (Hoyt and Schatten, 1998) and Wolf (Wolf, 1861) sunspot numbers.

    4. (d)

      Assimilating observed magnetograms directly into the flux transport simulations (Worden and Harvey, 2000; Schrijver et al., 2003; Durrant et al., 2004). Worden and Harvey (2000) use the flux transport model to produce evolving magnetic synoptic maps from NSO Kitt Peak data. They develop a technique where full-disk magnetograms are assimilated when available and a flux transport model is used to fill in for unobserved or poorly observed regions (such as the far side of the Sun or the poles). Schrijver et al. (2003) use a similar technique but with full-disk MDI observations, demonstrating it over a 5 yr period. For the near side of the Sun, the MDI observations are inserted every 6 h to reproduce the actual field over a 60° degree window where the measurements are most accurate. The authors show that the magnetic flux transport process correctly predicts the return of magnetic elements from the far side, except for the case of flux that emerged on the far-side. To account for this, the authors also include far side acoustic observations for the emergence of new regions (Lindsey and Braun, 2000; Braun and Lindsey, 2001). In contrast to the method of Schrijver et al. (2003), which inputs observations from daily MDI magnetograms, Durrant et al. (2004) inserted fields from synoptic magnetograms once per solar rotation for all latitudes between ± 60°. They used this to investigate the transport of flux poleward and the reversal of the polar field. Recently, the Worden and Harvey (2000) model has been incorporated in a more rigorous data assimilation framework to form the Air Force Data Assimilative Photospheric Flux Transport Model (ADAPT, Arge et al., 2010; Henney et al., 2012). An ensemble of model realisations with different parameter values allow both data and model uncertainties to be incorporated in predictions of photospheric evolution.

A common treatment of the source term is to include only large magnetic bipoles of flux exceeding 1020 Mx. Extensions of the model to include small-scale magnetic regions are described in Section 2.2.2. While many of the parameters for newly emerging bipoles may be determined observationally, or specified through empirical relationships, the parameter about which there is most disagreement in the literature is the variation of the tilt angle (γ) with latitude, or Joy's law. While it should in principle be observed directly for each individual magnetic bipole, this is possible only in recent cycles for which magnetogram observations are available. Traditionally the tilt angle was chosen to vary with latitude λ as γ ∼ λ/2, but more recent studies suggest that a much smaller variation with latitude is required (γ ∼ 0.15λ., Schüssler and Baumann, 2006). The tilt angle is a critical quantity as variations can have a significant effect on the net amount of magnetic flux pushed poleward (see Figure 2) and subsequently on the reversal times of the polar fields and amount of open flux. This will be discussed further in Sections 2.2.4 and 4.1.

Figure 3
figure 3

(a) Profile of differential rotation Ω versus latitude (Snodgrass, 1983) that is most commonly used in magnetic flux transport simulations. (b) Profiles of meridional flow u versus latitude used in various studies. The profiles are from Schüssler and Baumann (2006) (solid line), van Ballegooijen et al. (1998) (dashed line), Schrijver (2001) (dotted line), and Wang et al. (2002b) (dash-dot and dash-dot-dot-dot lines). This figure is based on Figure 3 of Schüssler and Baumann (2006).

2.2.2 Extensions

Since the early flux transport models were produced (DeVore et al., 1985b), new variations have been developed with new formulations and features added by a variety of authors. These include:

  1. 1.

    Extending the basic flux transport model for the large-scale field to include the emergence of small-scale fields down to the size of ephemeral regions (Worden and Harvey, 2000; Schrijver, 2001). Such small scale fields are a necessary element to model the magnetic network and simulate chromospheric radiative losses. Including the small scale emergences leads to a much more realistic description of the photospheric field where discrete magnetic elements can be seen at all latitudes (compare Figure 4 to Figure 2), however since they are randomly oriented small-scale regions do not have a significant effect on large-scale diffusion or on the polar field (Wang and Sheeley Jr, 1991; Worden and Harvey, 2000). To resolve such small-scales, Schrijver (2001) introduced a particle-tracking concept along with rules for the interaction of magnetic elements with one-another. Another important extension introduced by Schrijver (2001) was magneto-convective coupling. With its introduction the early decay of active regions could be considered, where strong field regions diffuse more slowly due to the suppression of convection.

  2. 2.

    Introducing an additional decay term into the basic flux transport model (Schrijver et al., 2002; Baumann et al., 2006). In the paper of Baumann et al. (2006) this term is specified as a linear decay of the mode amplitudes of the radial magnetic field. The term approximates how a radial structure to the magnetic field, along with volume diffusion, would affect the radial component that is constrained to lie on a spherical shell. The authors only include a decay term based on the lowest order radial mode, as it is the only radial mode with a sufficiently long decay time to affect the global field. Through simulating the surface and polar fields from Cycles 13 – 23 (1874 – 2005) and requiring the simulated polar fields to reverse at the correct time as given by both magnetic and proxy observations, the authors deduce a volume diffusivity of the order of 50 – 100 km2 s-1.

  3. 3.

    Coupling the magnetic flux transport model to a coronal evolution model so that both the photospheric and coronal magnetic fields evolve together in time. The first study to consider this was carried out by van Ballegooijen et al. (1998). Limitations of this initial model were that in the coronal volume the radial (B r ) and horizontal (B θ ,B φ ) field components evolve independently from one another, and that no force balance was considered for the coronal field. Later developments (van Ballegooijen et al., 1998, 2000; Mackay and van Ballegooijen, 2006, see Section 3.2.3) removed all of these restrictive assumptions so that a fully coupled evolution of the coronal field along with radial diffusion, radial velocities, and force balance occurred. In an additional study the technique was also extended to include a simplified treatment of the convective zone (van Ballegooijen and Mackay, 2007).

  4. 4.

    Reformulating the flux transport equation into a “synoptic” transport equation, that evolves synoptic magnetic field maps such as those produced by NSO/Kitt PeakFootnote 3 from one Carrington rotation to the next (McCloughan and Durrant, 2002). Their synoptic transport equation takes the form

    $$\frac{{\partial \mathcal{B}_r }} {{\partial \tau }} = - \frac{{f(\theta )}} {{\sin \theta }}\frac{\partial } {{\partial \theta }}(\sin \theta u(\theta )\mathcal{B}_r ) - \frac{{v_\varphi (\theta )f(\theta )}} {{\sin \theta }}\frac{{\partial \mathcal{B}_r }} {{\partial \varphi }} + \frac{{Df(\theta )}} {{\sin \theta }}\frac{\partial } {{\partial \theta }}\left( {\sin \theta \frac{{\partial \mathcal{B}_r }} {{\partial \theta }}} \right) + \frac{{Df(\theta )}} {{\sin ^2 \theta }}\frac{{\partial ^2 \mathcal{B}_r }} {{\partial \varphi ^2 }}, $$
    ((3))

    where B r represents the synoptic magnetogram, τ a time index such that integer values of 2π/Ω correspond to successive synoptic magnetograms, and f(θ) = 1/(1 + v φ (θ)/(Ω sin θ)). The individual terms in Equation (3) are similar to those of Equation (1) but modified to take into account that, due to differential rotation, different latitudes return to central meridian at different times.

  5. 5.

    Although the surface flux transport model is primarily considered in 2D, Leighton (1964) also introduced a reduced 1-D form where the large-scale radial field is spatially averaged in the azimuthal direction. Due to the form of Equation (2), the differential rotation Ω(θ) plays no role in such a model. Recent applications of this 1-D model have been to reconcile surface flux transport with results from axisymmetric, kinematic dynamo simulations in the r-θ plane. Cameron and Schüssler (2007), through describing the emergence of magnetic flux based on sunspot records in a similar manner to that used by Dikpati et al. (2006) and Dikpati and Gilman (2006), determine whether the 1-D surface flux transport model has any predictive properties. They show that in some cases a significant positive correlation can be found between the amount of flux canceling across the equator and the strength of the next cycle. However, when detailed observations of bipole emergences are added, they find that any predictive power may be lost. Cameron et al. (2012) have used the same model to derive constraints on parameter values in the dynamo models, based on observed flux evolution at the surface.

Figure 4
figure 4

Example of the magnetic flux transport simulations of Schrijver (2001) which apply a particle-tracking concept to simulate global magnetic fields down to the scale of ephemeral regions. Image reproduced by permission from Schrijver (2001), copyright by AAS.

The magnetic flux transport models described above have been widely used and extensively compared with observations of the Sun’s large-scale magnetic field. These comparisons are described next for time-scales of up to a single 11-year cycle (Section 2.2.3) and also for multiple solar cycles (Section 2.2.4). In addition, simulations where the profile of meridional flow is changed from that shown in Figure 3b are considered (Section 2.2.5) along with the application of these models to stellar magnetic fields (Section 2.2.6).

2.2.3 Short term applications of magnetic flux transport models

It is generally found that for simulations extending up to a full 11-year solar cycle (termed here short term simulations) magnetic flux transport models are highly successful in reproducing the key features of the evolution of the Sun’s radial magnetic field. These include:

  1. 1.

    The 27 day and 28 – 29 day recurrent patterns seen in plots of the Sun’s mean line-of-sight field. These patterns can be respectively attributed to differential rotation and the combined effect of differential rotation with magnetic flux emergence (Sheeley Jr et al., 1985).

  2. 2.

    Reproducing the rigid rotation of large-scale magnetic features at the equatorial rate, caused by a balance between differential rotation and a poleward drift due to the combined effects of meridional flow and surface diffusion (DeVore, 1987; Sheeley Jr et al., 1987; Wang et al., 1989b).

  3. 3.

    Reproducing the reversal time of the polar fields along with the strength and top-knot profile (DeVore and Sheeley Jr, 1987; Durrant and Wilson, 2003).

  4. 4.

    Matching qualitatively the strength and distribution of the radial magnetic field at the photosphere. Examples include the formation of switchbacks of polarity inversions lines at midlatitudes due to the dispersal of magnetic flux from low latitude active regions (Wang et al., 1989b; van Ballegooijen et al., 1998) and the return of magnetic elements from the far-side of the Sun (Schrijver and DeRosa, 2003).

A full discussion of these results along with their historical development can be found in the review by Sheeley Jr (2005).

Some short term studies have used magnetic flux transport models to predict the nature of the polar magnetic field. These studies assume that the polar fields are solely due to the passive transport of magnetic flux from low to high latitudes on the Sun. An example is Durrant et al. (2001, 2002) where the authors consider a detailed study of high latitude magnetic plumes. These plumes are produced from poleward surges of magnetic flux that originate from activity complexes (Gaizauskas et al., 1983). While a good agreement was found between observations and models, the authors did find some small discrepancies. They attributed these to small bipole emergences outside the normal range of active latitudes considered in flux transport models. This indicates that the field at high latitudes is not solely the result of magnetic fields transported from low latitudes and a local dynamo action may play a role. To date the only surface flux transport models to include such emergences are those of Schrijver (2001) and Worden and Harvey (2000).

In an additional study, Durrant (2002) and Durrant and Wilson (2003) investigated in detail the reversal times of the polar magnetic fields, comparing the reversal times deduced from KP normal component magnetograms with those found in synoptic flux transport simulations. To test whether magnetic observations or simulations gave the best estimate they deduced the locations of PILs at high latitudes from Hα filament data. They then compared these to (i) the observed PIL locations deduced from Kitt Peak data and (ii) the simulated PILs from the synoptic flux transport equations. The study found that above 70° latitude the location of the PIL — as given by the Hα filaments — was in fact better determined by the flux transport process rather than the direct magnetic field observations. This was due to the high level of uncertainty and error in measuring the field at high latitudes as a result of foreshortening and solar B angle effects (Worden and Harvey, 2000). However, the simulation had to be run for over 20 rotations so that the highlatitude fields in the polar regions were solely the product of flux transport processes. This ensured that any systematic errors in the observations at high latitudes were removed.

A recent application of magnetic flux transport models has been to improve forecasts of the solar 10.7 cm (2.8 GHz) radio flux, using an empirical relation between this quantity and total “sunspot” and “plage” fields in the simulation (Henney et al., 2012). The 10.7 cm radio flux is widely used by the space weather community as a proxy for solar activity, with measurements dating back to 1947. For other aspects of space weather, flux transport models need to be coupled to coronal models: these applications will be discussed in Section 4.

2.2.4 Multiple solar cycle applications of magnetic flux transport models

Agreement between observations and magnetic flux transport simulations is generally good over time-scales of less than 11 years. But when magnetic flux transport simulations are run for multiple 11-year solar cycles, some inconsistencies are found. One of the first long term simulations was carried out by Schrijver et al. (2002) who simulated the photospheric field over 32 cycles from the Maunder minimum (1649) to 2002. With detailed observations of magnetic bipoles available only for the last few solar cycles, the authors used synthetic data, scaling the activity level to the observed sunspot number. They also assumed that the flux transport parameters remained the same from one cycle to the next, and used a meridional flow profile that peaked at mid-latitudes with a value of 14 m s-1. Since polar field production varies with the amount of emerging flux, the authors found that if the high latitude field is solely described by the passive advection of magnetic flux from the active latitudes, then reversals of the polar field within each cycle may be lost. This is especially true if a series of weak cycles follows stronger ones. Although reversals were lost, they did return in later cycles, but the reversal in Cycle 23 did not match the observed time (see Figure 1 of Schrijver et al., 2002). To ensure that reversals occurred in every cycle, Schrijver et al. (2002) introduced a new exponential decay term for the radial field in Equation (1). This acts to reduce the strength of the large scale field over long periods of time and prevents excess polar fields from building up. They found that a decay time of 5 – 10 yr is optimal in allowing polar field reversals to occur from one cycle to the next. Later, Baumann et al. (2006) also found that a decay term was required in long term simulations to maintain the reversal of the polar field. They formulated the decay term in terms of 3D diffusion of the magnetic field constrained on the 2D solar surface.

While Schrijver et al. (2002) and Baumann et al. (2006) introduced a new physical term into the 2D flux transport model, Wang et al. (2002a) considered a different approach. Rather than keeping the advection due to meridional flow constant from one cycle to the next, they varied the strength of meridional flow such that stronger cycles were given a faster flow. If a factor of two change occurs between strong and weak cycles, then reversals of the polar field may be maintained (see Figure 2 of Wang et al., 2002a). This is because faster meridional flow leads to less flux canceling across the equator, and less net flux transported poleward. If correct, this suggests that in stronger cycles the magnetic fluxes in each hemisphere are more effectively isolated from one another. To complicate matters further, Cameron et al. (2010) put forward a third possibility for maintaining polar field reversals. This was based around an observed cycle-to-cycle variation in the tilt angle of sunspot groups (Dasi-Espuig et al., 2010). If there is an anti-correlation between cycle strength and tilt angles, where stronger cycles are assumed to have smaller tilt angles, then neither variable rates of meridional flow, or extra decay terms are required to maintain polar field reversals from one cycle to the next. To show this, Cameron et al. (2010) consider a simulation extending from 1913 – 1986 and show that with such a tilt angle variation, the polar field reversal may be maintained. Although a decreased tilt angle alone was sufficient to maintain polar field reversals in this study, a longer term study by Jiang et al. (2011b) found that the radial decay term had to be re-introduced when carrying out simulations from 1700 to present (Jiang et al., 2011a). The authors attributed this to inaccuracies in the input data of observed activity levels. However, it is not clear that observed tilt angles in the recent Cycle 23 were sufficiently reduced for this alone to explain the low polar fields in 2008 (Schrijver and Liu, 2008).

This discussion indicates that multiple combinations of model parameters may produce qualitatively the same result. This is illustrated by the study of Schrijver and Liu (2008) who consider the origin of the decreased axial dipole moment found in 2008 compared to 1997. They simulated the global field from 1997 – 2008 and deduced that their previously included decay term was insufficient by itself to account for the lower dipole moment. Instead, they found that to reproduce the observations a steeper meridional flow gradient is required at the equator. This steeper gradient effectively isolates the hemispheres, as in Wang et al. (2002a). However, in contrast to Wang et al. (2002a), in their study Schrijver and Liu (2008) did not have to increase the meridional flow rate, which was maintained at around 10 m s-1. In a similar study, Jiang et al. (2011c) illustrated two additional ways of reproducing the lower polar field strengths in 2008, namely (i) a 28% decrease in bipole tilt angles, or (ii) an increase in the meridional flow rate from 11 m s-1 to 17 m s-1. They, however, found that the first case delayed the reversal time by 1.5 yr, so was inconsistent with observations. Unfortunately, all of these possible solutions are within current observational limits, so none can be ruled out.

2.2.5 Variations in the meridional flow profile

While Wang et al. (2002a) and Schrijver and Liu (2008) introduce overall variations of the meridional flow profile and rate, they maintain a basic poleward flow profile in each hemisphere. Recent studies have considered more significant (and controversial) changes to the meridional flow, motivated by helioseismic observations. Two broad types of variation have been considered: a counter-cell near the poles, and an inflow towards activity regions.

Jiang et al. (2009) considered the possibility of the existence of a counter-cell of reversed meridional flow near the poles. This change was motivated directly by observations showing that in Cycles 21 and 22 the polar field was strongly peaked at the poles, while in Cycle 23 it peaked at a latitude of 75° and then reduced by nearly 50% closer to the pole (Raouafi et al., 2008). While the observed structure in Cycles 21 and 22 was consistent with a meridional flow profile that either extended all the way to the poles, or switched off at around 75°, the Cycle 23 structure was not. To investigate what form of meridional flow could reproduce this, Jiang et al. (2009) introduced a counter-cell of meridional flow beyond 70° latitude and showed that the observed distribution of magnetic flux in Cycle 23 could be achieved with a counter-cell rate of 2 – 5 m s-1. The authors showed that the counter-cell need not be constant and that similar results could be obtained if it existed only in the declining phase. To date existence of such a counter-cell has not been verified beyond doubt even though some results from helioseismology support its existence (Haber et al., 2002; González Hernández et al., 2006). Also, due to the difference in profile of polar magnetic fields in Cycles 21 and 22 it is unlikely that such a counter-cell existed in these cycles.

The second type of flow perturbation suggested by observations are activity-dependent inflows towards the central latitude of emergence of the butterfly diagram (Gizon and Rempel, 2008; González Hernández et al., 2008). At least in part, this modulation of the axisymmetric meridional flow component results from the cumulative effect of the sub-surface horizontal flows converging towards active regions (Hindman et al., 2004). In an initial study, DeRosa and Schrijver (2006) simulated the effect of (non-axisymmetric) inflows towards active regions by adding an advection term to the flux transport model, its velocity proportional to the horizontal gradient of radial magnetic field. They found that if surface speeds exceed 10 m s-1 then an unrealistic “clumping” of magnetic flux that is not observed occurs. More recently, Jiang et al. (2010b) have considered the effect of an axisymmetric flow perturbation towards active latitudes of speed 3 – 5 m s-1, similar to perturbations observed during Cycle 23. The authors consider how such a flow effects the polar field distribution though simulated solar cycles. The main effect is to decrease the separation of magnetic bipoles, subsequently leading to more cancellation in each hemisphere and less net flux transport poleward. This leads to an 18% decrease in the polar field strength compared to simulations without it. While this is a significant decrease, the authors note that it cannot (alone) account for the weak polar fields observed at the end of Cycle 23 as those fields decreased by more than a factor of 2. However, in conjunction with a variable meridional flow from one cycle to the next (Wang et al., 2002a) or the steepening gradient of the flow (Schrijver and Liu, 2008), the inflow may have a significant effect. Cameron and Schüssler (2010) show that one can incorporate such axisymmetric flow perturbations by setting the speed at a given time proportional to the latitudinal gradient of longitude-averaged B r . Not only does this generate appropriate inflows, but these can explain the observed solar cycle variation of the P 12 Legendre component of meridional flow (Hathaway and Rightmire, 2010) without changing the background flow speed.

2.2.6 Stellar applications of magnetic flux transport models

Due to advances in measuring stellar magnetic fields, magnetic flux transport models have recently been applied in the stellar context. One of their main applications has been to consider how polar or high latitude intense field regions or stellar spots may arise in rapidly rotating stellar systems (Strassmeier and Rice, 2001). Initial studies (Schrijver and Title, 2001; Mackay et al., 2004) used as a starting point parameters from solar magnetic flux transport simulations. These parameters were then varied to reproduce the key observational properties of the radial magnetic fields on stars as deduced through ZDI measurements. In the paper of Schrijver and Title (2001) the authors consider a very active cool star of period 6 days. A key feature of their simulations is that they fix all parameters to values determined for the Sun and only vary the emergence rate to be 30 times solar values. They show that even if solar emergence latitudes are maintained, then the flux transport effects of meridional flow and surface diffusion are sufficient to transport enough flux to the poles to produce a polar spot. In Figure 5, a comparison of typical solar (Figure 5a) and stellar (Figure 5b) magnetic field configurations from Figure 2 of Schrijver and Title (2001) can be seen. For the case of the Sun the magnetic fields saturate at ± 70 Mx cm-2 and for the star ± 700 Mx cm-2. Both images show the polar region from a latitude of 40°. A clear difference can be seen, where for the Sun the pole has a weak unipolar field. For the rapidly rotating star there is a unipolar spot with a ring of strong opposite polarity flux around it. The existence of this ring is partly due to the non-linear surface diffusion used in Schrijver and Title (2001) which causes intense field regions to diffuse more slowly (see Section 2.2.2).

Figure 5
figure 5

Comparison of (a) solar and (b) stellar magnetic field configurations from Schrijver and Title (2001). For the stellar system all flux transport parameters are held fixed to solar values and only the emergence rate is increased to 30 times solar values. For the case of the Sun the magnetic fields saturate at ± 70 Mx cm-2 and for the star ± 700 Mx cm-2. Image reproduced by permission from Schrijver and Title (2001), copyright by AAS.

In contrast to the unipolar poles modeled by Schrijver and Title (2001), some ZDI observations of rapidly rotating stars show intermingling of opposite polarities at high latitudes, where intense fields of both polarities lie at the same latitude and are not nested. An example of this can be seen in AB Dor (Figure 1b) which has a rotation period of 1/2 day. Mackay et al. (2004) showed that in order to produce strong intermingled polarities within the polar regions more significant changes are required relative to the solar flux transport model. These include increasing the emergence latitude of new bipoles from 40° to 70° and increasing the rate of meridional flow from 11 m s-1 to 100 m s-1. Increased emergence latitudes are consistent with an enhanced Coriolis force deflecting more flux poleward as it travels through the convection zone (Schüssler and Solanki, 1992). However, at the present time, the predicted enhanced meridional flow is still below the level of detection. While the simulations of Schrijver and Title (2001) and Mackay et al. (2004) have successfully reproduced key features in the magnetic field distributions of rapidly rotating stars, as yet they do not directly simulate the observations as has been done for the Sun. Presently, there is insufficient input data on the emergence latitudes of new bipoles in order to carry out such simulations. In the study of Ičsik et al. (2007) the authors used magnetic flux transport simulations to estimate the lifetime of starspots as they are transported across the surface of a star. The authors show that many factors may effect the star spot lifetime, such as the latitude of emergence and differential rotation rate. In particular the authors show that for rapidly rotating stars the lifetime of spots may be 1 – 2 months, however the lifetime is lower for stars with strong differential rotation (AB Dor) compared to those with weak differential rotation (HR 1099).

In recent years magnetic flux transport simulations have been combined with models for the generation and transport of magnetic flux in the stellar interior (Holzwarth et al., 2006; Ičsik et al., 2011) to predict interior properties on other stars. The first study to link the pre- and posteruptive transport properties was carried out by Holzwarth et al. (2006). In this study, the authors quantified what effect the enhanced meridional flow predicted by Mackay et al. (2004) would have on the dynamics and displacement of magnetic flux tubes as they rise though the convective zone (Moreno-Insertis, 1986; Schüussler et al., 1994). The authors found that the enhanced meridional flow leads to a non-linear displacement of the flux tubes, where bipoles were displaced to emerge at higher latitudes. This was consistent with the required higher latitude of emergence to produce the desired intermingling in the poles. However, in a more detailed model, Ičsik et al. (2011) combine a thin-layer α-Ω dynamo model with a convective zone transport model and surface flux transport model to provide a complete description of the evolution of magnetic field on stars. Through this combined model the authors show that for rapid rotators (period ∼ 2 days), due to the effect of the Coriolis force, the surface signature of the emergence and subsequent transport may not be a clear signal of the dynamo action that created it. This has important consequences for both the observation and interpretation of magnetic fields on other stars. Living Reviews

3 Coronal Magnetic Field Models

The magnetic fields observed in the solar photosphere extend out into the corona, where they structure the plasma, store free magnetic energy and produce a wide variety of phenomena. While the distribution and strength of magnetic fields are routinely measured in the photosphere, the same is not true for the corona, where the low densities mean that such measurements are very rare (Cargill, 2009). To understand the nature of coronal magnetic fields, theoretical models that use the photospheric observations as a lower boundary condition are required.

In this section, we survey the variety of techniques that have been developed to model the coronal magnetic field based on the input of photospheric magnetic fields. We limit our discussion to global models in spherical geometry, and do not consider, for example, models for the magnetic structure of single active regions. In future revisions of this review, we will extend our discussion to include observations of both coronal and prominence magnetic fields in order to determine the validity of these models. Applications of the global models are the subject of Section 4.

A difficulty arising with any such attempt at a global model based on observations of the photospheric magnetic field, is that we can only observe one side of the Sun at a single time (except for some limited observations by STEREO), but need data for all longitudes. This problem may be addressed either by compiling time series of full-disk observations into a synoptic observed magnetogram, or by assimilating individual active regions into a time-dependent surface flux transport model (Section 2.2). In this section, we will assume that the photospheric distribution of B r (RR, θ, φ) at a given time has already been derived by one of these techniques, and consider how to construct the magnetic field of the corona.

Virtually all global models to date produce equilibria. This is reasonable since the coronal magnetic evolution on a global scale is essentially quasi-static. The magnetic field evolves through a continuous sequence of equilibria because the Alfvén speed in the corona (the speed at which magnetic disturbances propagate) is of the order of 1000 km s-1, much faster than the large-scale surface motions that drive the evolution (1 – 3 km s-1). To model this quasi-static evolution, nearly all models produce a series of independent “single-time” extrapolations from the photospheric magnetic fields at discrete time intervals, without regard to the previous magnetic connectivity that existed in the corona. The exception is the flux transport magneto-frictional model (Section 3.2) which produces a continuous sequence of equilibria. In this model, the equilibrium assumption is violated locally for brief periods during dynamical events such as CMEs or flares. Individual CME events are becoming practical to simulate in fully time-dependent MHD codes, but equilibrium models are still the only feasible way to understand when and where they are likely to occur.

We divide the techniques into four categories. The first three are potential field source surface models (Section 3.1), force-free field models (Section 3.2), and magnetohydrostatic models (Section 3.3). All of these models output only (or primarily) the magnetic field. The fourth category are full MHD models (Section 3.4), which aim to self-consistently describe both the magnetic field and other plasma properties.

3.1 Potential field source surface models

The most straightforward and therefore most commonly used technique for modeling the global coronal magnetic field is the so-called Potential Field Source Surface (PFSS) model (Schatten et al., 1969; Altschuler and Newkirk Jr, 1969). We note that an implementation of the PFSS model by Luhmann et al. (2002), using photospheric magnetogram data from Wilcox Solar Observatory, is available to run “on demand” at http://ccmc.gsfc.nasa.gov/. The key assumption of this model is that there is zero electric current in the corona. The magnetic field is computed on a domain RrRss, between the photosphere, R, where the radial magnetic field distributionFootnote 4 is specified and an outer “source surface”, Rss, where the boundary conditions B θ = B φ = 0 are applied. The idea of the source surface is to model the effect of the solar wind outflow, which distorts the magnetic field away from a current-free configuration above approximately 2R (once the magnetic field strength has fallen off sufficiently). Distortion away from a current-free field implies the presence of electric currents at larger radii: a self-consistent treatment of the solar wind requires full MHD description (Section 4.4), but an artificial “source surface” boundary approximates the effect on the magnetic field in the lower corona. The PFSS model has been used to study a wide variety of phenomena ranging from open flux (see Section 4.1), coronal holes (Section 4.2), and coronal null points (Cook et al., 2009) to magnetic fields in the coronae of other stars (Jardine et al., 2002a).

The requirement of vanishing electric current density in the PFSS model means that ∇ × B = 0, so that

$$B = - \nabla \Psi ,$$
((4))

where Ψ is a scalar potential. The solenoidal constraint ∇ · B = 0 then implies that Ψ satisfies Laplace's equation

$$\nabla ^2 \Psi = 0,$$
((5))

with Neumann boundary conditions

$$\left. {\frac{{\partial \Psi }} {{\partial r}}} \right|_{r = R_ \odot } = - B_r (R_ \odot ,\theta ,\varphi ),\left. {\frac{{\partial \Psi }} {{\partial \theta }}} \right|_{r = R_{ss} } = \left. {\frac{{\partial \Psi }} {{\partial \varphi }}} \right|_{r = R_{ss} } = 0. $$
((6))

Thus, the problem reduces to solving a single scalar PDE. It may readily be shown that, for fixed boundary conditions, the potential field is unique and, moreover, that it is the magnetic field with the lowest energy (ƒ V B2/(2μ0)dV) for these boundary conditions.

Solutions of Laplace’s equation in spherical coordinates are well known (Jackson, 1962). Separation of variables leads to

$$\Psi (r,\theta ,\varphi ) = \sum\limits_{l = 0}^\infty {\sum\limits_{m = - l}^l {\left[ {f_{lm} r^l + g_{lm} r^{ - (l + 1)} } \right]} } P_l^m (\cos \theta )e^{im\varphi } ,$$
((7))

where the boundary condition at Rss implies that g lm = f lm R 2l+1ss for each pair l, m, and P m l are the associated Legendre polynomials. The coefficients f lm are fixed by the photospheric B r (R, θ, φ) distribution to be

$${f_{lm}} = - {\left( {lR_ \odot ^{l - 1} + (l + 1)R_ \odot ^{ - l - 2}R_{ss}^{2l + 1}} \right)^{ - 1}}{b_{lm}},$$
((8))

where b lm are the spherical harmonic coefficients of this B r (R, θ, φ) distribution. The solution for the three field components may then be written as

$${B_r}(r,\theta ,\varphi ) = \sum\limits_{l = 0}^\infty {\sum\limits_{m = - l}^l {{b_{lm}}{c_l}(r)P_l^m(\cos \theta ){e^{im\varphi }}} } ,$$
((9))
$${B_\theta }(r,\theta ,\varphi ) = - \sum\limits_{l = 0}^\infty {\sum\limits_{m = - l}^l {{b_{lm}}{d_l}(r)\frac{{dP_l^m(\cos \theta )}}{{d\theta }}{e^{im\varphi }}} } , $$
((10))
$${B_\varphi }(r,\theta ,\varphi ) = - \sum\limits_{l = 0}^\infty {\sum\limits_{m = - l}^l {\frac{{im}}{{\sin \theta }}{b_{lm}}{d_l}(r)P_l^m(\cos \theta ){e^{im\varphi }}} } , $$
((11))

where

$${c_l}(r) = {\left( {\frac{r}{{{R_ \odot }}}} \right)^{ - l - 2}}\left[ {\frac{{l + 1 + l{{(r/{R_{ss}})}^{2l + 1}}}}{{l + 1 + l{{({R_ \odot }/{R_{ss}})}^{2l + 1}}}}} \right], $$
((12))
$${d_l}(r) = {\left( {\frac{r}{{{R_ \odot }}}} \right)^{ - l - 2}}\left[ {\frac{{l - {{(r/{R_{ss}})}^{2l + 1}}}}{{l + 1 + l{{({R_ \odot }/{R_{ss}})}^{2l + 1}}}}} \right]. $$
((13))

In practice, only a finite number of harmonics l = 1, ..., L max are included in the construction of the field, depending on the resolution of the input photospheric B r distribution. Notice that the higher the mode number l, the faster the mode falls off with radial distance. This implies that the magnetic field at larger r is dominated by the low order harmonics. So calculations of the Sun’s open flux (Section 4.1) with the PFSS model do not require large Lmax. An example PFSS extrapolation is shown in Figure 10a.

The PFSS model has a single free parameter: the radius Rss of the source surface. Various studies have chosen this parameter by optimising agreement with observations of either white-light coronal images, X-ray coronal hole boundaries, or, by extrapolating with in situ observations of the interplanetary magnetic field (IMF). A value of Rss = 2.5R. is commonly used, following Hoeksema et al. (1983), but values as low as Rss = 1.3R have been suggested (Levine et al., 1977). Moreover, it appears that even with a single criterion, the optimal Rss can vary over time as magnetic activity varies (Lee et al., 2011).

A significant limitation of the PFSS model is that the actual coronal magnetic field does not become purely radial within the radius where electric currents may safely be neglected. This is clearly seen in eclipse observations (e.g., Figure 6a), as illustrated by Schatten (1971), and in the early MHD solution of Pneuman and Kopp (1971). Typically, real polar plumes bend more equatorward than those in the PFSS model, while streamers should bend more equatorward at Solar Minimum and more poleward at Solar Maximum.

Figure 6
figure 6

Comparison between (a) a drawing of the real eclipse corona, (b) the standard PFSS model, and (c) the Schatten (1971) current sheet model. The current sheet model better reproduces the shapes of polar plumes and of streamer axes. Image adapted from Figure 1 of Zhao and Hoeksema (1994).

Schatten (1971) showed that the PFSS solution could be improved by replacing the source surface boundary Rss with an intermediate boundary at Rcp = 1.6R, and introducing electric currents in the region r > Rcp. To avoid too strong a Lorentz force, these currents must be limited to regions of weak field, namely to sheets between regions of B r > 0 and B r < 0. These current sheets support a more realistic non-potential magnetic field in r > Rcp. The computational procedure is as follows:

  1. 1.

    Calculate B in the inner region rRcp from the observed photospheric B r , assuming an “exterior” solution (vanishing as r → ∞).

  2. 2.

    Re-orientate B(Rcp, θ, φ) to ensure that B r > 0 everywhere (this will temporarily violate ∇ · B = 0 on the surface r = Rcp).

  3. 3.

    Compute B in the outer region r > Rcp using the exterior potential field solution, but matching all three components of B to those of the re-orientated inner solution on Rcp. (These boundary conditions at Rcp are actually over-determined, so in practice the difference is minimised with a least-squares optimisation).

  4. 4.

    Restore the original orientation. This creates current sheets where B r = 0 in the outer region, but the magnetic stresses will balance across them (because changing the sign of all three components of B leaves the Maxwell stress tensor unchanged).

This model generates better field structures, particularly at larger radii (Figure 6c), leading to its use in solar wind/space weather forecasting models such as the Wang-Sheeley-Arge (WSA) model (Arge and Pizzo, 2000, available at http://ccmc.gsfc.nasa.gov/). The same field-reversal technique is used in the MHS-based Current-Sheet Source Surface model (Section 3.3), and was used also by Yeh and Pneuman (1977), who iterated to find a more realistic force-balance with plasma pressure and a steady solar wind flow.

By finding the surface where B r /|B| = 0.97 in their MHD model (Section 3.4), Riley et al. (2006) locate the true “source surface” at r ≈ 13R, comparable to the Alfvén critical point where the kinetic energy density of the solar wind first exceeds the energy density of the magnetic field (Zhao and Hoeksema, 2010). The Alfvén critical points are known not to be spherically symmetric, and this is also found in the Riley et al. (2006) model. In fact, a modified PFSS model allowing for a non-spherical source surface (though still at 2 – 3R) was proposed by Schulz et al. (1978). For this case, before solving with the source surface boundary, the surface shape is first chosen as a surface of constant |B| for the unbounded potential field solution.

3.2 Nonlinear force-free field models

While potential field solutions are straightforward to compute, the assumption of vanishing electric current density (j = 0) in the volume renders them unable to model magnetic structures requiring non-zero electric currents. Observations reveal such structures both in newly-emerged active regions, e.g., X-ray sigmoids, and outside active latitudes, including long-lived H.. filament channels and coronal magnetic flux ropes. A significant limitation of the potential field is that it has the lowest energy compatible with given boundary conditions. Yet many important coronal phenomena derive their energy from that stored in the magnetic field; this includes large-scale eruptive events (flares and CMEs), but also small-scale dynamics thought to be responsible for heating the corona. We still lack a detailed understanding of how these events are initiated. Such an understanding cannot be gained from potential field models, owing to the lack of free magnetic energy available for release. The models in this section, based on the force-free assumption, allow for electric currents and, hence, free magnetic energy.

In the low corona of a star like the Sun, the magnetic pressure dominates over both the gas pressure and the kinetic energy density of plasma flows. Thus to first approximation an equilibrium magnetic field must have a vanishing Lorentz force,

$$j \times B = 0,$$
((14))

where ∇ × B/μ0 is the current density. Such a magnetic field is called force-free. It follows that

$$\nabla \times B = \alpha B$$
((15))

for some scalar function α(r) which is constant along magnetic field lines (this follows from ∇ · B = 0). The function α(r) depends only on the local geometry of field lines, not on the field strength; it quantifies the magnetic “twist” (Figure 7). The potential field j = 0 is recovered if α = 0 everywhere.

Figure 7
figure 7

Twisted magnetic fields. Panels (a) and (b) show the direction of twist for a force-free field line with (a) positive α, and (b) negative α, with respect to the potential field α = 0. Panels (c) and (d) show X-ray sigmoid structures with each sign of twist, observed in active regions with the Hinode X-ray Telescope (SAO, NASA, JAXA, NAOJ). These images were taken on 2007 February 16 and 2007 February 5, respectively.

Unfortunately, for non-zero α, the extrapolation of coronal force-free magnetic fields is not straightforward. If α is constant everywhere, then taking the curl of (15) leads to the vector Helmholtz equation

$$(\nabla ^2 + \alpha ^2 )B = 0,$$
((16))

which is linear and may be solved analytically in spherical harmonics (Durrant, 1989). While some early attempts were made to model the global corona using such linear force-free fields (Nakagawa, 1973; Levine and Altschuler, 1974), their use is limited for two main reasons. Firstly, the constant α in such a solution scales as L-1, where L is the horizontal size of the area under consideration. So for global solutions, only rather small values of α, close to potential, can be used. Secondly, observations of twist in magnetic structures indicate that both the sign and magnitude of α ought to vary significantly between different regions on the Sun. In addition, a mathematical problem arises if one attempts to apply the same boundary conditions as in the potential field model (namely a given distribution of B r on r = R, and B θ = B φ = 0 on r = Rss). A strictly force-free field satisfying B θ = B φ = 0 and B r ≠ 0 on an outer boundary must have α = 0 on that boundary (Aly and Seehafer, 1993), so that α = 0 on all open magnetic field lines. Unless α = 0, this is incompatible with the linear force-free field.Footnote 5 From this, it is clear that the use of linear force-free fields is restrictive and will not be discussed further. Instead, we will consider nonlinear force-free fields where α is a function of position.

The problem of extrapolating nonlinear force-free fields from given photospheric data is mathematically challenging, with open questions about the existence and uniqueness of solutions. Nevertheless, several numerical techniques have been developed in recent years, though they have largely been applied to single active regions in Cartesian geometry (Schrijver et al., 2006; DeRosa et al., 2009). Applications to global solutions in spherical geometry are in their infancy: we describe here the three main approaches tried.

3.2.1 Optimisation method

Wiegelmann (2007) has developed a method for numerically computing nonlinear force-free fields in spherical geometry, based on the optimisation procedure of Wheatland et al. (2000) and using the vector magnetic field in the photosphere as input. The idea is to minimise the functional

$$L[B] = \int_V {\left( {\frac{{|(\nabla \times B) \times B|^2 }} {{B^2 }} + |\nabla \cdot B|^2 } \right)} dV, $$
((17))

since if L = 0 then j × B = 0 and ∇.B = 0 everywhere in the volume V. Differentiating (17) with respect to t, Wheatland et al. (2000) show that

$$\frac{1} {2}\frac{{dL}} {{dt}} = - \int_V {\frac{{\partial B}} {{\partial t}} \cdot FdV} - \oint_{\partial V} {\frac{{\partial B}} {{\partial t}} \cdot GdS} , $$
((18))

where

$$F = \nabla \times (\Omega \times B) - \Omega \times (\Omega \times B) - \nabla (\Omega \cdot B) + \Omega (\nabla \cdot B) + \Omega ^2 B,$$
((19))
$$G = n \times (\Omega \times B) - n(\Omega \cdot B),$$
((20))
$$\Omega = \frac{1} {{B^2 }}\left[ {(\nabla \times B) \times B - (\nabla \cdot B)B} \right].$$
((21))

Evolving the magnetic field so that

$$\frac{{\partial B}} {{\partial t}} = \mu F$$
((22))

for some μ > 0, and imposing B/∂t = 0 on the boundary ∂V, will ensure that L decreases during the evolution. To apply this procedure, Wiegelmann (2007) takes an initial potential field extrapolation, and replaces B θ and B φ on r = R. with the measured horizontal magnetic field from an observed vector magnetogram. The 3D field is then iterated with Equation (22) until L is suitably small.

Wiegelmann (2007) demonstrates that this method recovers a known analytical force-free field solution (Figure 8), although the most accurate solution is obtained only if the boundary conditions at both the photosphere and upper source surface are matched to the analytical field. The method is yet to be applied to observational data on a global scale, largely due to the limitation that vector magnetogram data are required for the lower boundary input. Such data are not yet reliably measured on a global scale, though SDO or SOLIS should be a key step forward in obtaining this data. Another consequence is that even where they are measured there is the complication that the photospheric magnetic field is not force-free. Pre-processing techniques to mitigate these problems are under development (Tadesse et al., 2011).

Figure 8
figure 8

Test of the optimisation method for nonlinear force-free fields with an analytical solution (figure adapted from Wiegelmann, 2007). Panel (a) shows the analytical solution (Low and Lou, 1990) and (b) shows the PFSS extrapolation used to initialise the computation. Panel (c) shows the resulting NLFFF when the boundary conditions at both r = R and r = Rss are set to match the analytical solution, while panel (d) shows how the agreement is poorer if the analytical boundary conditions are applied only at the photosphere (the more realistic case for solar application).

3.2.2 Force-free electrodynamics method

Contopoulos et al. (2011) have recently proposed an alternative method to compute global nonlinear force-free field extrapolations, adapted from a technique applied to pulsar magnetospheres. It has the advantage that only the radial component B r . is required as observational input at the photosphere. The computation is initialised with an arbitrary 3D magnetic field such as B = 0 or that of a dipole, and also with zero electric field E. The fields E and B are then evolved through the equations of force-free electrodynamics,

$$\frac{{\partial B}} {{\partial t}} = - c\nabla \times E,$$
((23))
$$\frac{{\partial E}} {{\partial t}} = c\nabla \times B - 4\pi j,$$
((24))
$$\nabla \cdot B = 0,$$
((25))
$$\rho _e E + \frac{1} {c}j \times B = 0,$$
((26))

where ρ e = ∇ · E/(4π) is the electric charge density. The force-free condition (26) enables j to be eliminated from (24) since it follows that

$$j = \frac{c} {{4\pi }}\nabla \cdot E\frac{{E \times B}} {{B^2 }} + \frac{c} {{4\pi }}\left( {\frac{{B \cdot \nabla \times B - E \cdot \nabla \times E}} {{B^2 }}} \right)B.$$
((27))

Equations (23) and (24) are integrated numerically with time-dependant E θ , E φ imposed on the lower boundary r = R. These are chosen so that B r gradually evolves toward its required distribution. The photospheric driving injects electrodynamic waves into the corona, establishing a network of coronal electric currents. The outer boundary is chosen to be non-reflecting and perfectly absorbing, mimicking empty space. As the photosphere approaches the required distribution, the charge density and electric fields diminish, but the coronal currents remain. As E → 0 a force-free field is reached.

The authors have tested their method with an observed synoptic magnetogram: an example force-free field produced is shown in Figure 9. They find that active region magnetic fields are reproduced quite rapidly, but convergence to the weaker fields in polar regions is much slower. Slow convergence is undesirable because numerical diffusion was found to erode the coronal currents before the photosphere reached the target configuration. However, the convergence rate could be improved by choosing the initial condition to be a dipole, approximating the average polar field in the magnetogram. This highlights an important feature of the method: non-uniqueness. The force-free field produced is not defined solely by the photospheric boundary condition, but depends both on (i) the choice of initialisation and (ii) the path followed to reach the final state. Contopoulos et al. (2011) suggest that one way to choose between possible solutions would be to incorporate measurements from vector magnetograms. The model described in the next section takes a different approach in treating the construction of the coronal magnetic field as an explicitly time-dependent problem.

Figure 9
figure 9

A force-free magnetic field produced by the force-free electrodynamics method for Carrington rotation 2009. Field line colours show the twist parameter α, coloured red where α > 0, blue where α < 0, and light green where α ≈ 0. Image reproduced by permission from Contopoulos et al. (2011), copyright by Springer.

3.2.3 Flux transport and magneto-frictional method

Recently, van Ballegooijen et al. (2000) and Mackay and van Ballegooijen (2006) have developed a new technique to study the long-term evolution of coronal magnetic fields, which has now been applied to model the global solar corona (Yeates et al., 2007, 2008b). The technique follows the build-up of free magnetic energy and electric currents in the corona by coupling together two distinct models. The first is a data driven surface flux transport model (Yeates et al., 2007). This uses observations of newly emerging magnetic bipoles to produce a continuous evolution of the observed photospheric magnetic flux over long periods of time. Coupled to this is a quasi-static coronal evolution model (Mackay and van Ballegooijen, 2006; Yeates et al., 2008b) which evolves the coronal magnetic field through a sequence of nonlinear force-free fields in response to the observed photospheric evolution and flux emergence. The model follows the long-term continuous build-up of free magnetic energy and electric currents in the corona. It differs significantly from the extrapolation approaches which retain no memory of magnetic flux or connectivity from one extrapolation to the next.

The photospheric component of the model evolves the radial magnetic field B r on r = R with a standard flux transport model (Section 2.2.1), except that newly emerging bipolar active regions are inserted not just in the photosphere but also in the 3D corona. These regions take an analytical form (Yeates et al., 2007), with parameters (location, size, flux, tilt angle) chosen to match observed active regions. An additional twist parameter allows the emerging 3D regions to be given a non-zero helicity: in principle this could be determined from vector magnetogram observations, but these are not yet routinely available.

The coronal part of the model evolves the large-scale mean field (van Ballegooijen et al., 2000) according to the induction equation

$$\frac{{\partial A}} {{\partial t}} = v \times B + \mathcal{E},$$
((28))

, where B = ∇ × A is the mean magnetic field (with the vector potential A in an appropriate gauge). The mean electromotive force ε describes the net effect of unresolved small-scale fluctuations, for example, braiding and current sheets produced by interaction with convective flows in the photosphere. Mackay and van Ballegooijen (2006) and Yeates et al. (2008b) assume the form

$$\mathcal{E} = - \eta j,$$
((29))

where

$$\eta = \eta _0 \left( {1 + 0.2\frac{{\left| j \right|}} {B}} \right)$$
((30))

is an effective turbulent diffusivity. The first term is a uniform background value η0 = 45 km2 s-1 and the second term is an enhancement in regions of strong current density, introduced to limit the twist in helical flux ropes to about one turn, as observed in solar filaments.

Rather than solving the full MHD equations for the velocity v, which is not computationally feasible over long timescales, the quasi-static evolution of the coronal magnetic field is approximated using magneto-frictional relaxation (Yang et al., 1986), setting

$$v = \frac{1} {\nu }\frac{{j \times B}} {{{\rm B}^2 }}.$$
((31))

This artificial velocity causes the magnetic field to relax towards a force-free configuration: it may be shown that the total magnetic energy decreases monotonically until j × B = 0. Here, the magneto-frictional relaxation is applied concurrently with the photospheric driving so, in practice, a dynamical equilibrium between the two is reached. At the outer boundary r = 2.5 R, Mackay and van Ballegooijen (2006) introduced an imposed radial outflow, rather than setting B θ = B φ = 0 exactly. This allows magnetic flux ropes to be ejected after they lose equilibrium (simulating CMEs on the real Sun), and simulates the effect of the solar wind in opening up magnetic field lines radially above this height.

Figure 10b shows an example nonlinear force-free field using this model, seen after 100 days of evolution. Figure 10a shows a PFSS extrapolation from the same photospheric B r . distribution. The coronal field in the non-potential model is significantly different, comprising highly twisted flux ropes, slightly sheared coronal arcades, and near potential open field lines. Since this model can be run for extended periods without resetting the coronal field, it can be used to study long-term helicity transport across the solar surface from low to high latitudes (Yeates et al., 2008a).

Figure 10
figure 10

(a) Example of a Potential Field Source Surface extrapolation. (b) Example of a non-potential coronal field produced by the global coronal evolution model of Mackay and van Ballegooijen (2006) and Yeates et al. (2008b). The grey-scale image shows the radial field at the photosphere and the thin lines the coronal field lines. Image reproduced by permission from Yeates et al. (2010b), copyright by AGU.

There are several parameters in the model. As in the PFSS model the location of the upper boundary Rss is arbitrary, although it has less influence since the non-potential magnetic field strength falls off more slowly with radius than that of the potential field. In the magneto-frictional evolution, the turbulent diffusivity η0 and friction coefficient v are arbitrary, and must be calibrated by comparison with observed structures or timescales for flux ropes to form or lose equilibrium. Finally a 3D model is required for newly-emerging active regions: the simple analytical bipoles of the existing simulations could in fact be replaced with more detailed extrapolations from observed photospheric fields in active regions.

3.3 Magnetohydrostatic models

Models in this category are based on analytical solutions to the full magnetohydrostatic (MHS) equations,

$$j \times B - \nabla p - \rho \nabla \psi = 0,$$
((32))
$$\nabla \times B = \mu _0 j,$$
((33))
$$\nabla \cdot B = 0,$$
((34))

where for the coronal application ψ = -GM/r is taken to be the gravitational potential. Three-dimensional solutions that are general enough to accept arbitrary B r (R, θ, φ) as input have been developed by Bogdan and Low (1986) and further by Neukirch (1995). The basic idea is to choose a particular functional form of j, and use the freedom to choose ρ and p to obtain an analytical solution. Note that, while the distributions of ρ, p, and B are self-consistent in these solutions, the 3D forms of ρ and p are prescribed in the solution process. In particular, they cannot be constrained a priori to satisfy any particular equation of state or energy equation. A realistic treatment of the thermodynamics of the plasma requires a full MHD model (Section 3.4).

The solutions that have been used to extrapolate photospheric data have a current density of the form

$$j = \alpha B + \xi (r)\nabla (r \cdot B) \times e_r,$$
((35))

where

$$\xi (r) = \frac{1} {{r^2 }} - \frac{1} {{(r + a)^2 }}$$
((36))

and α, a are constant parameters. Thus, the current comprises a field-aligned part and a part perpendicular to gravity. For this form of j, Neukirch (1995) shows that the solution takes the form

$$B_\theta = \sum\limits_{l = 1}^\infty {\sum\limits_{m = - 1}^l {\sum\limits_{j = 1}^2 {A_{lm}^{(j)} \frac{{l(l + 1)}} {r}u_l^{(j)} (r)P_l^m (\cos \theta )e^{im\varphi } } } } , $$
((37))
$$B_\theta = \sum\limits_{l = 1}^\infty {\sum\limits_{m = - 1}^l {\sum\limits_{j = 1}^2 {A_{lm}^{(j)} \left[ {\frac{1} {r}\frac{{d(ru_l^{(j)} (r))}} {{dr}}\frac{{dP_l^m (\cos \theta )}} {{d\theta }} + \frac{{\mu _0 \alpha mi}} {{\sin \theta }}u_l^{(j)} (r)P_l^m (\cos \theta )} \right]e^{im\varphi } } } } , $$
((38))
$$B_\varphi = \sum\limits_{l = 1}^\infty {\sum\limits_{m = - 1}^l {\sum\limits_{j = 1}^2 {A_{lm}^{(j)} \left[ {\frac{{im}} {{r\sin \theta }}\frac{{d(ru_l^{(j)} (r))}} {{dr}}P_l^m (\cos \theta ) - \mu _0 \alpha u_l^{(j)} (r)\frac{{dP_l^m (\cos \theta )}} {{d\theta }}} \right]e^{im\varphi } } } } , $$
((39))

where

$$u_l^{(1)} (r) = \frac{{\sqrt {r + a} }} {r}J_{l + 1/2} \left( {\alpha (r + a)} \right),$$
((40))
$$u_l^{(2)} (r) = \frac{{\sqrt {r + a} }} {r}N_{l + 1/2} \left( {\alpha (r + a)} \right).$$
((41))

Here, Jl+1/2 and Nl+1/2 are Bessel functions of the first and second kinds, respectively. The coefficients A (j) lm are determined from the boundary conditions, as in the PFSS model. The plasma pressure and density then take the forms

$$p(r,\theta ,\varphi ) = p_0 (r) - \frac{{\xi (r)}}, {2}(r \cdot B)^2 $$
((42))
$$\rho (r,\theta ,\varphi ) = \rho _0 (r) + \frac{{r^2 }} {{GM}}\left( {\frac{1} {2}\frac{{d\xi (r)}} {{dr}}(r \cdot B)^2 + r\xi (r)B \cdot \nabla (r \cdot B)} \right). $$
((43))

The functions p0(r) and ρ0(r) describe a spherically symmetric background atmosphere satisfying - ∇p0 - ρ0ψ = 0, and may be freely chosen.

Solutions of this form have been applied to the coronal extrapolation problem by Zhao et al. (2000), Rudenko (2001), and Ruan et al. (2008).Footnote 6 These solutions are of exterior type, i.e., with |B| → 0 as r → ∞. For α ≠ 0, the latter condition follows from the properties of both types of Bessel function.

There are two free parameters in the solution that may be varied to best fit observations: a and α. Broadly speaking, the effect of increasing a is to inflate/expand the magnetic field, while the effect of increasing |α| is to twist/shear the magnetic field. This is illustrated in Figure 11. Note that too large a value of α would lead to zeros of the Bessel functions falling within the computational domain, creating magnetic islands that are unphysical for the solar corona (Neukirch, 1995). The (unbounded) potential field solution j = 0 is recovered for α = a = 0, while if a = 0 then j = αB and we recover the linear force-free field case. If α = 0 then the current is purely horizontal and the solution reduces to Case III of Bogdan and Low (1986). Gibson and Bagenal (1995) applied this earlier solution to the Solar Minimum corona, although they showed that it was not possible to match both the density distribution in the corona and the photospheric magnetic field to observations. This problem is likely to be exacerbated at Solar Maximum. Similarly, Ruan et al. (2008) found that the strongest density perturbation in this model appears in active regions in the low corona. Preventing the density from becoming negative can require an unrealistically large background density ρ0 at these radii, particularly for large values of the parameter a.

Figure 11
figure 11

Effect of the parameters a and α in the MHS solution for a simple dipolar photospheric magnetic field. Closed magnetic field lines expand when a is changed (left column), while the shear increases as the field-aligned current parameter α is increased (middle column). If both current systems act together (right column) then the field expands and becomes twisted. Image reproduced by permission from Zhao et al. (2000), copyright by AAS.

Another limitation of this magnetic field solution is that the magnetic energy is unbounded, since the Bessel functions decay too slowly as r → ∞. But, as pointed out by Zhao et al. (2000), the model is applicable only up to the cusp points of streamers, above which the solar wind outflow must be taken into account. So this problem is irrelevant in practical applications. For the α = 0 solution of Bogdan and Low (1986), Zhao and Hoeksema (1994) showed how the model can be extended to larger radii by adding the “current sheet” extension of Schatten (1971) (described in Section 3.1). Instead of using an “exterior” solution in the external region they introduce an outer source surface boundary at Rss ≈ 14R, corresponding to the Alfvén critical point. The resulting model better matches the shape of coronal structures and the observed IMF (Zhao and Hoeksema, 1995) and is often termed the Current Sheet Source Surface (CSSS) model, though the authors term it HCCSSS: “Horizontal Current Current Sheet Source Surface”, to explicitly distinguish it from the Schatten (1971) model using potential fields. This CSSS model has also been applied by Schüssler and Baumann (2006) to model the Sun’s open magnetic flux (Section 4.1).

Finally, we note that Wiegelmann et al. (2007) have extended the numerical optimisation method (Section 3.2) to magnetohydrostatic equilibria, demonstrating that it reproduces the analytical solution of Neukirch (1995). This numerical technique offers the possibility of more realistic pressure and density profiles compared to the analytical solutions, although there is the problem that boundary conditions on these quantities must be specified.

3.4 Full magnetohydrodynamic models

In recent years, a significant advance has been made in the construction of realistic 3D global MHD models (Riley et al., 2006; DeVore and Antiochos, 2008; Lionello et al., 2009; Downs et al., 2010; Feng et al., 2012). Such models are required to give a self-consistent description of the interaction between the magnetic field and the plasma in the Sun’s atmosphere. They allow for comparison with observed plasma emission (Section 4.4), and enable more consistent modeling of the solar wind, so that the simulation domain may extend far out into the heliosphere (e.g., Riley et al., 2011; Tóth et al., 2012; Feng et al., 2012). On the other hand, additional boundary conditions are required (for example, on density or temperature), and the problem of non-uniqueness of solutions found in nonlinear force-free field models continues to apply here. While these models are sometimes used to simulate eruptive phenomena such as coronal mass ejections, in this review we consider only their application to non-eruptive phenomena. For applications to eruptive phenomena the reader is directed to the reviews of Forbes et al. (2006) and Chen (2011).

An example of the resistive MHD equations used in the paper of Lionello et al. (2009) for constructing a steady-state solution are:

$$\frac{{\partial \rho }} {{\partial t}} + \nabla .(\rho v) = 0,$$
((44))
$$\rho \left( {\frac{{\partial v}} {{\partial t}} + v.\nabla v} \right) = \frac{1} {c}J \times B - \nabla (p + p_w ) + \rho g + \nabla .(\nu \rho \nabla v), $$
((45))
$$\frac{1} {{\gamma - 1}}\left( {\frac{{\partial T}} {{\partial t}} + v.\nabla T} \right) = - T\nabla .v + \frac{m} {{k_\rho }}S,$$
((46))
$$S = - \nabla .q - n_e n_p Q(t) + H_{ch}, $$
((47))
$$\nabla \times B = \frac{{4\pi }} {c}J,$$
((48))
$$\nabla \times E = - \frac{1} {c}\frac{{\partial B}} {{\partial t}},$$
((49))
$$E + \frac{{v \times B}} {c} = \mu J,$$
((50))

where B, J,E, ρ, v, p, T, g, η, v,Q(t), n e , n p , γ = 5/3,Hch, q, p w are the magnetic field, electric current density, electric field, plasma density, velocity, pressure, temperature, gravitational acceleration, resistivity, kinematic viscosity, radiative losses, electron and proton number densities, polytropic index, coronal heating term, heat flux, and wave pressure, respectively. Different formulations of the MHD equations may be seen in the papers of DeVore and Antiochos (2008), Downs et al. (2010), and Feng et al. (2012) where key differences are the inclusion of resistive or viscous terms and the use of adiabatic or non-adiabatic energy equations.

To date, non-eruptive 3D global MHD simulations have been used to model the solar corona through two distinct forms of simulation. Firstly, there is the construction of steady state coronal solutions from fixed photospheric boundary conditions (Riley et al., 2006; Vásquez et al., 2008; Lionello et al., 2009; Downs et al., 2010). To compute these solutions, the system is initialised by (i) specifying the photospheric distribution of flux (often from observations), (ii) constructing an initial potential magnetic field and, finally, (iii) superimposing a spherically symmetric solar wind solution. Equations (44) – (50) (or their equivalents) are then integrated in time until a new equilibrium is found. A key emphasis of this research is the direct comparison of the resulting coronal field and plasma emission with that seen in observations (Section 4.4). In the paper of Vásquez et al. (2008) a quantitative comparison of two global MHD models (Stanford: Hayashi, 2005, and Michigan: Cohen et al., 2007), with coronal densities determined through rotational tomography was carried out. In general the models reproduced a realistic density variation at low latitudes and below 3.5R, however, both had problems reproducing the correct density in the polar regions. In contrast to this construction of steady state MHD solutions, advances in computing power have recently enabled global non-eruptive MHD simulations with time-dependent photospheric boundary conditions. These boundary conditions have been specified both in an idealised form (e.g., Lionello et al., 2005; Edmondson et al., 2010) and from synoptic magnetograms (Linker et al., 2011). An initial application of these models has been to simulate the Sun’s open flux and coronal holes: see Section 4.2.

To quantify the difference between the magnetic field produced in global MHD simulations and PFSS extrapolations, Riley et al. (2006) carried out a direct comparison of the two techniques. The results from each technique are compared in Figure 12 for periods of low activity (top) and high activity (bottom). An important limitation of this study was that, for both the PFSS model and MHD model, only the radial magnetic field distribution at the photosphere was specified in order to construct the coronal field. Under this assumption, a good agreement was found between the two fields where the only differences were (i) slightly more open flux and larger coronal holes in the MHD simulation, and (ii) longer field lines and more realistic cusp like structures in the MHD simulations. While the two set of simulations closely agree, as discussed by the authors it is important to note that if vector magnetic field measurements are applied on a global scale as an additional lower boundary condition, then significant differences may result.

Figure 12
figure 12

Comparison of PFSS extrapolations (right column) with steady state MHD solutions (left column) for two Carrington rotations. Image reproduced by permission from Figure 5 of Riley et al. (2006), copyright by AAS.

While global MHD models have mostly been applied to the Sun, Cohen et al. (2010) applied the models to observations from AB Doradus. Through specifying the lower boundary condition from a ZDI magnetic map and also taking into account the rotation (but not differential rotation) of the star the authors showed that the magnetic structure deduced from the MHD simulations was very different from that deduced from PFSS models. Due to the rapid rotation a strong azimuthal field component is created. They also considered a number of simulations where the base coronal density was varied and showed that the resulting mass and angular momentum loss, which are important for the spin-down of such stars, may be orders of magnitude higher than found for the Sun.

4 Application of Magnetic Flux Transport and Coronal Models

In this section, we survey the combined application of magnetic flux transport (Section 2.2) and coronal models (Section 3) to describe a number of phenomena on the Sun. These include the Sun’s open magnetic flux (Section 4.1), coronal holes (Section 4.2), and the hemispheric pattern of solar filaments (Section 4.3). Finally, the application of global MHD models to predict the plasma emission from the Sun is described in Section 4.4. It should be noted that when the magnetic flux transport and coronal models are combined, two distinct modeling techniques are applied, each of which have their own limitations. It is therefore important to consider how the limitations of each affect the results. In addition to this, some models also use observations as input. These observations also have many practical limitations related to calibration, resolution, instrumental effects and data reduction effects. It is currently beyond the scope of the present article to consider the full consequences of these combined uncertainties and limitations. However, all readers should keep this in mind when assessing results and outcomes. In each of the following sections a brief description of the main properties of each coronal phenomenon is first given, before discussing the application of the modeling techniques.

4.1 Open flux

The Sun’s open magnetic flux is the part of the Sun’s large-scale magnetic field that extends from the solar surface out into interplanetary space, where it forms what is known as the Interplanetary Magnetic Field (IMF). Understanding the origin and variation of the open flux is important as:

  • — It is the origin of the high speed solar wind and IMF, both of which directly interact with the Earth’s magnetosphere.

  • — Irregularities or curvatures in the IMF modulate the flux of high energy Cosmic Rays that impact the Earth’s upper atmosphere. Such impacts produce 14C (Stuiver and Quay, 1980) and 10Be (Beer et al., 1990) isotopes which are found in tree rings and ice cores, respectively. These isotopes may be used as historical data sets of solar activity since there is an anticorrelation between the strength of open flux and the abundance of the isotopes.

  • — The open flux has been correlated with cloud cover variations (Svensmark and Friis-Christensen, 1997; Svensmark, 1998), however, any relation to global warming is very uncertain.

  • — Strong correlations are found between the Sun’s open flux and total solar irradiance (TSI, Lockwood, 2002). Therefore historical data such as 14C and 10Be, which allow us to deduce the amount of open flux, may also be used as a proxy for past TSI variations.

Our present day understanding of the Sun’s open flux comes mainly from an important result from the Ulysses mission: the magnitude of the radial IMF component in the heliosphere is independent of latitude (Balogh et al., 1995). Thus, magnetometer measurements at a single location at 1 AU may be used to deduce the total open flux simply by multiplying by the surface area of a sphere. In Figure 13, the variation of open flux relative to sunspot number can be seen for the last 3.5 solar cycles. Over a solar cycle the open flux varies at most by a factor of two, even though the surface flux shows a much larger variation. However, this variation is not regular from one cycle to the next. It is clear from the graph that Cycles 21 and 22 have a much larger variation in open flux than Cycles 20 and 23, indicating a complex relationship between the surface and open flux. A key property of the open flux is that it slightly lags behind the variation in sunspot number and peaks 1 – 2 years after cycle maximum. The peak time of open flux therefore occurs around the same time as polar field reversal, indicating that contributions from both polar coronal holes and low-latitude coronal holes must be important.

Figure 13
figure 13

Graph of open flux variation (top) and sunspot number (bottom) for cycles 20 – 23. Image reproduced by permission from Lockwood et al. (2004), copyright by EGU.

While both direct and indirect observations exist for the Sun’s open flux, such measurements cannot be made for other stars. Due to this, theoretical models which predict the magnitude and spatial distribution of open flux are used. Within the stellar context, the distribution of open flux with latitude (McIvor et al., 2006) plays a key role in determining the mass and angular momentum loss (Cohen et al., 2007) and subsequently the spin down of stars (Weber and Davis Jr, 1967; Schrijver et al., 2003; Holzwarth and Jardine, 2007).

Over the last 20 years a variety of techniques has been developed to model the origin and variation of the Sun’s open flux. One category are the semi-empirical magnitude variation models. These follow only the total open flux, and are driven by observational data. Different models use either geomagnetic data from Earth (Lockwood et al., 1999, 2009), sunspot numbers (Solanki et al., 2000, 2002), or coronal mass ejection rates (Owens et al., 2011). Such models have been used to extrapolate the open flux back as far as 1610 (in the case of sunspot numbers), although the results remain uncertain. However, these models consider only the total integrated open flux, not its spatial distribution and origin on the Sun. To describe the latter, coupled photospheric and coronal magnetic field models have been applied. The photospheric field is specified either from synoptic magnetic field observations or from magnetic flux transport simulations (Section 2.2), and a wide variety of coronal models have been used. These range from potential field source surface (PFSS, Section 3.1) models to current sheet source surface (CSSS, Section 3.1) models and more recently to nonlinear force-free (NLFF, Section 3.2) models. In all of these coronal models, field lines reaching the upper boundary are deemed to be open, and the total (unsigned) magnetic flux through this boundary represents the open flux.

Wang and Sheeley Jr (2002) combined synoptic magnetograms and PFSS models to compute the open flux from 1971 to 1998. The synoptic magnetograms originated from either Wilcox Solar Observatory (WSO, 1976 – 1995) or Mount Wilson Observatory (MWO, 1971 – 1976, 1995 – 1998). The authors found that to reproduce a good agreement to IMF field measurements at 1 AU, they had to multiply the magnetograms by a strong latitude-dependent correction factor (4.5 – 2.5 sin2 λ). This factor, used to correct for saturation effects in the magnetograph observations, significantly enhanced the low latitude fields compared to the high latitude fields but gave the desired result.

The use of such a latitude dependent correction factor was recently questioned by Riley (2007), who pointed out that the correction factor used by Wang and Sheeley Jr (2002) was only correct for use on MWO data and did not apply to WSO data. Instead WSO data, which made up the majority of the magnetogram time series used, should have a constant correction factor of 1.85 (Svalgaard, 2006). On applying this appropriate correction and repeating the calculation, Riley (2007) showed that PFSS models gave a poor fit to observed IMF data (see Figure 4 of Riley, 2007). To resolve the difference, Riley (2007) put forward an alternative explanation. He assumed that the open flux has two contributions: (i) a variable background contribution based on the photospheric field distribution on the Sun at any time, such as that obtained from the PFSS model, and (ii) a short term enhancement due to interplanetary coronal mass ejections (ICMEs), for which he gave an order-of-magnitude estimate. Through combining the two a good agreement to IMF field observations was found.

As an alternative to synoptic magnetic field observations, magnetic flux transport models have been widely used to provide the lower boundary condition in studies of the origin and variation of the Sun’s open flux. A key element in using magnetic flux transport simulations is that the distribution and strength of the high latitude field are determined by the properties and subsequent evolution of the bipoles which emerge at lower latitudes. Hence, assuming the correct input and advection parameters, the model will produce a better estimate of the polar field strength compared to that found in magnetogram observations which suffer from severe line-of-sight effects above 60° latitude. Initial studies (using a PFSS coronal model) considered the variation in open flux from a single bipole (Wang et al., 2000; Mackay et al., 2002a) as it was advected across the solar surface. Through doing so, the combined effects of latitude of emergence, Joy’s Law and differential rotation on the open flux was quantified. Later studies extended these simulations to include the full solar cycle, but obtained conflicting results.

To start with, Mackay et al. (2002b), using idealised simulations and commonly used transport parameters, found that such parameters when combined with PFSS models produced an incorrect time-variation of open flux. The open flux peaked at cycle minimum, completely out of phase with the solar cycle and inconsistent with the observed 1 – 2 year lag behind solar maximum. The authors attributed this to the fact that in PFSS models only the lowest order harmonics (see Section 3.1), which are strongest at cycle minimum and weakest at maximum, significantly contribute to the open flux. They concluded that the only way to change this was to include non-potential effects (see below) which would increase the open flux contribution from higher order harmonics during cycle maximum. In response, Wang et al. (2002b) repeated the simulation and found similar results when using bipole fluxes and parameters deduced from observations. However, they were able to obtain the correct variation in open flux by both trebling all observed bipole fluxes and increasing the rate of meridional flow to 25 m s-1 (dash-dot-dot-dot line in Figure 3b). These changes had the effect of significantly enhancing the low latitude field at cycle maximum but weakening the polar fields at cycle minimum. As a result the correct open flux variation was found. Opposing such a strong variation in the input parameters, Schüssler and Baumann (2006) put forward another possibility. Through modifying the magnetic flux transport model to include a radial decay term for B r (Section 2.2), and changing the coronal model to a CSSS model with Rcp ∼ 1.7R and Rss ∼ 10R (Section 3.3), they found that the correct open flux variation could be obtained, but only if new bipole tilt angles were decreased from 0.5λ to 0.15λ. Later studies by Cameron et al. (2010) using the same technique showed that the radial decay term was not required if the tilt angles of the bipoles varied from one cycle to the next, with stronger cycles having weaker tilt angles.

The discussion above shows that, similar to the polar field strength (Section 2.2), the correct variation of the Sun’s open flux may be obtained through a variety of methods. These include variations in the bipole tilt angles, rates of meridional flow, and different coronal models. More recently, Yeates et al. (2010b) showed that standard values for bipole tilt angles and meridional flow may produce the correct variation of open flux if a more realistic physics-based coronal model is applied (see Section 3.2.3). By allowing for electric currents in the corona, a better agreement to the IMF measurements can be found compared to those from PFSS models. Figure 14a compares various PFSS extrapolations using different magnetogram data (coloured lines) to the measured IMF field (grey line). In this graph the discrepancy between all potential field models and the IMF is particularly apparent around cycle maximum. To resolve this, Yeates et al. (2010b) ran the global non-potential model described in Section 3.2.3 over four distinct 6-month periods (labeled A-D in Figure 14a). Two periods were during low solar activity (A and D) and two during high solar activity (B and C). The results of the simulation for period B can be seen in Figure 14b. In the plot, the dashed lines denote open flux from various PFSS extrapolations, the grey line the observed IMF field strength and the black solid line the open flux from the non-potential global simulation. This clearly gives a much better agreement in absolute magnitude terms compared to the PFSS models. Yeates et al. (2010b) deduced that the open flux in their model has three main sources. The first is a background level due to the location of the flux sources. This is the component that may be captured by PFSS models, or by steady-state MHD models initialised with PFSS extrapolations (Feng et al., 2012). The second is an enhancement due to electric currents, which results in an inflation of the magnetic field, visible in Figure 14b as an initial steady increase of the open flux curve during the first month to a higher base level. This inflation is the result of large scale flows such as differential rotation and meridional flows along with flux emergence, reconfiguring the coronal field. A similar inflation due to electric currents enhances the open flux in the CSSS model, although there the pattern of currents is arbitrarily imposed, rather than arising physically. Finally, there is a sporadic component to the open flux as a result of flux rope ejections representing CMEs similar to that proposed by Riley (2007). While this model gives one explanation for the open flux shortfall in PFSS models, we note that alternative explanations have been put forward by a variety of authors (Lockwood et al., 2009; Fisk and Zurbuchen, 2006). Again, although models have been successful in explaining the variation of the open flux, there is still an uncertainty on the true physics required. Free parameters that presently allow for multiple models must be constrained by future observations. At the same time, it is evident from using existing observations that care must be taken in their use. In particular, the effect on the modeling results of practical limitations such as calibration, resolution, instrumental effects, and data reduction effects must be quantified.

Figure 14
figure 14

Graph of (a) IMF variation (grey line) and PFSS approximations to the IMF field (collared lines). (b) Non-potential open flux estimate (thick black line) along with IMF variation (grey line) and PFSS approximations to the IMF (coloured lines) over a 6 month period. Image reproduced by permission from Yeates et al. (2010b), copyright by AGU.

Although they include the spatial distribution of open flux, the models described above have mainly focused on the total amount of open flux, because this may be compared with in situ measurements. However, these models can also be used to consider how open flux evolves across the solar surface, and its relationship to coronal holes. Presently, there are two opposing views on the nature of the evolution of open flux. The first represents that of Fisk and Schwadron (2001) who suggest that open flux may diffuse across the solar surface through interchange reconnection. In this model, the open flux may propagate into and through closed field regions. The second view represents that of Antiochos et al. (2007) who postulate that open flux can never be isolated and thus cannot propagate through closed field regions. Results of global MHD simulations testing these ideas are discussed in the next section.

4.2 Coronal holes

Coronal holes are low density regions of the corona which appear dark in X-rays or EUV, or as light, blurred patches in He i 10830 Å (Cranmer, 2009). It has been known since the Skylab era of the early 1970s that coronal holes are associated with regions of open magnetic field lines (Section 4.1), along which plasma escapes into the heliosphere (Zirker, 1977). They therefore play an important role in understanding the sources of different solar wind streams and the relationship between the Sun and the heliosphere.

The basic link between coronal holes and open magnetic field was found using PFSS extrapolations from synoptic magnetogram observations. While the open field regions in such extrapolations produce pleasing agreement with observed coronal hole boundaries (Levine, 1982; Wang et al., 1996; Neugebauer et al., 1998; Luhmann et al., 2002), our understanding has been enhanced by coupling coronal magnetic models with photospheric simulations, in the following ways:

  1. 1.

    Photospheric simulations can be used to produce a more accurate picture of the global magnetic field by assimilating observed data and following the evolution on timescales shorter than the 27 days of the synoptic charts. Schrijver and DeRosa (2003) found that some open field regions in such a model were quite stable for weeks or even months, whereas others fluctuated rapidly.

  2. 2.

    The photospheric evolution can be imposed or controlled, in order to test hypotheses for the behaviour of coronal holes, and to try to understand the origin of their observed properties.

We outline here the studies that fall into the second category. Firstly, a series of papers by Wang, Sheeley, and co-workers, coupling surface flux transport to PFSS extrapolations, has significantly added to our understanding of coronal hole dynamics. By isolating different terms in the flux transport model, Wang and Sheeley Jr (1990) demonstrate how each of the terms in the flux transport model play their roles in producing the global pattern of coronal holes. Supergranular diffusion spreads out active region flux so as to create unipolar areas containing open flux, and also facilitates the build up of trailing polarity holes in each hemisphere after Solar Maximum by annihilating leading polarity flux at the equator. Differential rotation helps to accelerate the formation of axisymmetric polar holes by symmetrising the active region flux distribution. Meridional flow (i) concentrates flux at the poles - preventing the polar holes from spreading to lower latitudes, (ii) hastens the decay of low-latitude holes by transporting them to mid-latitudes where differential rotation is more efficient, and (iii) impedes cancellation across the equator, thus reducing the flux imbalance in each hemisphere.

A striking observation arising from Skylab was that coronal holes rotate differently from the underlying photosphere. This is exemplified by the Skylab “Coronal Hole 1”. Nash et al. (1988) were able to reproduce and explain the behaviour of this hole using the coupled flux transport and PFSS model (see also Wang and Sheeley Jr, 1993). The mechanism is illustrated by their model of a single bipolar active region in a dipolar background field (Wang and Sheeley Jr, 1990; Wang et al., 1996, reproduced here in Figure 15). There are two important points. The first is that coronal hole boundaries in the PFSS model are determined by the global magnetic structure of the coronal field. This is because the source surface field depends only on the lowest order spherical harmonics. In particular, the open field regions in Figure 15 are not simply a superposition of those that would be obtained from the bipole and background fields individually (Wang et al., 1996). Hence, it need not be surprising that coronal holes are observed to rotate differently from the footpoints of individual field lines. The second point is that the rotation rate is determined by the non-axisymmetric component only. In other words, the “coronal hole extension” in Figure 15 must rotate approximately with the bipole giving rise to it. It follows that the rotation rate of the coronal hole matches that of the photospheric latitude where the bipole is located, not necessarily the latitude where the open field lines have their footpoints. This explains why the northern and southern coronal hole extensions in Figure 15 behave differently: the bipole flux causing the northern hole extension is located in a narrow band close to the equator, so that the hole rotates rigidly with the equatorial 27-day period. On the other hand, the flux causing the southern hole extension is spread between 0° and 40° latitude, so that this hole is considerably sheared by differential rotation. The same idea explains why rigidly rotating coronal holes are more prevalent in the declining phase of the solar cycle: during this phase, the non-axisymmetric flux is concentrated more at lower latitudes.

Figure 15
figure 15

Effect of differential rotation on open field regions for a bipole (black/white) located in a dipolar background field (red/dark blue). The open region (“coronal hole”) in the North (yellow) rotates almost rigidly, while that in the South (light blue) is sheared. Image reproduced by permission from Wang et al. (1996), copyright by AAAS.

More recently, new insights into coronal hole evolution have started to come from coupling global MHD models (Section 3.4) to surface flux transport. For example, an important implication of the PFSS model for coronal magnetic fields is that continual magnetic reconnection is necessary to maintain the current-free field (Wang and Sheeley Jr, 2004). This has been tested in a global MHD model by Lionello et al. (2005), who applied differential rotation to a configuration consisting of a single bipole in a dipolar background field. They confirm the results of Wang et al. (1996) that the coronal hole extension rotates nearly rigidly without significant change, even when the surface field is significantly sheared by differential rotation. The dominant reconnection process was found to be interchange reconnection, with continual reconnection opening field on the eastern hole boundary and closing it on the western boundary. In some cases, when field line footpoints pass over the coronal hole boundaries multiple times, multiple openings and closing occur. However, it should be noted that these simulations used a uniform resistivity orders of magnitude greater than that in the real corona, so that reconnection may be less efficient in reality. Fisk (1996) has proposed an alternative scenario where open field lines continually circulate in latitude and longitude as their footpoints rotate differentially. This process was applied to explain the origin of high latitude Solar Energetic Particles (SEPs) seen by Ulysses (Keppler et al., 1995; Maclennan et al., 1995), which are thought to originate in low latitude corotating interacting regions (CIRs, McDonald et al., 1976). However, in a test with their global MHD model and an initial tilted dipole field, Lionello et al. (2006) found that although field lines did move in latitude, the coronal hole boundaries rotated in a manner consistent with the extrapolation models, and not perfectly rigidly as predicted by Fisk (1996). This meant that the latitudinal excursion of around 25° found in the simulations was insufficient to explain that required by Fisk (1996).

Recently, Linker et al. (2011) have applied their global MHD model to test another hypothesis, namely, that the total amount of open flux is conserved. This has been proposed by Fisk and co-workers based on heliospheric observations (Fisk and Schwadron, 2001). They propose that open flux is transported by interchange reconnection (between open and closed field lines), but not created or destroyed through other forms of reconnection. A feature of this model is that open field lines may be transported through closed field regions in a random manner during the reversal of the polar field. In their simulations, Linker et al. (2011) start with a synoptic magnetogram for CR1913, but add two small bipoles of varying orientation inside one of the large equatorward coronal hole extensions. They then advect the bipoles through the coronal hole from the open to closed field regions. The idea is to test (a) whether open flux associated with the bipoles can be transported out of the coronal hole and into the closed field region (as required by the interchange model; Fisk and Zurbuchen, 2006), and (b) whether isolated coronal holes can exist in each hemisphere. The results (Figure 16) seem to contradict the requirements of the interchange model, and to uphold the conjecture of Antiochos et al. (2007), namely that what appear to be isolated coronal holes in each hemisphere, are in fact, always connected (by very thin corridors) of open flux to the polar holes in each hemisphere. These thin corridors may intersect with strong field regions that form part of the streamer belt. Thus, they may have a non-negligible contribution to the amount of open flux, which may help resolve the discrepancy between the amount of open magnetic flux obtained in models and that directly observed through IMF measurements (Section 4.1). Following this, Titov et al. (2011) constructed a topological model for the open flux corridors that connect coronal holes. Further studies by Edmondson et al. (2009, 2010) considered idealised case studies of how magnetic reconnection affects the boundaries of coronal holes when either a bipole is (i) advected across the coronal boundary (Edmondson et al., 2010) or (ii) stressed by rotational motions (Edmondson et al., 2009). In line with the conjecture of Antiochos et al. (2007) and quasi-steady model results of Linker et al. (2011), the authors see a smooth transition of the open flux where the open field regions are never isolated. Such evolution of coronal hole boundaries along with the translation and opening/closing of magnetic field lines has important consequences for the solar wind. A full discussion of this is beyond the scope of the present review but details can be found in Cranmer (2009) and Antiochos et al. (2012).

Figure 16
figure 16

Simulation results from the paper of Linker et al. (2011) where a bipole is advected from an open field region (shaded area on the solar surface) to a closed field region across the boundary of a coronal hole. Red denotes positive photospheric flux and blue negative flux. The green lines denote (a) initially open field lines of the bipole which then successively close down (b) - (d) as the bipole is advected across the coronal hole boundary. Image reproduced by permission, copyright by AAS.

4.3 Hemispheric pattern of solar filaments

Solar filaments (a.k.a. prominences) are cool dense regions of plasma (T ∼ 10, 000 K) suspended in the hotter (T ∼ 106 K) solar corona. Coronal magnetic fields are key to the existence of solar filaments, where they provide both support against gravity and insulation of the cool material from the surrounding hot corona. The main observational properties of filaments and the theoretical models used to describe them are described in the reviews by Labrosse et al. (2010) and Mackay et al. (2010).

In recent years, solar filaments and their birth grounds called filament channels (Gaizauskas, 1998), which always overlie photospheric polarity inversion lines (PILs), have been classified by the orientation of their main axial magnetic field. This orientation, named the filament chirality (Martin et al., 1994) may take one of two forms: dextral or sinistral. Dextral/sinistral filaments have an axial magnetic field that points to the right/left when the main axis of the filament is viewed from the positive polarity side of the PIL. In force-free field models (Aulanier and Démoulin, 1998; Mackay et al., 1999; van Ballegooijen et al., 2000; Mackay and van Ballegooijen, 2005) this chirality is directly related to the dominant sign of magnetic helicity that is contained within the filament and filament channel. A dextral filament will contain dominantly negative helicity, a sinistral filament positive helicity. As filaments and their channels form over a wide range of latitudes on the Sun, ranging from the active belts to the polar crown, they may be regarded as useful indicators of sheared non-potential fields and magnetic helicity within the solar corona.

A surprising feature of the chirality of filaments is that it displays a large-scale hemispheric pattern: dextral/sinistral filaments dominate in the northern/southern hemispheres, respectively (Martin et al., 1994; Zirker et al., 1997; Pevtsov et al., 2003; Yeates et al., 2007). This pattern is intriguing as it is exactly opposite to that expected from differential rotation acting on a North-South coronal arcade. Although dextral/sinistral filaments dominate in the northern/southern hemisphere, observations show that exceptions to this pattern do occur. Therefore any model which tries to explain the formation of filaments and their dominant axial magnetic fields must explain not only the origin of this hemispheric pattern but also why exceptions to it arise.

Since filament chirality is directly related to magnetic helicity (dextral ∼ negative, sinistral ∼ positive) the formation and transport of filaments across the solar surface is an indication of the global pattern of magnetic helicity transport on the Sun (Yeates et al., 2008a), a key feature in explaining many eruptive phenomena. Static extrapolation techniques (including MHD models of the relaxation type) cannot study the generation nor transport of this helicity across the surface of the Sun. This can only be achieved with models that couple the evolution of both photospheric and coronal fields over long periods of time.

The first attempt to explain the hemispheric pattern of filaments using global magnetic flux transport models was carried out by van Ballegooijen et al. (1998). By using observed magnetic flux distributions and initial coronal fields which were potential, they simulated the non-equilibrium evolution of the coronal field. They found that the flux transport effects acting alone on potential fields create approximately equal numbers of dextral and sinistral channels in each hemisphere, in contradiction with observations. Following this, Mackay and van Ballegooijen (2005) re-considered the origin of the hemispheric pattern through combined flux transport and magneto-frictional relaxation simulations (Section 3.2.3), where the coronal field responds to surface motions by relaxing to a nonlinear force-free equilibrium. In the simulations the authors considered an idealised setup of initially isolated bipoles. Therefore, in contrast to the study of van Ballegooijen et al. (1998), they did not use initial potential fields. The authors demonstrate that the hemispheric pattern of filaments may be explained through the observational properties of newly-emerging bipoles such as (i) their dominant tilt angles (-10°:30°, Wang and Sheeley Jr, 1989), and (ii) the dominant helicity with which they emerge in each hemisphere (Pevtsov et al., 1995). A key feature of these simulations was that the possible occurrence of exceptions to the hemispheric pattern was quantified for the first time, arising from large positive bipole tilt angles and minority helicity.

In a more conclusive study, the results of Mackay and van Ballegooijen (2005) have been tested through a direct comparison between theory and observations by Yeates et al. (2007, 2008b). First, they used Hα observations from BBSO over a 6 month period to determine the location and chirality of 109 filaments (Yeates et al., 2007) relative to the underlying magnetic flux. In the second stage they used combined magnetic flux transport and magneto-frictional simulations (Section 3.2.3), based on actual photospheric magnetic distributions found on the Sun. Unlike previous studies, these were run for the whole six month period without ever resetting the surface field back to that found in observations or the coronal field to potential. To maintain accuracy over this time, 119 bipoles had to be emerged with properties determined from observations. Hence, the simulations were able to consider long term helicity transport across the solar surface from low to high latitudes. Yeates et al. (2008b) carried out a direct one-to-one comparison of the chirality produced by the model with the observed chirality of the filaments. An example of this can be seen in Figure 17, where Figure 17a shows the global field distribution after 108 days of evolution. The enlargement in Figure 17b shows a simulated flux rope structure with an axial field of dextral chirality. The real Hα filament that formed at this location can be seen in Figure 17c. Through studying the barbs and applying a statistical technique, the filament was determined to be of dextral chirality so the chirality formed in the simulation matches that of the filament.

Figure 17
figure 17

Example of the comparison of theory and observations performed by Yeates et al. (2008b). (a) Magnetic field distribution in the global simulation after 108 days of evolution, showing highly twisted flux ropes, weakly sheared arcades, and near potential open fields. On the central image white/black represents positive/negative flux. (b) Close up view of a dextral flux rope lying above a PIL within the simulation. (c) BBSO H.. image of the dextral filament observed at this location.

Through varying the sign and amount of helicity emerging within the bipoles, Yeates et al. (2008b) (see their Figure 5b) show that by emerging dominantly negative helicity in the northern hemisphere and positive in the southern, a 96% agreement can be found between the observations and simulations. An important feature is that the agreement is equally good for minority chirality filaments as well as for dominant chirality filaments. Therefore, the global coronal model applied was equally good in producing positive and negative helicity regions at the correct times and spatial locations to represent solar filaments. Another key feature of the simulations is that a better agreement between the observations and simulations is found the longer the simulations are run. This indicates that the Sun has a long term memory of the transport of helicity from low to high latitudes. The reason for this high agreement is described in the paper of Yeates and Mackay (2009b) where seven different mechanisms are involved in producing the observed chiralities. The results demonstrate that the combined effects of differential rotation, meridional flow, supergranular diffusion, and emerging twisted active regions are sufficient to produce the observed hemispheric pattern of filaments. While the model of Yeates et al. (2008b) obtained an excellent agreement with the observed chirality of filaments, the 6 month simulation was unable to reproduce dextral/sinistral chirality along the polar crown PILs in the northern/southern hemispheres. The cause was that the simulation was not run for long enough to allow meridional flow, which acts on a time scale of 2 yr (see Section 2.2.1), to transport the helicity produced at low latitudes up into the polar regions. A later simulation running for a full 11-year solar cycle (see Figure 4a of Yeates and Mackay, 2012) was able to obtain the correct dextral/sinistral chirality along the polar crown PILs in each hemisphere, once the duration of the simulation exceeded the meridional flow timescale. This shows that, at least in the context of solar filaments, the flux transport and magneto-frictional model adequately describes the build-up and transport of helicity from low to high latitudes on the Sun. The model has subsequently been applied to the formation of magnetic flux ropes and their subsequent loss of equilibrium (Yeates and Mackay, 2009a; Yeates et al., 2010a), a possible mechanism for the initiation of coronal mass ejections.

4.4 Plasma emission

Although simple models of the coronal plasma emission may be derived by populating PFSS magnetic field lines with hydrostatic atmospheres (Jardine et al., 2002b), self-consistent models of plasma density and temperature require full-MHD models (Section 3.4). Indeed, a key emphasis of global MHD models has been direct comparison of the steady state solutions with either white-light (Rušin et al., 2010) or multi-spectral EUV and X-Ray observations (Lionello et al., 2009; Downs et al., 2010). Early global MHD models using a simplified polytropic energy equation were able to produce a good representation of the Suns large-scale magnetic field (Mikić et al., 1999), however they were unable to reproduce realistic emission profiles in EUV and X-rays. This was attributed to the fact that they did not produce a sufficiently high density and temperature contrast between open and closed field regions. To improve the models, a more realistic energy equation has been incorporated, including the effects of thermal conduction, radiative losses, and coronal heating, along with modeling the upper chromosphere and transition region (Lionello et al., 2001, 2009; Downs et al., 2010).

Figure 18 shows a comparison of the global MHD model of Lionello et al. (2009) with observations for CR 1913 (August – September 1996). The left hand column shows three EUV pass bands (171, 195, 284 Å) along with a Yohkoh/SXT X-ray image. The other three columns show synthetic emission profiles constructed from density and temperature distributions in the global MHD simulation. Each column uses exactly the same MHD model: the only difference is the form of coronal heating applied (Hch in Equation (47)). While each solution produces roughly the same magnetic structure, the most important factor in reproducing the emissions is the coronal heating. The second column, which uses an exponential form falling off with height, gives the worst comparison. When spatially varying heating is applied (either from Schrijver et al., 2004 or the composite form of Lionello et al., 2009), the simulation captures many of the features of the observed emission. Similar results were also found in the paper of Downs et al. (2010). In an alternative comparison, Rušin et al. (2010) compared the output from the global MHD model shown in Figure 18 with white light eclipse observations from 1 August 2008. To compare the model and observations, the authors simulated the white light emission of the corona by integrating the electron density along the line-of-sight. The results of the comparison can be seen in Figure 5 of Rušin et al. (2010) where the model successfully reproduced the overall shape of the corona along with the shape and location of three helmet streamers.

Figure 18
figure 18

Comparison of observations and simulated emissions from a global steady-state MHD model. The observations are shown in the first column, while the emissions due to different coronal heating profiles are shown in the other columns. Image reproduced by permission from Figure 8 of Lionello et al. (2009), copyright by AAS.

5 Conclusions

In this review, our present day understanding of global magnetic fields on the Sun and other stars has been discussed from both observational and theoretical viewpoints. For the Sun, we now have long-running data sets of global magnetic activity. For other stars, we are just beginning to learn of the wide variety of magnetic morphologies that exist. In terms of theoretical models, recent years have seen a significant advance in global modeling techniques. Global magnetic fields may be modeled under the nonlinear force-free approximation for long periods of time, or for short periods of time using highly detailed MHD models. A key feature of these models is that they have been validated through various comparisons with observations. While our understanding has significantly increased over the last decade, there are several immediate outstanding issues. Five of these are discussed below, where the list is not exhaustive and the order in no way reflects their importance.

  1. 1.

    While we have a detailed understanding of the normal magnetic field component on the solar photosphere, the same is not true for the vector field. Between SDO and SOLIS we now have daily full disk vector magnetograms. Analysis of these will gain us an understanding of the emergence and transport of vector fields across the Sun, as well as the origin and evolution of magnetic helicity, a key component in eruptive phenomena. Additionally, vector magnetic field measurements are need not just in strong field regions, but also in weak field regions.

  2. 2.

    For theoretical models, regular observations of the vector fields mean that more constraints can be applied to input data of bipole emergences that are used to drive these models. New techniques must be developed for the incorporation of vector data into simulations. Only through this can theory and observations be truly unified.

  3. 3.

    A number of inconsistencies with the application of surface flux transport models need to be resolved. These relate to the use of additional decay terms, variations in the rate or profile of meridional flow, and the exact tilt angle applied for new emerging bipoles. Currently various combinations of all of these have been shown (when coupled with different coronal models) to give similar results. It would be useful for the various authors to resolve these issues through comparing results from their codes, driven by the variety of input data sets available. In addition, more constraints on the input data are required from observations.

  4. 4.

    While our understanding of stellar magnetic fields has greatly increased, observations are still sporadic. Long term data sets of individual stars across different spectral types are required. From this we may deduce different magnetic flux emergence and transport parameters which are critical for the next generation of dynamo models.

  5. 5.

    Global time-dependent MHD models with evolving lower boundary conditions must be developed, to provide a self-consistent model for the evolution and interaction of magnetic fields with plasma.

These issues may be addressed with the current suite of ground and space-based observatories along with developing new theoretical techniques. A key element in achieving this goal is the increasing computing power available both for real-time data reduction and for theoretical modeling. Future revisions of the article will hopefully describe answers to some of these questions. Along with answering these questions, the scope of the review in future revisions will also be increased to include additional topics such as: the connection between the Sun and Heliosphere, observations of coronal and solar prominence magnetic fields, global eruptive models and a fuller discussion of the consequences of observational limitations on the models which directly employ observations.