The upper atmosphere of the Sun is governed by the complex structure of the magnetic field. This controls the heating of the coronal plasma to over a million kelvin. Numerical experiments in the form of three-dimensional magnetohydrodynamic simulations are used to investigate the intimate interaction between magnetic field and plasma. These models allow one to synthesize the coronal emission just as it would be observed by real solar instrumentation. Large-scale models encompassing a whole active region form evolving coronal loops with properties similar to those seen in extreme ultraviolet light from the Sun, and reproduce a number of average observed quantities. This suggests that the spatial and temporal distributions of the heating as well as the energy distribution of individual heat deposition events in the model are a good representation of the real Sun. This provides evidence that the braiding of fieldlines through magneto-convective motions in the photosphere is a good concept to heat the upper atmosphere of the Sun.
Soon after the discovery in the late 1930s that the corona of the Sun is much hotter than its surface, the first proposals were made on how to heat the upper atmosphere. For a long time, models were restricted to analytical estimations or comparatively simple one-dimensional descriptions. With the Solar and Heliospheric Observatory, SOHO , a whole new world of observations became available that answered some questions, and inspired far more new problems to be solved. With the turn of the century, the first large-scale three-dimensional full magnetohydrodynamic (MHD) models of a whole active region became available, albeit scaled down in size so that the active region could fit into the computational domain .
In these and later related models, granular motions of the near-surface convection in the photosphere braid the magnetic field. At least on the scales resolved by the full three-dimensional MHD models, the driving is slower than the Alfvén speed, so the currents in the corona will change only slowly. Consequently, this process is considered direct current (DC) heating. This is in contrast to models of reduced MHD, where also faster changes can be considered, leading to alternating current (AC) heating. The price paid by the reduced MHD models is that they cannot be used to synthesize the coronal emission.
In this review, three-dimensional full MHD models will be considered that fully account for the mass, momentum and energy balance, including heating, heat conduction and radiative losses. These models allow one to synthesize observables that are directly comparable with observations. Including a full active region into the computational domain requires the computational box to have a side length of the order of 100 Mm. Currently, with reasonable computational effort, one can afford 500–1000 grid points in each direction, which then limits the grid spacing to about 100 km. This is just sufficient to resolve the driving granular motions in the photosphere and is clearly sufficient to match current extreme UV observations, which provide a spatial resolution of about 1/3′′, corresponding to about 250 km.
While the numerical models provide a synthetic corona that resembles the real solar corona in many aspects, still quite a few characteristic parameters do not match the solar case. In particular, just as in most astrophysical numerical simulations, this applies to the magnetic Reynolds number, which is because of the very small magnetic resistivity of the solar coronal plasma (see end of §2 and §6). As long as the processes converting magnetic field into heat are local, this mismatch of parameters might not be too severe. Because the mean-free-path length is still below the typical grid spacing of the models (§2, figure 1), this seems to be satisfied. By contrast, the physical processes in flares constitute a significant violation of the MHD assumptions. Thus, in this article, the emphasis will be on the active regions outside flares.
2. Coronal heating: a multi-scale multi-process phenomenon
The question of heating the corona is a problem involving multiple scales and multiple processes. While the structures we observe in the corona typically have length scales best described in units of megametre, the ultimate dissipation process might well operate on scales best described by units of metre, or even smaller. Obviously, the structures that govern the processes of coronal heating have to span many orders of magnitude. Across this huge span of length scales, many different processes will govern the relevant physics, and certainly not all can be captured by an MHD model.
On the large scales, extrapolations of the magnetic field provide a valuable tool to investigate the magnetic structures along with the electric currents that can provide the energy to heat the plasma and drive the dynamics. Such force-free models are a reasonable description in all regions where plasma-β, i.e. the ratio of thermal to magnetic pressure, is small. They have been applied to describe the corona as a whole  down to scales covering only part of an active region . When carefully examined, even potential field extrapolations of small regions, below the scale of super-granules, might be employed to study the connectivity and energetics of the (smaller scale) corona . While these force-free (and potential) models describe only the magnetic field, but not the plasma, they can be used to estimate the free energy of the magnetic field available to heat the plasma. However, different extrapolation codes do not always give the same results . In a force-free model, one can describe only a snapshot of the corona, because one assumes that the magnetic field is in an equilibrium state. By contrast, magnetofrictional models [8,9] account for the history of the changing magnetic field. With some assumptions, one can even deduce a proxy for the emissivity based on the currents along the fieldlines and produce some synthetic images of the corona . One important application often used to field extrapolations is the investigation of the topology of the magnetic field in terms of (quasi) separatrix layers, etc. This is an extended topic in its own right and here a reference to a recent textbook has to suffice .
An efficient tool to study the evolution of the magnetic field, including the interaction with the plasma, is reduced MHD. Such models are well suited to study the original fieldline braiding concept of Parker . After initial numerical experiments [13,14], MHD turbulence simulations basically confirmed that the braiding concept can provide sufficient energy to heat the corona . More recent braiding experiments will be discussed extensively in this issue . Because high resolution can be achieved in these reduced MHD experiments, one can also study the transport of Alfvén waves. These experiments show that this AC heating process supplies enough energy to the corona, too . In a more recent study, it is even claimed that the Alfvén wave mechanism is the only one to provide sufficient energy, which the fieldline braiding DC process does not . This conclusion seems premature, though, because they assume that the braiding would only happen within a single footpoint in the photosphere. By contrast, in the braiding concept, the motions of the neighbouring magnetic field concentrations in the intergranular lanes would also produce braiding, and will deliver sufficient energy to produce a hot and dense corona, as outlined in previous models (and detailed in the subsequent sections). The flux-tube tectonics model  can be considered as the basis for this concept, and indeed three-dimensional DC models produce coronae matching observed quantities in terms of energy, temperature and pressure.
By their very nature, the models discussed above have to be large scale in the sense that they can be applied only for structures where the MHD assumptions are fulfilled. One important length scale is the mean free path of the electrons, λ=vth/νee, with the thermal speed vth and the collision frequency νee. For typical values of the temperature T=106 K and the electron density n=109 cm−3 one finds λ=100 km . To illustrate this further, figure 1 shows the distribution of the electron mean-free-path lengths derived from a three-dimensional MHD model of a small coronal active region . This confirms the above rough estimate to be an upper limit. In a large part of the coronal volume, the mean free path is only of the order of 1 m, but in particular in the less dense regions its value can reach 100 km. This underlines that care has to be taken with the MHD assumptions. This also sheds light on the modelling efforts using adaptive grids to better resolve the structures that appear, e.g. in the process of current-sheet formation [21,22]. When refining the grid, one might quickly run into scales that become smaller than the mean free path, the MHD assumptions break down, and one has to consider ways to implement new physics into the models.
In particular, one-dimensional coronal models run into this difficulty with the electron mean-free-path length. The heat conductivity according to Spitzer , q∝T5/2∇T, is very sensitive to the temperature T. Known at least since the first one-dimensional coronal models including the transition region [24–26], as a consequence of the heat conductivity in the transition region the temperature gradient becomes very steep, which requires a very small grid spacing. In more recent models, the grid spacing is as small as 1 km  or even smaller when an adaptive grid is used [28,29]. As pointed out early in [30,31], other (turbulence) processes might take over well before such small scales for the temperature gradients are reached. This could then serve as a way to solve the problem of the increasing emission measure towards low temperatures that exist in one-dimensional models . Naturally, three-dimensional models such as the ones discussed below will not reach the spatial resolution of the one-dimensional models, but are limited to a vertical grid spacing of above 10 km, depending on the size of the computational domain and the number of grid points. Therefore in larger-scale three-dimensional models, the transition region cannot be resolved using Spitzer heat conduction. To cope with this problem, one can correct for the impact on the coronal pressure [33,27]. However, in the light of the above discussion comparing the electron mean-free-path length with the grid spacing, it is not clear that this is really desirable. Instead, processes other than Spitzer heat conduction will operate if the temperature gradient becomes too steep, effects not described by MHD.
The magnetic field is ultimately dissipated on very small length scales. According to classical transport theory , the magnetic resistivity in the corona should be of the order of η=1 m2 s−1. In the induction equation, the dissipative term can be of importance only if the magnetic Reynolds number Rm=Uℓ/η is of order unity or smaller. Thus for magnetic dissipation to be efficient, the length scale of the dissipation is roughly ℓ=η/U, where U is the characteristic velocity. Assuming a typical velocity of 1 km s−1 (which will be a lower limit), this yields a length scale of only ℓ=1 mm (see also §6).
So the discussion finally comes to scales where the actual velocity distribution function of ions (and electrons) has to be considered. From measurements in the solar wind, it is well known that the velocity distribution function is not a Maxwellian and that one has to fully account for plasma kinetic effects to understand its evolution in time and as a function of distance from the Sun . For a review see . In the solar wind, the plasma is basically non-collisional, and usually it is argued that in the corona the densities are high enough to ensure that collisions quickly lead to Maxwellian velocity distributions. However, the discussion of the electron mean free path already showed that assumption starts breaking down when considering spatial scales of approximately 1–10 km, which will be reached soon by numerical experiments of increasing spatial resolution. In addition, one-dimensional loop models including ion kinetics show the effects of departures from Maxwellians and the importance of non-MHD effects even in a normal coronal loop with a length of 60 Mm [36,37]. Therefore, future three-dimensional models of the corona with improved spatial resolution will have to start thinking about how to implement non-MHD effects. Models of the chromosphere already did so by considering an extended Ohm's law accounting, for example, for Hall currents .
3. Three-dimensional magnetohydrodynamic coronal models: philosophy
When describing the corona on observable scales, i.e. on scales well above 100 km and thus much larger than, for example, the dissipation lengths mentioned above, one has to and can employ MHD. For full kinetic models, in particular, when thinking of the three-dimensional structure, a description of even the smallest observable scales is not feasible.
The following discussion will focus on models that allow one to synthesize coronal emission. For this a realistic temperature and pressure of the corona are required. To ensure this, the interplay of heat input, radiative losses and heat conduction back to the Sun is essential. While this has been possible in one-dimensional models since the 1980s [26,39], it is still computationally quite expensive in three-dimensional models. Therefore, many three-dimensional MHD models use a zero-β approximation, building in that the energy density of the plasma is negligible, but sill solving for the mass and momentum balance. These models can give a realistic description of the evolution of the magnetic field in complex situations, e.g. for the formation of eruptions [40,41], but will not provide (reliable) information on the pressure, which would be essential to synthesize the coronal emission.
For the numerical experiments with realistically looking coronae, one has to solve the induction equation together with the conservation of mass, momentum and energy. In the case of the energy equation, terms for heat conduction and radiative losses in the corona are required to ensure a realistic pressure in the corona. In most of the models that will be discussed below, the heating is due to Ohmic dissipation, i.e. a term proportional to the currents squared. Some models for large-scale structures employ a parametrization of the heat input that depends on magnetic field, density, driving velocity, etc. [42,43].
The energy supply to the three-dimensional models discussed below is basically due to the braiding of fieldlines , or due to the similar process of flux-tube tectonics , which can be considered as an evolution of the braiding concept. Here typically the driving in the photosphere is at sub-Alfvénic speeds and the resistivity in these three-dimensional models is too high to properly resolve the Alfvén wave-type AC heating. Consequently, these three-dimensional full MHD models are dominated by DC-type heating (but see discussion in §6 for AC versus DC heating).
In three-dimensional models covering smaller regions, typically parts of plage or network elements, with a typical box size of some 10 Mm, one can include the near-surface convection in the photosphere and the processes in the chromosphere along with the coronal evolution [44,45]. This allows a self-consistent description from the photosphere all the way into the corona, at the price of the size of the computational domain.
Instead of the self-consistent description of the convective motions in the photosphere, one can simply prescribe the horizontal driving motions in the photosphere using a velocity driver that shares the same properties as the real solar photosphere in terms of velocity and divergence of the flow, as done by the first three-dimensional models [2,46,47]. This allows one to run the models with a coarser grid, because the convection does not have to be fully resolved, and with the same grid size one can cover a much larger computational domain, encompassing a whole active region. Furthermore, in these models one can easily use observations of the velocity and the magnetic field in the photosphere to drive the structure, dynamics and heating of the corona .
Of course, there is a price to pay for not fully self-consistently including the magneto-convection in the photosphere and the corona, including the energy transport due to radiation. While the chromosphere in these large-scale models still acts as a reservoir for mass and energy, the dynamic interaction between the chromosphere and the corona cannot be captured. In particular, spicules will not be captured in the large-scale models, and there are only very few cases where they have been found in the simulations with smaller box size and higher spatial resolution . Also many of the details of the small-scale injection and heating of cool material into the upper atmosphere as deduced from observations will not be captured [50,51]. However, even the large-scale models include the mass exchange between the chromosphere and corona, albeit at a less dynamic level than the real Sun.
(a) Inversions of observations
To get a feeling for the plasma conditions in the corona, one can use spectroscopic or imaging observations in the extreme UV (EUV, ≈150–2000 Å). Spectrographs such as SUMER , EIS  or IRIS  provide not only Doppler shifts to deduce the velocities but also ratios of emission lines that can be used to deduce temperature, density or even elemental abundances, as noted in textbooks, e.g. . To follow larger regions in time, one can also use imaging instruments such as AIA  or XRT . If the filtergrams span a large enough range of peak contribution temperatures, one can make a good guess at the thermal structure [58,59]. However, even when many spectral lines (or bands) are used, the digital elevation model (DEM) inversion might fail . Still, these inversions are often convenient, because they are relatively quick, and they always give a result (as long as they converge). However, this latter feature is dangerous, because structures along the line of sight make a reliable interpretation difficult: if there is a compact dense and a more extended less dense structure along the line of sight, the resulting density will be somewhere in the middle, not necessarily reflecting the properties of any of the structures along the line of sight . This simply reflects the fact that these inversions often are an ill-posed problem .
The inversions make implicit assumptions that are often forgotten when using an inversion code as sort of a black box, e.g. for the emission measure inversion or density analysis with tools such as Chianti . Examples for the implicit assumptions are ionization equilibrium or absence of deviations from the line formation temperature due to flows and fast temporal changes. This is not to say that the forward models discussed in the next section do not apply simplifying assumptions, but there they are more obvious and can be better controlled.
Still, inversions provide very valuable information to get estimates for, for example, density and temperature, in particular in structures that are (still) too complicated to be sensibly described by a reliable forward model. One example is the high densities found at transition region temperatures in Ellerman-bomb-type events .
(b) Compare synthesized emission to observations: forward modelling
An MHD model provides the vectors of the magnetic field and the velocity, and the scalars of temperature and density at each grid cell of the computational domain as a function of time. The same is true for one-dimensional loop models (except for the magnetic field, of course). With this information on the thermal structure and assuming that the plasma is optically thin and in ionization equilibrium, it is easy to calculate the profiles of EUV emission lines  with tools such as Chianti . Integrating along a line of sight, one can produce spectral maps very similar to what is obtained when a spectrograph would raster a region on the Sun, only with full temporal resolution. One can then treat these spectral maps in the same way as observations to produce intensity maps, Doppler maps or perform an emission measure analysis [66,65]. When synthesizing imaging observations as obtained, for example, by AIA, one can take an even simpler route and use the temperature response function for the various filters sampling different temperatures. These have been made available through SolarSoft and are documented, for example, for AIA .
These synthetic observations can be compared directly with solar observations. One can compare the structures found in the model. Examples are the intensity distribution along and across loops, the expansion of coronal structures, or the spatial structures and average value of the Doppler shifts (most of which will be discussed briefly below). If the synthetic and observed properties match, one can say at least that the forward model provides one possible explanation for the observed feature. Of course, other explanations might exist. However, when comparing complex models to complex observations, there is a fair chance that, if the synthetic and real observations match, the model is in fact a good description of the real Sun. In this context, ‘complex’ means that not just one single aspect is matching, but that many observables (in the best case any one can think of) do match. To illustrate such synthesized observations, figure 2 shows an emerging active region model covering 150×75×75 Mm3 after two sunspots appeared. Many loops have formed and evolve quickly on a time scale of 10–30 min.
If a match is found, the real challenge is to isolate the process(es) that are responsible for creating the respective features in the model. These processes might then be the ones that are responsible for creating the features on the real Sun, too.
4. Energy deposition in three-dimensional magnetohydrodynamic models
The following discussion applies to three-dimensional MHD models of active regions outside flares. An MHD model would not capture the physical processes in a flare adequately. Thus also the conclusions on the temporal and spatial distribution of the energy input apply only to, for example, active region loops.
(a) Average heat input, temperature and pressure
Already three decades ago, observations suggested that the energy flux density into an active region is of the order of 1000–10 000 W m−2 (106–107 erg cm−2 s−1) and that temperatures of 3–5 MK are reached . As a rule of thumb these numbers still hold.
In early three-dimensional models [2,3], the average energy input into the corona was of the order of 100 W m−2 and the average temperatures were around or below 1 MK. However, these are averages over the whole computational box, including a significant amount of quiet Sun area with a filling factor of the active region part being less than 50%. The energy flux into the hot active region parts in the three-dimensional models is indeed of the order of the values observationally deduced and temperatures of 3 MK and more are reached in the hot bright loops .
This mitigates the criticism by van Ballegooijen et al. , who claimed that the three-dimensional models with DC heating cannot produce sufficiently hot coronae. Indeed, despite the moderate resolution of the order of 200 km, these three-dimensional models do reach temperatures as observed and provide the energy input required to heat the corona. This is because the driving of the flux tubes rooted in the photosphere is not only done within the flux tube . Instead, also braiding of neighbouring flux tubes contributes to the energy input, just as outlined in the flux-tube tectonics model . Of course, it is correct that a three-dimensional model capturing a whole active region cannot (yet) provide the spatial resolution to fully account for the driving within individual flux tubes in the photosphere, e.g. in terms of internal twisting . However, still the total energy flux density into the corona is consistent with observations and sufficient to sustain a corona of more than 3 MK (see also §6).
In order to describe the corona properly, the pressure also has to be reproduced by the three-dimensional models. Scaling laws based on a match between observations and static loop models connect temperature, pressure, energy input and loop length. They exist for a spatial uniform heat input  and for a modification with a heat input dropping exponentially with height . Investigating individual loops (and fieldlines) shows that the three-dimensional models indeed reproduce these scaling laws, albeit with a large scatter . This applies not only to the general trend, i.e. the slope of the power law, but also to the absolute value. The huge scatter in the three-dimensional model is expected because of the large spatial and temporal fluctuations of the energy input.
(b) Heat input with height and along individual fieldlines
In the large-scale three-dimensional simulations driven by magnetic braiding (or flux-tube tectonics), the energy input per volume drops roughly exponentially in the corona. This is mainly because in the coronal part the magnetic field very roughly drops exponentially and thus in a near-force-free environment also the current and the (Ohmic) heating rate will drop exponentially . The scale height of the heating (in the coronal part) is found in various models to be in the range of 5–15 Mm. Most probably this heating scale height is related to the typical distance of the main polarities in the photosphere forming an active region. A systematic study to investigate this dependence of the heating scale height on the model set-up has not yet been performed, though.
While for most of the three-dimensional models only the horizontally averaged heating rate is shown to drop with height, this is also true for individual fieldlines . Therefore, the exponential drop of the heat input is really also applicable to one-dimensional loop models, as has often long been assumed . Another factor for the average drop of the heat input is that most fieldlines emerging from the photosphere will close to nearby opposite polarities, and will thus reach only to small heights, while only fewer fieldlines will make it up into the upper atmosphere. However, this alone does not account for the average drop of the heat input but has to be considered along with the drop of the energy input along individual fieldlines.
The fact that the heating scale height is just between the pressure scale height of the chromosphere (well below 1 Mm) and the corona (50 Mm at 1 MK) leads to the consequence that the energy input per particle will peak somewhere in the transition region from the chromosphere to the corona (figure 3). This has been found to be the case in models where the near-surface convection is treated self-consistently  as well as in cases where the horizontal velocities in the photosphere are simply prescribed by a convection-like driver .
This location of the maximum energy input per particle in the middle transition region has a major consequence on the dynamics of the transition region plasma. As pointed out in a three-dimensional model  and hinted at before by a one-dimensional work , this leads to the cooler material being pushed down and the hotter material expanding into the upper atmosphere. Consequently, this provides an explanation for the observed transition region redshifts and the coronal blueshifts [75,76]. Looking at flow profiles along fieldlines or loops in three-dimensional models, this feature can also be isolated for individual structures [3,68]. In figure 4, right panel, such a flow structure is visible. This is an important indication that the flux braiding and flux-tube tectonics underlying the three-dimensional models do not only provide a sufficient amount of energy to heat the corona, but that this process also provides the proper spatial distribution of heat input.
(c) Spatio-temporal variability of the heat input
Many one-dimensional loop studies have been conducted prescribing different distributions of the energy input in space and time. Already early loop models studied the thermal response to increasing and decreasing the heating rate in time [26,77]. This showed that asymmetric heating leads to a difference in the coronal base pressure, which drives a coronal flow. It has to be emphasized that a pressure imbalance at the loop footpoints in the photosphere alone will not drive the flow . Later studies showed that the spatial distribution of heating along a loop can produce a significantly different temperature structure and hence a different emission pattern along a loop, e.g. in X-rays . Likewise, one-dimensional loop models assuming different distributions of the heating in time could reproduce the emission from a loop as a function of time and account for the lifetime of loops being longer than the cooling time scale [80,81]. These one-dimensional studies are well suited to infer a temporal and spatial variation (along the loop) of the heat input through a comparison with actual loop observations. However, the fundamental problem here is that the spatial and temporal variation of the energy input has to be assumed ad hoc and then iterated until a match with the observation is achieved.
By contrast, three-dimensional models including the dissipation of magnetic energy will self-consistently account for the temporal and spatial variation of the energy input. In particular, these models will also describe the variability of the heating in the direction perpendicular to the fieldlines. Solving a simplified problem, it has been shown that current sheets are roughly aligned with the magnetic field and result in a heat input intermittent in space and time [13,15]. These current sheets are also seen in three-dimensional models in a more realistic atmospheric set-up (albeit at a reduced resolution) as narrow fingers of enhanced heat input; see figs. 9 and 10 of .
The heat input occurs on a wide range of time scales, from below minutes to fractions of an hour. In the three-dimensional models, there are even fieldlines that are heated on one leg in a steady way and on the other leg in an intermittent fashion (figure 4). Here ‘steady’ and ‘intermittent’ correspond to time scales longer or shorter than the cooling time, i.e. if the atmosphere has time to adjust to the changed heat input or not. This shows that the distinction between these two extreme cases probably is not meaningful, and sheds new light on the long-standing debate between the camps favouring intermittent heating  and steady heating . In reality, there will be a continuous distribution of time scales for the energy release, so that this clear distinction between steady and intermittent seems to be somewhat arbitrary, if these two extremes can be found in a three-dimensional model even on the same fieldline.
(d) Footpoint heating or coronal heating? Both!
The three-dimensional models show that the distinction between footpoint- or chromosphere- dominated heating and corona heating is of limited use. If the stressing of the magnetic field leads to strong currents in the chromosphere, in the same process also the currents in the upper parts along that fieldline will be increased, and hence the heating increases not only in the chromosphere but also in the coronal part of the respective magnetic flux tube. In a recent study, it is argued based on one-dimensional models that if the heating is very much concentrated in the chromosphere, the typical signatures of coronal loops are not reproduced . While within the framework of their model, their conclusion is correct, the underlying assumption of their model might be flawed, namely that one can limit the energy deposition to a very narrow stretch along a given fieldline (in some of their models some 500 km at the very end of a 150 Mm long fieldline). However, at least according to the existing three-dimensional MHD models this does not seem realistic. If there is a very strong magnetic disturbance in the chromosphere, that will naturally also propagate up along the fieldline and lead to an increased energy input in the upper coronal part of a loop. This can be seen by examples of the heating rate along individual fieldlines in several three-dimensional models [3,68,73].
In these three-dimensional models, the energy deposition per volume is enormous in the chromosphere; typically there is a drop of the volumetric heating rate by more than five orders of magnitude from the photosphere to the base of the corona. This might mislead one to state that the energy deposition is totally controlled in the chromosphere. The density drops even more dramatically than the energy input, so that the heating rate per particle (or unit mass) has a peak in the transition region above the chromosphere. In particular, the energy input into the coronal volume is small compared with the chromosphere, but non-negligible and still significant! Still, certainly much more energy is deposited in the chromosphere than in the corona, so in this sense the corona is the tail of the dog.
(e) Nanoflares, nanoflare storms and braiding of fieldlines
In a three-dimensional MHD model, the energy deposition is highly intermittent in space and time. Its distribution in energy roughly follows a power law , which is consistent with observations over a large range of energies and with the concept of self-organized criticality; see e.g. the review . In the coronal part of the volume, the typical energy released in a three-dimensional model is of the order of 1017 J (or 1024 erg) , which is in the same range as the nanoflares originally predicted by Parker .
Investigating small volumes of the three-dimensional model (or individual grid cells), the heat input is enhanced sometimes for only a few minutes, depositing the typical energy of a nanoflare, or many events cluster and form a longer-lasting energy deposition covering often 20 min or more ; see also figure 4. These clusters of events can be considered to be the same as ‘nanoflare storms’  that typically last several hundred to several thousand seconds. These were introduced (even though not coining this name) to match the light curves of warm loops in one-dimensional models . The three-dimensional MHD models show that not only nanoflares but also the storms are a natural consequence of the fieldline braiding and flux-tube tectonics governing the energy deposition .
The three-dimensional MHD models also modify the interpretation of the original Parker cartoon for the flux braiding (see fig. 1 of ). In that original work, the magnetic field is (roughly) constant in space, which is similar to three-dimensional experiments without gravity . Therefore, also the currents and thus the Ohmic heating will be (roughly) constant in space. By contrast, in the more realistic set-up of the three-dimensional models, the opposite polarities are in the same plane, the photosphere and thus the field expand with height, leading to a drop of the heating rate with height. In addition, in Parker's original work, gravity was not considered. This, together with the drop of the heat input per volume, leads to the maximum energy deposition per particle at the base of the corona (see §4b and figure 3). In consequence, Parker's nanoflares will not deliver the same energy throughout the atmosphere, but they will be stronger in lower heights, as found in the more realistic three-dimensional models .
(f) Heating, heat conduction and radiation
In a traditional static model for a coronal loop, there is a balance between heat input, H, heat conduction back to the Sun, ∇⋅q, and (optically thin) radiative losses, Lrad. The physics of the radiation imply that it is most effective only at temperatures below about 0.1 MK. The energy flux through heat conduction, q, is highly dependent on temperature, being most efficient at high temperatures . Consequently, in a model with heating being constant in space, the energy input locally in the transition region will play a minor role and the radiative losses there will be balanced mainly by the heat conduction back from the corona . This is part of the problem why the temperature gradient becomes so steep.
In a dynamic set-up with a concentration of the energy input towards the chromosphere (and the heating per particle having a peak in the transition region), the situation can be quite different. It has been shown that then in the energy balance the heating term H clearly dominates in the phase when the loop forms, and both ∇⋅q and Lrad play a minor role . So locally the plasma in the transition region is not heated through heat conduction from above, but through (local) Ohmic dissipation (see fig. 4a of ). Only in the later stage, after the heating has ceased, is the traditional distribution reached, where at transition region temperatures the radiative losses are mainly balanced by heat conduction from above and the local energy input plays a minor role (see fig. 4c of ).
This shows that one has to be careful in carrying over traditional concepts from (static) one-dimensional loops when investigating more complex (realistic) set-ups where the heat input is self-consistently driven by the footpoint motions in a three-dimensional model.
(g) Flux emergence
The Poynting flux into the upper atmosphere heating the plasma and driving the dynamics might be caused by several processes. Besides the driving of the fieldlines in the photosphere and chromosphere in the horizontal directions in an AC or DC fashion, emergence of magnetic flux can play a major role.
The emergence of new magnetic flux has two major consequences, namely to increase the amount of (free) magnetic energy and to alter the configuration of the magnetic field—both will have a significant impact on the coronal heating and dynamics. The consequences of the flux emergence on the coronal dynamics are many. The process can drive the formation of stable and unstable (erupting) flux ropes [91,92] or lead to recurrent eruptions if the emerging flux tube is twisted . This is a whole topic of its own covered in several reviews [94,95].
Most flux-emergence models were mainly concerned with the evolution of the magnetic field, and thus paid less (or no) attention to the thermal evolution of the plasma. However, one can couple a flux emergence model forming an active region [96,97] that covers the solar interior and the near-surface magneto-convection to a coronal model . Some results of this model have already been discussed in §4f. Here the magnetic field is basically emerging into an almost field-free region, so there is no eruptive interaction with the pre-existing magnetic field. In this case, the main effect is to provide a higher level of (free) magnetic energy, which together with the horizontal motions leads to a Poynting flux that is higher than in an evolved active region with no newly emerging flux . Future studies will have to separate the effects of adding magnetic field due to flux emergence and of horizontal motions leading to AC or DC heating in a braiding-type motion.
5. Observables from three-dimensional magnetohydrodynamic models
As detailed in §3b, one of the more interesting aspects of the three-dimensional MHD models is that they allow one to synthesize the transition region and coronal emission. A prerequisite for this is that the energy equation is treated in sufficient detail and contains not only the heating term, but also radiation and heat conduction along the magnetic field. Then the temperature and pressure (and thus the density) derived from the model will be trustworthy, which is essential to get the emission correct.
What do we learn from the emission synthesized from the MHD models? If we reproduce (at least certain aspects of) observed structures, we can at least say that the processes governing the simulation are consistent with what is happening in the Sun. Of course, one has to keep in mind that other processes might also lead to the same observables.
(a) Synthetic observations
The synthetic observables of the large-scale three-dimensional MHD models show the typical features also found in observations. Already the first three-dimensional models showed a coronal volume dominated by loops , and later models improved on this. Because, in general, the heat input is not symmetric on a fieldlline, flows along the loops will be initiated that in general look similar to observed active region siphon-flow features (cf. figure 5).
Average quantities have been synthesized to check for a ‘statistical’ match between three-dimensional MHD models and observations. These successful tests include the differential emission measure in the low transition region , the transition region redshifts [65,66,99] and coronal blueshifts , or the temporal variability of intensity and Doppler shift . However, we are still far from understanding the widths of transition region emission lines .
Comparing individual loop structures that appear in emission also shows a match in terms of the variation of the emission along the loop, and lifetime of the loop . Loops in the model match observations in terms of diameter, or cross section [102,68]. In general, loop structures found in observations, such as thin or thick loops at constant cross section, or expanding envelopes of several thin loops, are also seen in the three-dimensional MHD models . A three-dimensional MHD model that was run to match an actual observation did even reproduce the location of a set of loops in three-dimensional space, as revealed by a comparison to stereoscopic observations .
Together this provides ample evidence that the three-dimensional MHD models do a good job in reproducing the solar corona, despite all the trade-offs these models necessarily had to apply. Still, further testing needs to be done. Because of the high computational demand, the number of three-dimensional MHD models to be compared with real observations is rather limited. One can easily set up a one-dimensional loop model and run it in less than a day on a single processor (depending on the problem), so in essence one can use trial and error to tailor a one-dimensional model to match a particular observation. By contrast, the three-dimensional models run for months on highly parallel machines. So already for practical reasons fine tuning is not an option for three-dimensional models, which underlines the success these models have had so far in matching observations.
(b) Structure of coronal loops
What can the three-dimensional MHD models tell us about the formation of coronal loops? By synthesizing the coronal emission, one can locate the exact position of a given loop at any given time, and using the output from the MHD model, one can then study why that particular loop appeared and why it looks the way it does.
One of the puzzling questions on loops concerns their diameter and constant cross section. Observations show that loops in EUV and in X-rays have a roughly constant cross section [104,105], and that this diameter seems to be of the order of 1.5–2 Mm for isolated loops [106,103]. So modelling is confronted to explain not only the constant cross section, but also their finite width. The constant cross section seems to contradict the fact that the magnetic field is expanding, and the finite width raises the question as to what fundamental width across the magnetic field would determine this.
Even in a potential field extrapolation, where in general the magnetic field expands with height, the constant cross section could be an artefact of the observing geometry because the shape of a magnetic (flux) tube varies along the loop . However, then at some places one should see also strong expansion, which does not seem to be the case.
In a three-dimensional MHD model, the shape of the resulting loop in coronal emission depends not only on the magnetic field structure, but also on the thermal properties (temperature and density), which will change across the fieldlines, although no material is transported across the field [43,102]. The loop in a particular EUV band reflects only a small range in temperature, so that a very high temperature at the top of an expanding magnetic tube will effectively cut off the top part and thus lead to a loop that looks like it has a constant cross section . This scenario is also found with loops in a more dynamic situation of a model of an emerging active region . While questions remain, e.g. concerning the appearance of the model loops in X-rays, these three-dimensional MHD models underline that the combination of magnetic and thermal effects will most probably be responsible for the constant cross section.
The current models also produce EUV loops with roughly the correct width of about 2 Mm [102,68]. At this point, it is not yet clear if these models have a sufficient spatial resolution and a small enough (numeric) magnetic resistivity for solid conclusions. (The grid spacing in these simulations in the corona is 300–400 km, so the 2 Mm wide loops are barely resolved by the sixth-order finite-difference scheme employed.) Nonetheless, these models hint at the role that the change of temperature across loops  and the temporal evolution of the temperature on expanding fieldlines play .
In contrast to this process located within the corona, the magnetic field structure near the footpoints of the loop also could play a crucial role. It has been argued that the cross section can be governed by the reconnection-enabled cross-field plasma dispersal in response to the braiding motions in the photosphere . Alternatively, one could go one step further and argue that small-scale structures within magnetic concentration in the intergranular lanes could govern the loop width. Very small kilogauss flux tubes have been observed even in the quiet Sun . These would provide enough magnetic flux for a whole coronal loop. The vortex motions within the flux tubes in the photosphere [110,111] could then provide the substructure that is sort of mirrored by the loop width in the corona. Future high-resolution simulations will be needed to investigate this issue.
The intimate interaction of the three-dimensional magnetic structure with the thermal properties of the atmosphere show that one-dimensional loop models have to be treated with care. While they often give a good representation of the appearance of a coronal loop (along the magnetic field), they have to prescribe heating rates as a function of time (and space) and by design cannot account for the complex structures we see. Reconnection plays a major role in changing the connectivity and identity of fieldlines at the critical interface from the chromosphere into the corona where the plasma-β changes from above to below unity. Only in a three-dimensional model can one properly describe this place of strong interaction between plasma and magnetic field. In the end, three- and one-dimensional models are complementary: the three-dimensional models are essential to capture the complex structure and the one-dimensional model can be used for a detailed resolution of the thermal variability of the corona along selected fieldlines.
6. Summary and conclusion
In general, three-dimensional MHD models with an Ohmic-type heating induced by fieldline braiding and flux-tube tectonics produce realistic-looking coronae. The emission synthesized from these models gives a good representation of average quantities, e.g. for differential emission measure or Doppler shifts. At the same time, critical quantities of coronal loops are reproduced, e.g. width, constant cross section or intensity distribution along the loop. Still, more work is needed to compare the synthesized emission from the three-dimensional models to further aspects of observations in a similar way as one-dimensional models have been compared.
In these three-dimensional MHD models, the energy is injected into the corona by a net upward-directed Poynting flux. The magnetic field is driven by the motions in the photosphere and through this disturbances of the magnetic field are induced that heat the upper atmosphere through Ohmic dissipation. This leads to heating that is strongly concentrated in the chromosphere, but the heating in the corona is non-negligible (and important to sustain the corona), with the peak of the heating per particle located in the transition region.
Current models encompassing a whole active region do not provide the spatial resolution to capture all the fast small-scale motions that will lead to the generation of all kinds of MHD waves propagating into the corona. So consequently, in current large-scale three-dimensional models mainly DC heating is in operation. However, the energy flux provided by the braiding and flux-tube tectonics at the current resolution is already sufficient to power the corona and to sustain a temperature of several MK and a pressure compatible with observations in coronal loops.
Increasing the resolution will allow one to resolve more waves and will increase the role of AC heating. Given the estimates from reduced MHD models for wave heating, the energy flux through AC heating could be as strong as in the DC case. Assuming that the total energy input would increase by a factor of 2 by combining AC and DC, according to the scaling laws this would imply a modest increase of the temperature by some 20% and of the density by some 40%. In general, we can expect stronger small-scale activity where we already now have increased Poynting fluxes. So combining AC and DC heating, one might expect a slight overall increase in intensity and density, but not a fundamentally different appearance of the corona. Of course, this prediction will have to be checked, and increasing computer power will soon allow one to treat the full MHD problem at a resolution sufficient to resolve also the AC process. Most importantly, one should not make the clear distinction between AC and DC heating, because there is a continuous transition. Still, there are fundamental differences between AC and DC heating, in particular concerning the problem for the AC process crossing the strong gradients in the chromosphere and transition region to make it into the corona. Future models will have to show at which length and time scales most of the power is injected into the corona.
This discussion raises the question if details of the heat input matter at all. Let us assume that the (upward) Poynting flux in the photosphere sets how much energy is injected into which bundles of fieldlines, and that heat conduction is efficiently redistributing this energy. Then maybe in the end the actual small-scale process for the heating does not matter. At least, this could be the case if considering the corona on scales that we can currently observe, i.e. well above 100 km. This is way above the kinetic scales, and maybe we will lack the observable quantities to then distinguish between different microscopic models for the actual energy dissipation.
This indifference to the details of the energy dissipation might well be the reason why current three-dimensional MHD models do a surprisingly good job. Some of the parameters used in the models, foremost the magnetic resistivity, are off by a large margin when compared with classical transport theory. Even worse, one might ask the question why the actual (microscopic) dissipation of the energy should be well parametrized on the resolved (macroscopic) scales by Ohmic heating in the form of η j2. This applies to the three-dimensional full MHD models as well as to the reduced MHD models for AC heating. Again, the processes governing the (macroscopic) appearance of the corona might be quite robust and indifferent to the details. Also in weather prediction quite a few heuristic fixes are applied to the equations, often with surprising success.
Being a review this paper does not contain supporting data.
Conflict of interests
The author has no competing interests.
Sincere thanks are due to S. Bingert, F. Chen, P. Bourdin, P. Zacharias and T. Van Wettum. It has been a great pleasure working with them in the past years on the active region coronal models.
One contribution of 13 to a Theo Murphy meeting issue ‘New approaches in coronal heating’.
- Accepted March 4, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.