# A Comprehensive Radiative Magnetohydrodynamics Simulation of
Active Region Scale Flux Emergence from the Convection Zone to the Corona

Feng Chen School of Astronomy and Space Science, Nanjing University, Nanjing
210023, China<EMAIL_ADDRESS>Key Laboratory for Modern Astronomy and
Astrophysics (Nanjing University), Ministry of Education, Nanjing 210023,
China  Laboratory for Atmospheric and Space Physics, University of Colorado
Boulder, Boulder, CO 80303, USA  High Altitude Observatory, National Center
for Atmospheric Research, Boulder, CO 80307, USA  Matthias Rempel High
Altitude Observatory, National Center for Atmospheric Research, Boulder, CO
80307, USA  Yuhong Fan High Altitude Observatory, National Center for
Atmospheric Research, Boulder, CO 80307, USA Feng Chen<EMAIL_ADDRESS>
###### Abstract

We present a comprehensive radiative magnetohydrodynamic simulation of the
quiet Sun and large solar active regions. The 197 Mm wide simulation domain
spans from the uppermost convection zone to over 100 Mm in the solar corona.
Sophisticated treatments of radiative transfer and conduction transport
provide the necessary realism for synthesizing observables to compare with
remote sensing observations of the Sun. This model self-consistently
reproduces observed features of the quiet Sun, emerging and developed active
regions, and solar flares up to M class. Here, we report an overview on the
first results. The surface magnetoconvection yields an upward Poynting flux
that is dissipated in the corona and heats the plasma to over one million K.
The quiescent corona also presents ubiquitous propagating waves, jets, and
bright points with sizes down to 2 Mm. Magnetic flux bundles generated in a
solar convective dynamo emerge into the photosphere and gives rise to strong
and complex active regions with Over $10^{23}$ Mx magnetic flux. The coronal
free magnetic energy, which is about 18% of the total magnetic energy,
accumulates to about $10^{33}$ erg. The coronal magnetic field is not
forcefree, as the Lorentz force needs to balance the pressure force and
viscous stress as well as to drive magnetic field evolution. Emission measure
from $\log_{10}T=4.5$ to $\log_{10}T>7$ provides a comprehensive view on
structures and dynamics in the active region corona, such as coronal loops in
various lengths and temperatures, mass circulation by evaporation and
condensation, and eruptions from jets to large-scale mass ejections.

Radiative magnetohydrodynamics (2009), Solar magnetic fields (1503), Solar
magnetic flux emergence (2000), Solar extreme ultraviolet emission (1493),
Solar corona (1483), Solar active regions (1974)

##  1 INTRODUCTION

The magnetic field of the Sun emerges from its interior to its surface (Zwaan,
1985; Harvey & Zwaan, 1993) and further to the Sun’s atmosphere (Low, 1996).
The magnetic field plays an important role in the coupling between Sun’s
interior and atmosphere and even beyond (Solanki et al., 2006; Wedemeyer-Böhm
et al., 2009). Active regions are formed through emergence of large-scale
strong magnetic flux and give rise to solar eruptions such flares and coronal
mass ejections (van Driel-Gesztelyi & Green, 2015; Toriumi & Wang, 2019).

Magnetohydrodynamic (MHD) simulations have been the most common tool for
modeling flux emergence in particular through the upper convection zone into
the solar atmosphere (e.g., Fan, 2001; Archontis et al., 2004; Manchester et
al., 2004; Magara, 2006; Isobe et al., 2006; Arber et al., 2007; Archontis &
Török, 2008; Hood et al., 2009; Fan, 2009a; MacTaggart & Hood, 2009;
MacTaggart & Haynes, 2014; Pariat et al., 2009; Leake et al., 2013, 2014; Fang
& Fan, 2015; Takasao et al., 2015; Toriumi & Takasao, 2017). More
comprehensive reviews on numerical models of flux emergence and their
connection with observations are given by Fan (2009b, on flux emergence in the
interior), Hood et al. (2012); Cheung & Isobe (2014), and Schmieder et al.
(2014).

The aforementioned examples considered the background stratification in the
volume spanning from the upper convection zone to the solar atmosphere but
employed idealized treatments on the thermodynamics of plasma. In addition,
convective motions that pose non-negligible impacts on the evolution of
magnetic field (Nordlund et al., 2009; Stein, 2012) were also missing. A fully
self-consistent model of the interaction between magnetic field and plasma
requires a more realistic representation of plasma thermodynamics, which is
also the basis for forward synthesis of plasma emission that can be
quantitatively compared to remote-sensing observations. For this purpose, it
becomes necessary to account for the important physical processes that
dominate the evolution of thermodynamic properties of plasma in the real Sun,
such as radiative and conductive transport, energy dissipation, and ionization
effects.

A fully self-consistent and realistic model for the whole domain from the
convection zone to the corona is a very challenging task. Tremendous efforts
have been made to tackle the problem from different aspects. A very common
approach was still diving the problem into two part, which are the convection
zone to the photosphere and the photosphere to the corona.

Radiative MHD (RMHD) simulations that solve radiative transport problem for
the optically thick radiation in together with the MHD equations became
available for realistic simulations flux emergence from the convection to the
photosphere. For example, Cheung et al. (2007, 2008) simulated emergence to
small-scale flux tubes, and Tortosa-Andreu & Moreno-Insertis (2009) and
Moreno-Insertis et al. (2018) studied flux emergence on the scale of
granulations. These simulations showed a good agreement with observations.
RMHD simulations of the uppermost convection zone were extended for problems
on larger scales, such as formation of large active regions through emergence
of coherent flux tubes (Cheung et al., 2010; Rempel & Cheung, 2014) or by
imposing a uniform flux sheets (Stein et al., 2011; Stein & Nordlund, 2012).
These studies pointed the important role of convection in the evolution of
emerging magnetic flux. As attempts to study the relation between the solar
surface and the deeper convection zone, Chen et al. (2017) studied the
coupling between a convective dynamo simulation down to the bottom of the
convection zone and a simulation of the upper 30 Mm. Recently, Toriumi & Hotta
(2019) and Hotta & Iijima (2020) demonstrated simulations of flux emergence in
an unprecedented domain that covers the full depth of the convection zone.

For the domain spanning from the photosphere to the corona, the primary
problem in a realistic model is the generation and transport of energy flux
into the upper atmosphere and the dissipation of energy, which can maintain
the high temperature of the corona. Gudiksen & Nordlund (2002); Gudiksen &
Nordlund (2005a, b) presented the first realistic model of the solar corona
that considers the balance between the energy losses through optically thin
radiation, thermal conduction, and the heating through dissipation of magnetic
energy. The Poynting flux resulting from photospheric magneto-convection has
been found to be sufficient to heat the corona to over one MK. Forward
synthesis of extreme ultraviolet (EUV) emission (Peter et al., 2004, 2006)
demonstrated that coronal loops and plasma dynamics in this realistic model
are consistent with observations. This method was further extended by
incorporating magnetic and velocity field from observations. Bourdin et al.
(2013) and Warnecke & Peter (2019) showed simulations that can reproduce
specific features in observed active regions.

Combining the modeling methods that have been successfully applied on the
upper convection zone and the atmosphere in a single computational domain is
very non-trivial step. A larger vertical extension is a clear factor to
increase the computation expense. Some realistic treatments on the physical
processes, such as solving radiative transfers through multiple rays, are
intrinsically computational-expensive. The situation becomes worse, because
more stringent constraints on time-stepping are imposed by including more
physics, for example, thermal conduction in the corona.

As alternative approaches, Chen et al. (2014, 2015) tried to couple a
realistic coronal simulation to a realistic flux emergence simulation (Rempel
& Cheung, 2014), and studied the evolution of active region scale magnetic
flux into the corona. The behavior of coronal loops are found to be closely
related to the upward energy flux generated by the forming sunspots. Abbett
(2007) presented a model of the convection zone and corona in a single volume.
Using artificial source terms for the optically thick radiation in the lower
atmosphere and an empirical heating model for the corona, the model yields
thermodynamic structures and convective patterns similar to those observed on
the Sun. Fang et al. (2010, 2012) applied the same approximation to study flux
emergence and investigated the impact of convective flows on the emerging
magnetic flux tubes.

Fully self-consistent and realistic simulations arose over the past decade.
Many of these simulations further considered non-local thermal equilibrium
(NLTE) effects that are important for modeling of the chromosphere (Carlsson
et al., 2019). Wedemeyer-Böhm et al. (2012) showed that simulated vortex
motions driven by magneto-convection can reproduce and explain observations
through the layers of the atmosphere. Iijima & Yokoyama (2015), De Pontieu et
al. (2017), and Martínez-Sykora et al. (2017) studied small-scale activities
and their role in mass and energy transport though the atmosphere. Hansteen et
al. (2010, 2015) presented the coronal heating models, in which the generation
of energy flux by surface convection and the heating of corona through
dissipation is self-consistently considered in a single volume. Martínez-
Sykora et al. (2019) investigated the transport of magnetic energy through the
chromosphere. The emergence of magnetic flux tubes/sheets in small scales were
simulated by Martínez-Sykora et al. (2008), Archontis & Hansteen (2014), and
Hansteen et al. (2017). There realistic simulations are also developing with
including the effects of partial/non-equilibrium ionization and ion-neutral
interactions (Cheung & Cameron, 2012; Khomenko et al., 2018; Martínez-Sykora
et al., 2012, 2020).

Due to computational expense, the aforementioned (fully self-consistent and
realistic) simulations were limited to setups with a relatively small spatial
domain extent more representative of the quiet Sun or plage regions. A major
difficulty for extending these simulations to active region scales (for both,
the spatial domain extent and the amount of magnetic flux) is the Alfvén
velocity that can reach unrealistically high values in the atmosphere above
strong sunspots. Rempel (2017) extended the MURaM code (Vögler et al., 2005)
with coronal physics, using a computational efficient treatment that limits
the Alfvén velocity through the Boris correction and treats heat conduction in
a hyperbolic approach (see Section 2 for further detail. Using this approach,
Rempel (2017) simulated a stable active region comprised of a pair of strong
sunspots. Synthetic EUV emission shows loops-like features arching over the
corona between the two sunspots, which is inline with previous realistic
simulations that consider the domain above the photosphere. This method also
allows numerical simulations to tackle more extreme conditions during solar
eruptions with affordable computation expense. Cheung et al. (2019) presented
the first realistic simulation of a full solar flare lifecycle, in which the
simulated flare is triggered by the emergence of a small flux tube next to a
large sunspot in a pre-existing bipolar active region.

This recent progress in realistic simulations of active regions and flares and
in combination with ever increasing computational resources allow us to expand
this work by developing comprehensive active-region scale flux emergence from
the the uppermost convection zone into the corona. The new simulation
considers an unprecedented combination of a large spatial domain, a long
temporal evolution, and a high degree of realism, and provides a comprehensive
model of the quiet Sun, flare productive active regions, and solar eruptions
(flares and CMEs) in a continuous evolution. As compared with previous
simulations, quiet Sun models are extended to a much larger and higher domain;
active region models are extended to cover the whole formation stage through
flux emergence and to include the further evolution into the corona; and flare
models are extended from one single event to over 100 events with largely
different behaviors. The aim of this paper is to provide an introduction to
the new simulation and present some highlights from a preliminary analysis of
the data. Detailed analysis of the dataset for specific topics and in-depth
discussions with existing theories and observations are beyond the scope of
this paper and are more appropriate for more focused studies in the future.
The rest of the paper is organized as follows. The numerical method and setup
of the simulation are described in Section 2. Section 3 presents the general
properties of the quiet Sun and discusses the energy flux for coronal heating.
Section 4 presents the properties of the coronal magnetic field during flux
emergence into the solar atmosphere. Section 5 presents the dynamics of the
active region corona in three evolution stages. We summarize the results in
Section 6 and conclude in Section 7.

##  2 NUMERICAL SIMULATION

Figure 1: A 3D visualization of the AR run at $t{=}$27 h 21 m 38 s. The red,
green and blue arrows indicate the $x$, $y$, and $z$ directions, respectively.
The red-blue image at the bottom of the frame shows $B_{x}$ at the bottom
boundary, which presents the emerging magnetic flux bundles formed in the
solar convective dynamo simulation of Fan & Fang (2014). The gray-scale image
covering half of the horizontal domain shows $B_{z}$ in the photosphere, which
is comparable to an observed magnetogram. The short vertical $x$–$z$ cut at
the middle of the $y$ axis shows the magnetic field strength beneath the
sunspot. The green-blue colored features spreading in the volume illustrate
plasma at temperatures of 2 to 3 MK, with green for the lower and blue for the
higher temperatures. They correspond to coronal loops that arch over active
regions as seen in EUV images of the Sun. This figure is produced with VAPOR
(Li et al., 2019). An animation available in the electronic version of the
paper shows an overview of the domain with the point-of-view moving from the
top to the side of the domain.

###  2.1 Physics in the Simulation

In this study, we conduct the numerical simulation with the MURaM radiative
MHD code (Vögler et al., 2005; Rempel, 2014), which has been extended for
modeling the solar corona (Rempel, 2017). The code solves in a Cartesian
domain a set of equations for the induction of magnetic field and the
conservation laws of mass, momentum, and plasma energy, which is the sum of
the internal and kinetic energy of plasma. The energy equation also considers
sophisticated treatments on the physics that are important for thermodynamics
of the plasma in the solar atmosphere. These includes

  * •

Radiative transfer in the optically thick regime solved along 24 rays under
local thermal equilibrium approximation(Vögler et al., 2005),

  * •

Energy loss by optically thin radiation in the transition region and corona
(CHIANTI database Landi et al., 2012),

  * •

Anisotropic thermal conduction along magnetic field lines for fully ionized
plasma (Spitzer, 1962),

  * •

Tabular equation-of-state (provided by the BIFROST team Gudiksen et al.,
2011), which is a merger between OPAL (Rogers et al., 1996) and Uppsala
(Gustafsson et al., 1975) equation-of-state tables.

The chromosphere in the simulation presented here is treated in a simplified
manner by assuming local thermal equilibrium (LTE). A more sophisticated
treatment that considers non-LTE effects in the chromosphere, such as in the
BIFROST code, would be computationally too demanding for simulations on the
spatial and temporal scales presented here.

The other major challenges in simulations that contain a corona are the
stringent constraints on the time-step from large Alfvén velocities,
particularly in the low density areas above sunspots, and the large field-
aligned thermal conductivity in the corona.

The extended version of MURaM (Rempel, 2017) introduces the so-called _Boris
correction_ (Boris, 1970; Gombosi et al., 2002) to limit the Alfvén
velocities. This approach uses a semi-relativistic form of MHD equations,
where the speed of light imposes an upper limit to the Alfvén velocity. By
setting the speed of light to an artificially reduced, but still sufficiently
large value (in relation to the sound speed and advection velocities), we can
control the Alfvén velocity with negligible impacts on the dynamics of the
system.

An approach with a similar concept is also applied to the anisotropic thermal
conduction. In addition to the MHD equations, the code solves a hyperbolic
equation (i.e., non-Fickian) for the field-aligned thermal conduction flux
(Gombosi et al., 1993; Snodin et al., 2006). This approach filters out
unphysical transport velocities that are present in the parabolic system and
can limit the time-step substantially. Based on the expression for flux-
limiters such as Fisher (1986), it can be estimated that the free-streaming
limited transport velocity reaches about twice the ion speed of sound,
$c_{s}$. In order to ensure a proper treatment of MHD and free-streaming
limited heat conduction the maximum Alfvén velocity is chosen as
$\max(2c_{s},3|\mathbf{v}|)$. The implementation and choice of parameters for
the Boris correction and hyperbolic thermal conduction are discussed in Rempel
(2017).

In summary, the extended version of the MURaM code considers a compromise
between physical complexity and numerically efficient treatment that enables
comprehensive simulations in a single domain that spans from the upper-most
convection zone into the corona.

###  2.2 Simulation Setups

We conduct the simulation with two stages: quiet Sun (QS) and active region
(AR). The quiet Sun stage aims to establish a solar atmosphere that reaches a
dynamical equilibrium. It is not only a stand-alone model for studies of
coronal heating, energy and mass circulation, and small-scale dynamics in the
quiet Sun, but also the initial condition for the active region stage. The
active region stage considers a large amount magnetic flux emerging in a
spatial scale of a hundred Mm and focuses on the evolution of magnetic field
and emission structures in the active region corona, and in particular, solar
eruptions that likely occur in the complex and strong active regions in the
simulation.

For the description of the model setup and results, we first define the
notations used throughout the rest of the paper. $x$ and $y$ are the two
horizontal directions of the domain, and $z$ is the vertical direction.
$L_{x}$, $L_{y}$, and $L_{z}$ are the lengths of the domain in each direction.
$\Delta x$, $\Delta y$, $\Delta z$ are the grid spacing. In particular,
$z={0}$ is defined by the position of the photosphere, the layer in which the
average (Rosseland) optical depth ($\tau$) equals unity.

####  2.2.1 Quiet Sun Stage

The domain for the quiet Sun stage (hereafter, QS run) covers an area of
$L_{x}\times L_{y}{=}196.608^{2}$ Mm2. The total vertical extent ($L_{z}$) is
131.072 Mm, with about 18 Mm beneath the photosphere. The domain is resolved
by $1024\times 1024\times 2048$ grid points with a grid spacing of $\Delta
x\times\Delta y\times\Delta z{=}192^{2}\times 64$ km3.

The lateral boundaries of the domain are periodic for all variables. The
bottom boundary is symmetric for mass flows and all components of the magnetic
field. The top of the domain is set at about 113 Mm above the photosphere. We
use a damped boundary for the velocity, such that the velocities in the first
and second ghost cells are reduced to 50% and 25% of values in the last two
cells of the computation domain, respectively. The magnetic field in the ghost
cells is a potential field calculated based on $B_{z}$ at the top of the
domain.

The QS run is constructed by adding an isothermal corona to a surface magneto-
convection simulation with a domain that extends from the convection zone to
the upper photosphere. In order to limit transients in the initialization
process this is done in two steps. In the first step, a layer of 10 Mm is
added and the simulation is evolved for a few hours until the transition
region develops. Then the rest of the corona (about 103 Mm) is added. The
initial magnetic field is a small-scale mixed polarity field generated by a
small-scale dynamo operating in the turbulent convection zone (Vögler &
Schüssler, 2007; Rempel, 2014).

The simulation is evolved from the initial condition to a state in which the
mean temperature and density in the corona reaches an equilibrium. Then we
define this time as the beginning of the production run ($t_{QS}{=}0$), and
evolve the QS run further for about 12 solar hours.

####  2.2.2 Active Region Stage

The active region (hereafter, AR run) has a domain that covers an area of
$196.608^{2}$ Mm2 (identical to that of the QS run). The vertical extent is
122.88 Mm reaching 9.6 Mm beneath the photosphere. The domain is resolved by
$1024\times 1024\times 1920$ grid points, which leads to the same grid spacing
as for the QS run.

The boundary conditions for the lateral and top boundaries for the AR run are
identical to those of the QS run. The bottom boundary is a time-dependent
boundary coupled to a solar convective dynamo simulation by Fan & Fang (2014).
The implementation of the coupling is described in detail in Chen et al.
(2017, hereafter, Paper I). In brief, we extract a time series of magnetic and
velocity fields in a horizontal layer (a constant radius surface in spherical
coordinates) in the dynamo simulation. This time series captures the evolution
of several super-equipartition magnetic flux bundles that are self-
consistently generated in the convective dynamo. The extracted data are used
to fill the ghost cell values the AR run in this study. The benefits of
coupling are twofold. It solves the problem that global scale convective
dynamo simulations under the anelastic assumption cannot reach the solar
surface. And it provides a more realistic and physically consistent lower
boundary flux transport from the deep interior for the near-surface flux
emergence simulation.

We used a snapshot of the QS run at $t_{\rm QS}={6}$ h as the initial
condition for the AR run. When converting the QS run snapshot to the initial
condition, we cut off the lowest 8 Mm (128 grid points with $\Delta z{=}64$
km), such that the bottom boundary of the AR run is at 9.6 Mm beneath the
photosphere and the vertical domain is resolved by 1920 grid points (compared
with 2048 grid points for the QS run). In Paper I we studied the effect of
domain depth on the properties of the emerged flux. Using a more shallow
domain provides the following two advantages: (1) The time required for the
flux emergence through the convection zone is significantly shorter, which
provides a significant saving in computing in a setup with an expensive
corona. (2) A more shallow convection zone leads to a more complex active
region that offers the prospect of a more dynamic corona.

The AR run is evolved for about 48 solar hours, which is actually limited by
our computational resources. The analysis of a simulation with a similar depth
for the convection zone in Paper I showed that the peak of the unsigned
magnetic flux is reached at about 50 hours after the start of flux emergence.
Therefore this time period is sufficient to cover the emergence of magnetic
flux from the convection zone to the photosphere, the formation and growth of
active regions in the photosphere. This is the most dynamical stage during the
evolution of active regions.

##  3 RESULTS I: THE QUIET SUN CORONA

Figure 2: (a): Horizontally and temporally averaged temperature and density in
the QS run. Dotted lines around the solid curves indicate the maximum and
minimum values of the horizontally averaged temperature and density over the
temporal evolution of about 11.5 hours. The dashed line at $z{=}0$ Mm
indicates the location of the photosphere in the simulation, where the average
optical depth $\tau{=}1$, and the dashed line at $z{=}5.28$ Mm is the base of
the corona, as the corona filling factor $f_{c}>0.99$ (Equation (1)). (b):
Mean intensity of the (top-view) synthetic AIA 171, 193, and 211 channels
images shown in Figure 3. Red and blue lines are the mean temperature and
density in the corona (averaged above $z{=}5.28$ Mm), respectively, and are
plotted according to the axes in Panel (a).  Figure 3: Synthetic AIA 171, 193
and 211 channels images. The first row shows the view from the top of the
domain, which corresponds to observing the solar disk center. The middle row
is a side view along the $y$-axis. The bottom row show running difference
ratio of consecutive images in the middle row. An animation of this figure is
available in the electronic version of this paper. The animation covers an
temporal evolution of 50 min in the simulation and illustrate features like
coronal bright points, jets, campfire-like transient brightenings (Section
3.2), an eruption of magnetic flux rope in the quiet Sun (Section 3.3), and
ubiquitous upward propagating sound waves and Alfvénic disturbances (Section
3.4).  Figure 4: Examples of campfire-like small scale transient brightenings.
(a): Synthetic AIA 171 image from the top view of the whole horizontal domain.
Red boxes (10 Mm wide each) mark four examples shown in a detailed view in
Panels (b) – (e), respectively. These images are dislayed with the same color
scale (0 – 1000 DN/pixel/s) as in Figure 3.  Figure 5: Evolution of a magnetic
flux rope spontaneously formed in the quiet Sun. The blue lines are magnetic
field lines that outline the flux rope. The gray-scale images in the 3D views
show $B_{z}$ in the photosphere ($+/-100$ G). Red-frame panels present $B_{z}$
(gray-scale image, $+/-100$ G) and the horizontal velocity field (arrows) at
the (positive polarity) leg of the flux rope. Location of the red-frame panels
(about 1 Mm above the photosphere and fixed in time) is marked by the small
red frame in Panel (d). The red, green and blue arrows indicate the $x$, $y$,
and $z$ directions, respectively. The 3D visualization is produced with VAPOR
(Li et al., 2019).  Figure 6: Propagating disturbances in the quiet Sun
corona. Displayed quantities are in a $x-z$ cut through the center of domain
($y=L_{y}/2$). (a): Temporally smoothed AIA 193 intensity as defined by
Equation (2), showing evolution of the background media on time scales longer
than the 60 s averaging window. (b): Disturbances in the AIA 193 intensity as
defined by Equation (3), showing fluctuations on time scales shorter than the
60 s averaging window. The disturbance is displayed relative to the mean. The
blue dash line indicates the slit for the time-distance diagrams in Figure 7.
(c) & (d): Temporally smoothed and disturbance components of plasma pressure,
respectively. An animation available in the electronic version of the paper
shows the evolution of the background and the upward propagating disturbances
in a time period of about 70 min. The same time series is used to construct
the time-distance diagrams in Figure 7.  Figure 7: Time-distance diagrams of
the disturbances in (a) AIA 193 intensity, (b) plasma pressure, (c) $v_{x}$,
and (d) $B_{x}$. The diagrams are constructed along the slit shown in Figure
6(b) and for a time period of 70 min. Blue lines in each panel illustrate the
travel paths of sound waves (starting from 5 different time instances)
according to Equation (4). Similarly, green lines show the paths of Alfvénic
disturbances as given by Equation (5).

###  3.1 Stable and Self-Maintained Quiet Sun

After the initial relaxation stage that has been excluded from the production
run, the QS run generate a hot and tenuous corona that remains in a fully
self-maintained and extremely stable dynamical equilibrium.

We calculate the horizontally averaged temperature and density through the
whole vertical domain and the temporal averages of these vertical profiles.
Figure 2(a) shows the temporally and horizontally averaged of density and
temperature. The dotted lines next to the red/blue lines indicate the largest
deviations of the horizontally averaged quantities. Furthermore, we plot in
Figure 2(b) the temporal variation of the mean temperature and density
averaged in the whole corona. Here, we employ the following factor ($f_{c}$)
to determines if a gridpoint belongs to the corona:

$f_{c}(x,y,z)=\begin{cases}1,~{}~{}\tau<10^{-8}~{}\mathrm{and}~{}T>2\times
10^{5}\,K\\\ 0,~{}~{}\tau>10^{-8}~{}\mathrm{or}~{}~{}T<2\times
10^{5}\,K.\end{cases}$ (1)

The actual corona base is a very corrugated surface. Hereafter, for all
analysis of the QS run, we define the coronal base as the layer at $z{=}5.28$
Mm, where the horizontally averaged $f_{c}$ first reaches 0.99.

In order to illustrate the structures and activities in the simulated quiet
Sun corona, we synthesize images corresponding to observations obtained from
the Atmospheric Imaging Assembly (AIA, Lemen et al., 2012) on the Solar
Dynamics Observatory (SDO). Figure 3 shows synthetic images in the 171, 193,
and 211 channels from the top and side view. In a direct comparison of the
synthetic image with an actual AIA image, we find that the features seen in
the synthetic image are clearly sharper. This is because the numerical
simulation has a substantially higher spatial resolution (in terms of cell
size compared to AIA pixel size), and we do not account for the point spread
function of the instrument.

The simulated quiet Sun corona is much more homogeneous than the corona above
active regions (as shown later in this paper) and comprises of large-scale
structures that evolves slowly, short loops in 20–30 Mm length, and small-
scale transient brightenings. In Figure 2(b), we plot the mean intensities of
the 3 channels shown in Figure 3111The averaging of AIA intensities are not
limited to the volume of $z{>}5.28$ Mm. The AIA light curves exhibit
fluctuations on time scales of several hours and a few minutes, which
corresponds to evolution of the large-scale and small-scale structures,
respectively.

To summarize, the general thermodynamics properties quiet Sun corona is stable
and robust. This suggests a very robust balance between the heating and energy
losses in the corona, as discussed later in this paper.

###  3.2 Bright Points, Jets, and Campfires

Small-scale explosive events and high-energy jets has been observed in the
quiet Sun since decades ago (e.g., Habbal & Withbroe, 1981; Brueckner &
Bartoe, 1983; Dere et al., 1991; Aschwanden et al., 2000). The simulated quiet
Sun corona also hosts ubiquitous small-scale activities. By examining the
temporal evolution shown by the animation of Figure 3, we find numerous small-
scale brightening of any shape that appear and disappear on short time scales
of a few minutes or tens of seconds. These features resemble coronal bright
points that are commonly seen in EUV observations of the corona (see e.g.,
McIntosh & Gurman, 2005; Alipour & Safari, 2015; Madjarska, 2019, and
references therein). The bright points seen in this simulation are closely
related to the small scale mixed polarity field in the quiet Sun, which is
inline with the finding in observations (e.g., Pérez-Suárez et al., 2008; Mou
et al., 2016). This also implies that these bright points are likely a
manifestation of magnetic energy release, as suggested from both observational
and theoretical basis (e.g, Habbal, 1992; Tian et al., 2008; Zhang et al.,
2012b; Parnell & Galsgaard, 2004; Galsgaard et al., 2019).

Some of the bright points are accompanied by upsurges of plasma (i.e., jets
Raouafi et al., 2016; Shen, 2021), which are illustrated more clearly by the
synthetic AIA images from the side view in Figure 3. This is similar to the
relation between coronal jets and bright points found in observations (see
e.g., Doschek et al., 2010; Zhang et al., 2016; Mou et al., 2018; Panesar et
al., 2020, and references therein). Most of the jets in the simulation are
confined below 10 Mm in height, while some may reach a much greater height.
For example the one occurs at about $t_{\rm QS}{=}$2 h 32 m 15 s at
$(x,y){\approx}(140,45)$ (best seen in the 171 channel in the animation)
reaches about a height of 40 Mm with an average upward velocity of about 300
km/s.

Recently, ”campfires” as tiny brightening much smaller than typical bright
points are reported by Berghmans et al. (2021), and they are suggested to be
episodes of coronal heating events. Chitta et al. (2021) found that similar
features, albeit in slightly larger sizes, are also common in AIA
observations. Figure 4(a) shows at an arbitrarily chosen instance a top-view
synthetic AIA 171 image, in which many bright dots on spatial scales of a few
Mm can be easily identified, and Panels (b)–(e) highlight a few examples among
these features. They display different shapes, such as a small dot (less than
2 Mm wide), a larger dot (over 3 Mm wide), irregular dot with elongated
branches, and an elongated dot (less than 2 Mm wide but over 4 Mm long), which
are similar to the (larger events of) observed campfires.

This simulation was not designed with the particular purpose to study small
events on scales of a few Mm. Nevertheless, given the five-point-stencil of
the numerical algorithm and 192 km cell size of our simulation, the campfire
events are resolved by sufficient number of grid points. In a MURaM simulation
of a quiet Sun setup with a significantly higher resolution (48.8 km cell
size), Chen et al. (2021) found features similar to those in our simulation,
which strongly support the validity and robustness of the events we see.

###  3.3 Magnetic Flux Rope eruption

A quiet Sun flux rope eruption is made visible by the running difference ratio
of consecutive images ($(A_{i+1}-A_{i})/A_{i}$, where $A$ denotes a time
series of images) in the bottom row of Figure 3. The animation associated to
Figure 3 covers the whole evolution of this event. The selected snapshot shows
the ejection of balloon-like structure in the running difference ratio images.
It also shows a tendency of rotation and unwinding, which is clearly
noticeable in the animation. Meanwhile, a dimming region is visible the 193
channel images from both the top and side views.

We examine the magnetic field in this region and confirm the balloon-like
structure is a magnetic flux rope that is generated spontaneously by the
small-scale flux concentrations and vortex motions in the photosphere. The
evolution of the magnetic flux rope from build up to eruption is presented in
Figure 5. The structure exists long before its eruption seen in the animation
of Figure 3 and appears to be an already twisted flux bundle connecting two
magnetic flux concentrations with a separation of about 20 Mm, which is much
larger than the size of granules and more close to that of supergranules. A
persistent counterclockwise rotation (when seen from top, i.e., along $-z$) at
the positive polarity (also the more coherent) flux concentration constantly
adds twist to the flux bundle. At $t_{\rm QS}{\approx}$2 h30 m the flux bundle
becomes unstable and forms a kink shape with a dip in the middle, as shown in
Figure 5. The vortex motion of the positive polarity flux concentration
continuously stretches the dip upward, and meanwhile when seen from the side
view in Figure 3 (i.e., along the $y$-axis) the flux bundle rolls clockwise as
the animation of Figure 3 illustrates.

Vortex motions are a common pattern of the flow field in the solar
photosphere, as revealed in observations (see e.g., Brandt et al., 1988; Attie
et al., 2009; Bonet et al., 2008, 2010; Requerey et al., 2018, and references
therein) and in magneto-convection simulations (Moll et al., 2011; Shelyag et
al., 2011; Liu et al., 2019; Canivete Cuissa & Steiner, 2020). Studies also
showed that vortex flow fields may be coupled with magnetic fields and give
rise to rotating magnetic structures that extend to the atmosphere on various
spatial scales (Wedemeyer-Böhm et al., 2012; Su et al., 2012, 2014; Wedemeyer
et al., 2013; Kato & Wedemeyer, 2017; Giagkiozis et al., 2018). The flux rope
formed in this simulation is a particularly intriguing case, because it is a
coherent magnetic structure with its both legs anchoring in the photosphere.
Although each of its leg is like a magnetic swirl/tornado as shown in the
aforementioned studies, the shape and behavior of the whole structure, albeit
with a much smaller size and less magnetic flux, is more similar to flux ropes
in major solar eruptions. It leaves interesting questions of how frequent
these structures can form and if they can be detected.

###  3.4 Propagating Disturbances

The running difference ratio images in Figure 3 also reveal ubiquitous
propagating disturbances. They are most likely waves that are excited in the
lower atmosphere and propagate upward. We briefly examine the properties of
the waves in the vertical cut through the center of the domain (the $x-z$
plane at $y{=}L_{y}/2$). To better extract the wave-related disturbances from
the evolution of the background media, we decompose a quantity $u$ into a
temporal averaged part $u_{\rm mean}$ and a disturbance part $u_{\rm dist}$,
such that

$\displaystyle u_{\rm mean}(t)$
$\displaystyle=\frac{1}{a}\int_{t-a}^{t}udt^{\prime},$ (2) $\displaystyle
u_{\rm dist}(t)$ $\displaystyle=~{}u~{}-~{}u_{\rm mean},$ (3)

where $a$ represents the width of the time window for averaging and for the
analysis presented in this paper $a=60\,s$. As an example, Figure 6 displays
the mean and disturbance of the synthetic AIA 193 intensity ($I_{193}$) and
pressure ($p$). An evolution of about 70 min starting from $t_{\rm QS}{=}$1 h
51 m 26 s is shown by the associated animation, which makes the wave
disturbances much clearer for eyes.

There are in general two types of propagating disturbances. One is seen right
at the beginning of the animation as a single pulse of large dome-like
structure driven by a jet, and the disturbance amplitude is more than 10% of
the mean intensity and pressure. This is more similar to the large-scale waves
driven by major solar eruptions, which are also present in the AR run (Wang et
al., 2021). Here we are more interested in the other type, which is constantly
excited from the lower atmosphere and identified as ubiquitous weak
enhancement ($<5\%$) in the EUV intensity and pressure. We place a slit at
$x={113}$ Mm and extract the disturbances of $I_{193}$, $p$, $v_{x}$, and
$B_{x}$ along this slit for 70 min, yielding the time-distance (T-D) diagrams
shown in Figure 7. The slit is located in above a small magnetic flux
concentration, above which the magnetic field is dominated by the vertical
component. This offers a simpler background for propagation of the waves.

We can see that the T-D diagrams of $I_{193}$ and $p$ show many straight
ridges, corresponding to the upward propagating disturbances along the slit as
in the animation of Figure 6. The T-D diagrams of $v_{x}$ and $B_{x}$ display
more complex features, but the most evident ones are curved ridges indicating
a decrease of the propagation speed of the velocity and magnetic field
disturbances. To understand these features on the T-D diagrams, we calculate
the paths of disturbances traveling with the speed of sound and the Alfvén
speed by

$\displaystyle P_{cs}$ $\displaystyle=z_{0}+\int c_{s}dt,$ (4) $\displaystyle
P_{A}$ $\displaystyle=z_{0}+\int v_{A}dt.$ (5)

During this time period, the (time-dependent) largest wave speed allowed by
the Boris correction (see Section 2.1 and Rempel (2017)) is in the range of
700 to 1300 km/s, and the largest $v_{A}$ unaffected by the Boris correction
is less than 900 km. Moreover, the filling factor of the grid points with
$v_{A}>700(800)$ km/s is only about 1% (0.5%) and these points are mostly
found in the atmosphere below 5 Mm. Therefore, the propagation of Alfvén waves
in the most part of the corona, as we are interested here, is not affected by
the Boris correction. We plot in each panel of Figure 7 $P_{cs}$ and $P_{A}$
starting from $z_{0}=1$ Mm222To avoid the layers where $c_{s}$ and $v_{A}$ are
small. and 5 time instances that are about 10–15 min apart.

The disturbances of $I_{193}$ and $p$ that appear as straight ridges in the
T-D diagrams are clearly in-line with the paths of sound waves. Moreover, the
disturbances of $I_{193}$ and $p$ ($\rho_{\rm dist}$, which is not shown here,
is highly similar to $p_{\rm dist}$) coincide in time and position, which can
already be seen in Figure 6. Thus, the waves mode is compressible and is
consistent with behavior of sound waves. The constant slope of the ridges
corresponds to a speed of about 140 km/s, which just matches the speed of
sound in the corona.

Similarly, we find that most of the curved ridges in the T-D diagrams of
$v_{x}$ and $B_{x}$ disturbances show a good match with the paths of Alfvénic
disturbances. As shown by the curved travel path, the Alfvén speed is about
500 – 600 km/s in the low corona below 10 Mm, gradually decreases to about 150
km/s at about 40 Mm height, and becomes about 30 km/s in the high corona above
100 Mm. As a side note, the ratio of $c_{s}$ and $v_{A}$ along the slit also
reflect the change of plasma $\beta$ with height, as

$\beta=\frac{p}{B^{2}/(2\mu_{0})}=\frac{2}{\gamma}\frac{c_{s}}{v_{A}}\sim\frac{c_{s}}{v_{A}}.$

Waves are commonly detected in observations of the lower atmosphere (e.g.,
Jess et al., 2009; Centeno et al., 2009; Grant et al., 2018; Schmit et al.,
2020) and of the transition region and corona (see e.g., DeForest & Gurman,
1998; Krijger et al., 2001; De Pontieu et al., 2003; Tomczyk et al., 2007;
Srivastava et al., 2008, and references therein). Heggland et al. (2011) and
Kato et al. (2016) studied the propagation of self-consistently excited waves
in 2D simulations that account for the coupling between the solar atmosphere
and convection zone. As a comparison, by trading off the high spatial
resolution used in the two studies above, the simulation in the paper has a
much wider and higher domain, which allows waves to further propagate into the
corona. The analysis presented here is mostly for illustration of the presence
and basic properties of the waves. Detailed investigation is beyond the scope
of this paper.

###  3.5 Energy Fluxes for the Quiet Sun Corona

Figure 8: Balance of energy fluxes in the quiet Sun corona. (a): Conservation
of magnetic energy. The black solid line is the mean vertical Poynting flux
($S_{z}$) on the left hand side of Equation (9). The red dotted line shows the
right hand side in Equation (9), which is the sum of equivalent fluxes for the
Lorentz force work ($F_{\rm wlr}$) and resistive heating ($F_{\rm res}$) in
the corona above a certain height. The red and black lines are more or less
overlapped, demonstrating the expected balance between the input and release
of magnetic energy. The individual contribution of $F_{\rm wlr}$ and $F_{\rm
res}$ are shown by the orange dashed and dotted lines, respectively. The
orange dashed line is also very close to the red dotted line, which means that
$F_{\rm wlr}$ is much larger than $F_{\ res}$ and dominates the release of
magnetic energy. (b): Balance of energy fluxes (also as functions of time) at
the base of the quiet Sun corona ($z{=}5.28$ Mm) as given by Equation (10).
The energy deposition in the corona by the Lorentz force work ($F_{\rm wlr}$)
and resistive heating ($F_{\rm res}$) is balanced by the energy losses through
advection ($F_{\rm adv}$), optically thin radiation ($F_{\rm loss}$) and
thermal conduction ($F_{\rm cond}$). The viscous heating flux ($F_{\rm vis}$),
which does not explicitly appears in Equation (10), and measures the
conversion of $F_{\rm wlr}$ into plasma heating. Lines are noted by the names
with corresponding colors. See Section 3.5 for details.

We investigate the energy flux that provides the heating of the quiet Sun
corona from a view of the conservation law. The conservation law of magnetic
energy reads

$\frac{\partial e_{B}}{\partial
t}+\nabla\cdot\mathbf{S}=-~{}Q_{\mathrm{wlr}}-Q_{\mathrm{res}},$ (6)

where $e_{B}$ is magnetic energy density and $\mathbf{S}$ is the Poynting
flux. On the right hand side,

$Q_{\mathrm{wlr}}=\mathbf{v}\cdot(\mathbf{j}\times\mathbf{B})$

is the mechanical work done by the Lorentz force, with the current
$\mathbf{j}$ given by

$\mu_{0}\mathbf{j}=\nabla\times\mathbf{B}.$

$Q_{\mathrm{res}}$ is the dissipation that converts magnetic energy into
internal energy. In the MURaM code, $Q_{\mathrm{res}}$ is implemented by
accurately accounting for the magnetic energy that is dissipated by the
numerical scheme and by adding this amount into the plasma energy.

We integrate Equation (6) over the volume above a certain height ($V(z)$). One
the left hand side, because the quiet Sun in this simulation has reached an
equilibrium, the time derivative vanishes. The integral of the remainder
yields

$\int\limits_{V(z)}\nabla\cdot\mathbf{S}\,dV=\oint\limits_{\Sigma(z)}\mathbf{S}d\mathbf{\Sigma}=-S_{z}(z)\,L_{x}\,L_{y},$
(7)

where $\Sigma(z)$ is the surface of $V(z)$ and $S_{z}$ is the mean vertical
Poynting flux at height $z$. Because the side boundaries are periodic and the
top boundary is closed for vertical fluxes, we only need to account for the
vertical flux through the horizontal surface at height $z$.

For an energy deposition rate $Q$ on the right hand side of Equation (6), we
can define a flux $F$, such that the integral of $Q$ in this volume is equal
to the energy carried by the equivalent flux $F$ through the surface of the
volume. Given the boundary condition discussed above, only the vertical flux
through bottom of the volume matters, which leads to the relation

$\int\limits_{V(z)}Q\,dV=F(z)\,L_{x}\,L_{y}.$ (8)

With the relation in Equation (8), Equation (6) is reduced to

$S_{z}(z)=F_{\rm wlr}(z)+F_{\rm res}(z),$ (9)

where $F_{\rm wlr}$ and $F_{\rm res}$ are the equivalent energy fluxes for the
Lorentz force work and resistive heating, respectively. This equation
describes the conservation of magnetic energy when the 3D domain collapses
into the vertical direction.

We note that $S_{z}$, $F_{\rm wlr}$ and $F_{\rm res}$ are functions of height.
These vertical profiles are averaged over the 11.5 hours of evolution and
plotted in Figure 8(a). We find that $S_{z}$ that measures the input of
magnetic energy is about identical to $F_{\rm wlr}$ \+ $F_{\rm res}$ that
measures the release of magnetic energy (which is mostly converted to plasma
heating as shown later). The balance is consistent and robust over a height
range of more than 100 Mm, where the amplitudes of these fluxes change by
almost 4 orders of magnitude. The relation in Equation (9) is well satisfied
through the whole corona in the quiet Sun simulation, except near the very top
of the domain. We also plot $F_{\rm wlr}$ and $F_{\rm res}$ separately in
Figure 8(a). $F_{\rm wlr}$ clearly dominate in the release of magnetic energy
through the whole vertical domain.

We also investigate the balance of the energy fluxes that evolve the plasma
energy in the corona. When integrated over the volume V(z), the energy
equation solved by the simulation can be reformed as

$\frac{\partial E_{\rm HD}}{\partial t}=(F_{\rm adv}+F_{\rm wlr}+F_{\rm
res}+F_{\rm cond}+F_{\rm loss})L_{x}L_{y},$ (10)

where

$E_{\rm HD}=\int\limits_{V(z)}(e_{\rm int}+\frac{1}{2}\rho v^{2})dV$

is the sum of internal energy and kinetic energy, i.e., the plasma energy, in
the volume above a given height $z$.

Each term on the right hand side of Equation (10) represents a horizontally
averaged energy flux at a given height $z$. $F_{\rm adv}=(e_{\rm
int}+\frac{1}{2}\rho v^{2}+p)v_{z}$ accounts for the vertical advective flux,
$F_{\rm cond}$ is the thermal conduction flux, and $F_{\rm loss}$ is the
equivalent flux for the energy losses through optically thin radiation. We
plot the temporal evolution of the these fluxes at the corona base in Figure
8(b). In addition, we also show $F_{\rm vis}$, which is the equivalent energy
flux for the viscous heating $Q_{\rm vis}$ in the corona. The viscous heating
accounts for the dissipation of kinetic energy that contributes to plasma
heating. This part of energy is implicitly added to the internal energy,
because the simulation solves the equation of the plasma energy with a
conservative treatment 333A conservative treatment implies that a deduction in
the kinetic energy is considered as an increase in the internal energy..

Because the quiet Sun corona reaches a dynamical equilibrium and the time
derivative of $E_{\rm HD}$ vanishes, the energy deposition and losses in the
corona need to be balanced. In the analysis presented here, this is evaluated
by the energy fluxes through the corona base. Figure 8(b) demonstrate such a
balance that the energy input though $F_{\rm wlr}$ and $F_{\rm res}$ is
equalized by the losses through $F_{\rm cond}$, $F_{\rm loss}$, and $F_{\rm
adv}$. Interestingly the advective flux contributes as an energy loss. This is
due to that the total energy deposition exceeds the energy losses by thermal
conduction and radiation, and hence causes a downward flux of enthalpy plus
kinetic energy. The sum of all fluxes ($F_{\rm sum}$) fluctuates slightly but
is much smaller than the dominating energy deposition ($F_{\rm wlr}$) and loss
($F_{\rm cond}$) terms through the evolution. The temporal average of $F_{\rm
sum}$ is about $7\times 10^{3}$ erg/cm2/s, which is less than 1% of $F_{\rm
wlr}$ and hence suggests a very robust balance.

To summarize, the upward Poynting flux $S_{z}{=}9.03\times 10^{5}$ erg/cm2/s
provides the magnetic energy input and is balanced by $F_{\rm
wlr}{=}7.90\times 10^{5}$ erg/cm2/s and $F_{\rm res}{=}1.25\times 10^{5}$
erg/cm2/s, which release magnetic energy to plasma energy. $F_{\rm wlr}$ is
deposited first in the kinetic energy, and then quickly converted to plasma
heating, as $F_{\rm wlr}$ is mostly balanced with $F_{\rm vis}{=}6.75\times
10^{5}$ erg/cm2/s and their temporal evolutions are highly similar. The
magnetic energy flux into the corona contributes a direct heating of $F_{\rm
vis}+F_{\rm res}{=}8\times 10^{5}$ erg/cm2/s to coronal plasma. The residual
of $F_{\rm wlr}-F_{\rm vis}{=}1.15\times 10^{5}$ erg/cm2/s is balanced by the
work against pressure forces and eventually taken by the downward enthalpy
flux.

###  3.6 Implications on the Coronal Heating

This Simulation successfully produces a stable and self-maintained a million K
hot corona, which is a fundamental aspect of the coronal heating problem
(Klimchuk, 2006, 2015). We note that the heating of the quiet Sun in this
simulation are essentially identical to the Quiet Sun and Coronal Arcade
simulations in Rempel (2017) that has presented a thorough analysis on the
coronal heating process. Therefore, we do not repeat the analysis that has
been done in that study, and the analysis presented in this paper yields
consistent results in comparison to the simulations in Rempel (2017). We
briefly discuss the following questions related to the coronal heating.

Despite of the debate on whether the energy for the coronal heating is
provided by waves (Kudoh & Shibata, 1999; Moriyasu et al., 2004; van
Ballegooijen et al., 2011; Arregui, 2015) or fieldline braiding (Parker, 1983;
Galsgaard & Nordlund, 1996; Priest et al., 2002; Rappazzo et al., 2008;
Wilmot-Smith, 2015), it is generally accepted that in both scenarios the
energy is generated by the interaction of the magnetic field and granular
motions in the photosphere. In this simulation, we do not presume on what time
scale and spatial scale this interaction happens, instead, an upward energy
flux ($S_{z}$) of about $3\times 10^{7}$ erg/cm2/s in the photosphere, which
is similar to observed values (e.g., Yeates et al., 2014; Welsch, 2015), is
self-consistently generated by the magneto-convection. $S_{z}$ drops to about
$9\times 10^{5}$ erg/cm2/s at $z{=}5.28$ Mm, which provides sufficient energy
budget for the heating of the quiet Sun corona, as compared with the
requirement estimated from observations (Withbroe & Noyes, 1977).

As shown in Figure 8(b), the input energy fluxes are a composition of the high
and low frequency components that are supposed to represent the wave and
braiding contributions, respectively. Rempel (2017) analyzed the contribution
to the Poynting flux into the corona by dynamics on different time scales and
showed that on one hand, the Poynting flux distributes smoothly over a wide
range of time scales and there is no hard separation between the contributions
from waves and braiding. One the other hand, the energy flux into the corona
above 4.8 Mm (similar height to the coronal base defined in this study) is
mainly generated by the velocity and magnetic field on time scales larger than
10 minutes (to be compared with the Alfvén time scale of 2.5 minutes), which
clearly suggests that the energy flux is dominated by slow braiding of the
magnetic field.

In addition to direct energy deposition in the corona, heating events deep
down in the lower atmosphere can also generate upflows of heated plasma that
feed into the corona (e.g., De Pontieu et al., 2011; Samanta et al., 2019).
The jets and surges shown in Figure 3 indeed carry plasma of coronal
temperatures and their contribution can be measured as enhanced upward
enthalpy flux. However, the downward directed average enthalpy flux at the
coronal base demonstrates that upflows of heated plasma associated with these
jets and surges is not a primary source of hot plasma in the corona (at least
in this simulation). The same assessment is also applicable to coronal bright
points and campfire-like events. These heating episodes generally occur at
lower heights (see also Chen et al., 2021) and can indeed heat local plasma,
however, do not significantly impact the large scale corona above. The
dominating energy deposition is the dissipation of the Poynting flux in
extensive the corona.

The partition of the magnetic energy release and plasma heating shown in
Figure 8 indicates that the Lorentz force work is much more efficient than
resistive heating for releasing magnetic energy and consequently the viscous
dissipation contributes on average five times more than the resistive
dissipation to the heating of the coronal plasma. This is determined by the
high magnetic Prandtl number of this simulation (Rempel, 2017; Brandenburg &
Rempel, 2019), which is actually a reasonable representation of situation of
the real corona. As estimated by the ratio of the viscosity and resistivity
given by Spitzer (1962), the corona is of extremely high magnetic Prandtl
number ($10^{10}$). Rempel (2017) showed that the ratio of viscous and
resistive heatings strongly depends on the magnetic Prandtl number, however,
the sum of both are insensitive to the magnetic Prandtl number. The latter is
a reflection the energy conservation, which means that in a quiescent corona
the total heating (plus the partition into the kinetic energy, if any) is
determined by the Poynting flux that provides the available energy budget.

Furthermore, the energy conservation also implies that whether the corona can
be heated to a sufficiently high temperature is not determined by dissipation
coefficients used in a model. Even though the dissipation coefficients
(particularly the resistivity in those models based on the Ohmic heating) used
in large scale 3D simulations (Gudiksen & Nordlund, 2005a, b; Bingert & Peter,
2011; Chen et al., 2014; Hansteen et al., 2010, 2015) are supposed to be much
larger than the values in the real corona, the reason why the corona or
coronal loops in these models can be successfully heated to high temperatures
is a high upward Poynting flux that are found either over the whole coronal
base or at the footpoints of coronal loops.

The numerical scheme of the MURaM code uses numerical diffusivities that are
highly non-uniform in space. Nevertheless, Rempel (2017) estimated that the
effective viscosity and resistivity are on the order of $10^{13}$ and
$10^{11}$ cm2/s, respectively. Although the numerical resistivity is still
significantly larger than the resistivity given by Spitzer (1962) for typical
coronal conditions, the viscosity in this simulation is in average similar to
the Spitzer value of $10^{14}$ cm2/s.444The viscosity can also be anisotropic
with presence of magnetic field (Braginskii, 1965; MacTaggart et al., 2017),
which is not considered in present simulation. The dissipation of magnetic
energy is primarily by the viscous dissipation of the flows driven by the
Lorentz force on grid-spacing scale, while the role of the resistivity is to
allow for topology changes through reconnection, but direct resistive heating
is small. We suggest this is also likely the process in effect in the high
magnetic Prandtl Number corona. This scenario of energy conversion is similar
to that in reconnection models (Petschek, 1964; Longcope et al., 2009),
although the resolution used in our simulation is not sufficient to resolve
microscopic processes of magnetic reconnection.

##  4 RESULTS II: EMERGENCE OF ACTIVE REGIONS FROM THE CONVECTION ZONE TO THE
CORONA

We present an overview of the active region evolution from the aspects of
magnetic field and emission measure.

###  4.1 Formation of Active Regions in the Photosphere

Figure 9: Evolution of the photospheric and coronal magnetic field as active
regions emerge. Gray-scale images in square panels present $B_{z}$ at the
$\tau{=}1$ surface at six time instances. Color-scale panel below each gray-
scale image shows the magnetic field strength in the first 45 Mm of the
vertical domain (from 10 Mm beneath to 35 Mm above the photosphere), overlaid
by two contour lines that illustrate the asymmetry of coronal magnetic field
strength. The magnetic field strength is averaged in the $y$ direction between
the two dashed lines in the $B_{z}$ maps. The selected region for averaging
corresponds to the flux bundle that gives rise to the strong sunspot pair in
the center of the domain.

We use the QS run snapshot at $t_{QS}{=}6$ h as the initial snapshot for the
AR run after cutting out the lowermost 8 Mm of the domain and switching on the
data-driven boundary condition that imposes the velocity and magnetic field
adapted from the solar convective dynamo simulation. The super-equipartition
magnetic flux bundles generated in the solar convective dynamo start to emerge
to the photosphere. The rise of these flux bundles is predominately controlled
by convective upflows. The evolution of the flux bundles in the uppermost
convection zone and the formation of active regions in the photosphere have
been presented in detail in Paper I.

Though the depth of the convective layer is slightly deeper, the evolution of
photosphere and the convection zone in the AR run is basically the same as the
“$196\times 8$ run” in Paper I. A very brief recap of the active region
evolution is presented by the magnetograms ($B_{z}$ at the photosphere) in
Figure 9. Two major flux bundles reach the photosphere and form the largest
two sunspots pairs. The first flux bundle reaching the photosphere appears at
$t{\approx}10.5$ h in the lower half of the horizontal domain. After
$t{\approx}16$ h, the magnetic field related to the strongest flux bundle
appears in the center of the domain, and develops into a strong sunspot at
$t{\approx}22$ h. Over the whole evolution of about 50 hours, a few other flux
bundles also emerge into the photosphere and create smaller sunspot pairs, for
example, at $y{\approx}150$ Mm. We also note that the convection can destroy
the structure of a coherent flux bundle during the course of emergence. The
magnetic field removed from a strong flux bundle can be transport to the
surface by convective upflows. This process gives rise to many magnetic
concentrations in spatial scales of about 10 Mm. The small sunspot pairs and
magnetic concentrations add a significant complexity to the magnetic field in
the simulation, and play an crucial role in producing solar eruptions in the
active regions.

The small panel below each magnetogram in Figure 9 display the development of
the strong sunspot pair in the center of the domain. The quantity shown here
is the magnetic field strength averaged between $y{\approx 80}$ and 115 Mm,
i.e., between the two dotted white lines drawn on the magnetograms. The 45 Mm
vertical field-of-view of the small panels spans from 10 Mm beneath to about
35 Mm above the photosphere. The higher domain is not shown here. Two contours
are drawn in each panel to illustrate the asymmetry in the magnetic field
strength in the corona above the asymmetric photospheric sunspots.

In the first 10 hours when this flux bundle has not yet entered the domain
through the bottom boundary, the convection zone is filled with the small-
scale magnetic flux that is maintained by a small-scale dynamo ( also known as
local dynamo). When the flux bundle emerges close to the photosphere and give
rise to a few small and fragmented flux concentrations, it already leads to a
significant increase in the mean coronal magnetic field strength from about 10
G at $t{\approx}10.5$ h to over 60 G at $t{\approx}16$ h. The flux
concentration of the leading polarity located at $(x,y){\approx}(120,95)$
becomes a strong and coherent sunspot at $t{\approx}22$ h. At this moment, we
can see beneath the photospheric sunspot a strong vertical magnetic flux tube,
which significantly expands in the corona above the sunspot and creates a
dome-like area where the mean coronal field strength is about 300 G. After the
sunspot pair become fully developed, the mean field strength in the lowermost
35 Mm above the photosphere can reach about 500 G.

The sunspot pairs formed in this simulation show clear asymmetries that are
similar to those found in observed sunspots, in the sense of a stronger and
more coherent leading spot and weaker and less coherent trailing spot. These
properties are inherited from the flux bundles generated in the solar
convective dynamo, which were demonstrated in Paper I. It is not surprising
that the asymmetry extends to, if not the whole vertical domain, at least the
lower corona. The coronal field strength on the leading spot side is clearly
stronger than that on the trailing spot side. The difference is illustrate in
the small panels by the contour lines of coronal field strength that are
inclined between the two sunspots. At the three time instances between
$t{\approx}22$ and 32 h when the asymmetry is the most evident, the mean
magnetic field strength (between $z{=}0$ and 35 Mm) above the leading sunspot
(about 230, 450, 600 G) is about two times larger than that above the trailing
sunspot (about 130, 200, 300 G). Later when the trailing sunspot becomes
stronger, the asymmetry becomes less significant but is still present.

We plot in Figure 10(a) the total unsigned magnetic flux at the photosphere,
as well as the flux emergence rate, which is given by the time derivative of
the magnetic flux. The total unsigned magnetic flux for the 1972 Mm2 quiet Sun
is $2\times 10^{22}$ Mx. After 10 hours, as emerging flux bundles arrive in
the photosphere, the magnetic flux begin to increase at relative stable rate
of about $4\times 10^{21}$ Mx/hour. The flux continues to grow to about
$1.3\times 10^{23}$ Mx at end of the AR run. In the later half of the
evolution, the mean emergence rate also decreases to about $2\times 10^{21}$
Mx/hour. The evolution from $t{=}8$ to 48 hours leads to a mean flux emergence
rate of $2.5\times 10^{21}$ Mx/hour.

Considering the domain as a whole active region, its magnetic flux is
comparable to the largest active regions observed on the Sun, for example NOAA
active region 12192 that contains $1.6\times 10^{23}$ Mx unsigned flux (Sun et
al., 2015). This big simulated active region is accompanied by a faster flux
emergence than smaller active regions (Otsuji et al., 2011; Toriumi et al.,
2014; Norton et al., 2017)555When compared with results shown in signed flux
and signed emergence rate, the unsigned flux and unsigned emergence rate shown
in Figure 10(a) is divided by 2, as the flux in the photosphere is balanced in
this simulation.. Sun & Norton (2017) reported the fastest flux emergence ever
observed in an active region with about half of the peak flux of this
simulation. The observed mean (peak) emergence rate is about $4.9\times
10^{20}$ ($1.1\times 10^{21}$) Mx/hour, which is about $1/5$ ($1/4$) of
emergence rate in the simulation. The numerical experiments in Paper I shows
that when the bottom boundary is placed deeper (e.g., 18 Mm and 30 Mm beneath
the photosphere), it takes longer for flux emergence and less unsigned flux
emerges to the photosphere, yielding a reduced flux emergence rate (about
$1/3$ of the current simulation) that is more close to observations.

###  4.2 Evolution and Distribution of Magnetic Energy

Figure 10: Temporal evolution of magnetic flux emerged in the photosphere and
the magnetic energy contained in domain above the photoshere. (a): The solid
line shows the total unsigned vertical magnetic flux. The dotted line is flux
emergence rate, which is evaluated by the time derivative of the magnetic
flux. (b): The black and red lines plot the total magnetic energy ($E_{\rm
tot}$) and free magnetic energy ($E_{\rm free}$), respectively. In particular,
the solid red line is for $E_{\rm free}$ in the domain above the photosphere,
and the dotted red line only accounts for $E_{\rm free}$ in the domain above
$z{=}2$ Mm, which corresponds to the low $\beta$ corona. Note that the range
of axis for $E_{\rm tot}$ (left) is five times larger than that for $E_{\rm
free}$ (right). A zoomed view for the evolution of the free magnetic energy in
a time period of two hours is presented in the insert. This highlights the
decay of free magnetic energy during a solar eruption. The solid and dotted
red lines have the same meaning as in Panel (b). (c): Rate of change of the
total and free magnetic energies from panel (b). Note again that the range of
the axis on the left is three times larges than that on the right.  Figure 11:
Distribution of the free magnetic energy with height ($z{>}0$ Mm). (a) Time-
height diagram showing the evolution of the horizontally averaged free
magnetic energy density ($e_{\rm free}$). The contour lines of $f_{\rm cum}$
indicate the percentage of the free energy contained in the volume below a
certain height, as given by Equation (14). The dashed lines mark three
evolution stages, over which temporal averages of $e_{\rm free}$ are
calculated, respectively. (b): Horizontally and temporally averaged magnetic
energy density in the quiet Sun (QS). The temporal average is done for a
period of 5.6 hours starting from $t{=}0$ h, when the magnetic flux of
emerging active regions has not reached the photosphere. The solid black line
plots the total magnetic energy density ($e_{\rm tot}$) and the solid red line
is $e_{\rm free}$. The dotted black line is the ratio between the red and
black lines, i.e., the percentage of the free energy among the total magnetic
energy. (c): The same analysis as in Panel (b) for the early stage of active
region evolution (AR early). The temporal average is done for a period of 11.4
hours starting from $t{\approx}10$ h. The notation of lines are identical to
those in Panel (b). (d): The same analysis as in Panel (b) for the late stage
of active region evolution (AR late). The temporal average is done for a
period of about 14.4 hours starting from $t{\approx}32$ h. The notation of
lines are identical to those in Panel (b).

In this section, we study the temporal evolution and spatial distribution of
the magnetic energy in the atmosphere.

The total magnetic energy in the atmosphere, $E_{\mathrm{tot}}$ is given by

$E_{\rm tot}=\int\limits_{V(z{>}0)}\frac{\mathbf{B}^{2}}{2\mu_{0}}\,dV.$ (11)

Then we obtain the free magnetic energy, $E_{\mathrm{free}}$, by

$E_{\rm free}=\int\limits_{V(z{>}0)}\frac{\mathbf{B}^{2}-\mathbf{B}^{2}_{\rm
pot}}{2\mu_{0}}\,dV,$ (12)

where $\mathbf{B}_{\rm pot}$ the potential magnetic field that is computed
from $B_{z}$ at $z{=}0$ Mm and assuming $B_{z}$ vanishes at infinity ,while
horizontal boundaries are periodic as in the simulation. The temporal
evolution of the magnetic energy is shown in Figure 10(b). In particular, we
plot two lines for the free magnetic energy. One line (solid) is for the
domain above the photosphere, as defined in Equation (12). The other line
(dotted) only accounts for the domain of $z{>}\,2$Mm when evaluating the
integral in Equation (12). This part of the domain belongs to the low plasma
$\beta$ regime as shown later in this paper in Section 4.3. The change rates
of the total and free magnetic energy are evaluated by their time derivatives
and are plotted in Figure 10(c) with the same line style as in Panel (b). We
also note that the curves of the free energy / change rate are scaled
according to the axis on the right of Panel (b)/(c), which has a smaller range
compared to the axis on the left.

In the volume of 197${}^{2}\times{113}$ Mm3, the quiet Sun atmosphere contains
$8.5\times 10^{30}$ erg total magnetic energy. During the first few hours when
the emerging flux bundles are mostly in the convection zone, $E_{\rm tot}$
increases very slowly to $3.5\times 10^{31}$ erg at about $t{=}10$ h. As
active region magnetic field reaches into the photosphere and further expands
into the corona, $E_{\rm tot}$ starts to increase rapidly and grows to
$7.8\times 10^{32}$ erg at $t{=}20$ h, with an average growth rate of
$7.4\times 10^{31}\,$erg/hour. Thereafter, $E_{\rm tot}$ grows at the
relatively stable rate of $1.38\times 10^{32}\,$erg/hour and reaches about
$4.7\times 10^{33}$ erg by the end of the simulation.

The curve of $E_{\rm free}$ follows the general trend of $E_{\rm tot}$, and
the overall temporal evolution of both the total and free magnetic energy is
clearly correlated with the photospheric magnetic flux (linear correlation
coefficient $>0.95$). This is a natural consequence of the fact that in this
simulation the pre-existing magnetic field in the quiet Sun is weak and the
magnetic energy in the atmosphere is dominated by the newly-emerged strong
active region magnetic field.

The similarity between the curves of $E_{\rm free}$ and $E_{\rm tot}$ as shown
in Figure 10(b) implies that the percentage of the free energy, i.e., the
ratio of $E_{\rm free}$ and $E_{\rm tot}$, is a relatively constant value
during the evolution. In the time period between $t{=}10$ and 48 h, $E_{\rm
free}/E_{\rm tot}$ varies between 15% and 20% with a mean value of 17.9%,
which is equivalent to $E_{\rm tot}/E_{\rm pot}=1.22$, where $E_{\rm pot}$ is
the energy of the potential field that has been used to evaluate $E_{\rm
free}$ in Equation (12). In investigations of real solar active regions, it
seems that there is no general agreement on the percentage of free energy in
the total magnetic energy. In addition to the fact that this ratio is likely
to be highly case-dependent, the estimate of the magnetic energy in a
particular active region relies on force-free modeling of the coronal magnetic
field (Wiegelmann & Sakurai, 2012; Wiegelmann et al., 2014). The magnetic
energy estimated by different extrapolation methods are noticeably different,
and values of $E_{\rm tot}/E_{\rm pot}$ are usually between 1.0 and 1.3 but
may also be found below unity or over 1.5 (Schrijver et al., 2008; De Rosa et
al., 2009). Aschwanden et al. (2014) obtained $E_{\rm tot}/E_{\rm
pot}{=}1.106$ from 172 samples, which means $E_{\rm free}/E_{\rm
tot}{\approx}9.6\%$.

During certain times $E_{\rm free}$ shows a sudden decrease. One of such cases
is highlighted by the insert in Figure 10(b), where the free energy drops by
more than $2\times 10^{31}$ erg in about 20 minutes with peak energy decay
rate exceeding $10^{32}$ erg/hour. These sudden decreases are caused by
eruptions (flares and coronal mass ejections) occurring in the simulation. The
events already visible in the evolution of the total free magnetic energy
integrated over the whole domain are very energetic with a decrease of
magnetic energy on the order of $10^{31}$ erg. Many smaller events are present
as well, they are hardly discerned in Figure 10(b), but are clearly
illustrated by the change rate of the free magnetic energy shown in Figure
10(c). The behavior of gradual built up and sudden releases of magnetic (free)
energy is highly consistent with observations of emergence of active regions
leading to solar eruptions, for example, NOAA active regions 11072 (Liu &
Schuck, 2012), 11158 (Sun et al., 2012; Liu & Schuck, 2012; Tziotziou et al.,
2013), 11283 (Jiang et al., 2014, 2016), and 11817 (Liu et al., 2016).

In addition to the temporal evolution of the magnetic energy in the whole
domain, it is also interesting to investigate the distribution of magnetic
energy, in particular the free energy, within the domain. For this purpose, we
evaluate the horizontally averaged free energy density,
$e_{\mathrm{free}}(z)$, by

$e_{\rm
free}(z)=\frac{1}{L_{x}L_{y}}\int\limits_{0}^{L_{y}}\int\limits_{0}^{L_{x}}\frac{\mathbf{B}^{2}-\mathbf{B}_{\rm
pot}^{2}}{2\mu_{0}}\,dx\,dy.$ (13)

The horizontally averaged total magnetic energy, $e_{\rm tot}(z)$, can be
estimated in the same manner. Furthermore, we evaluate the following factor,
$f_{\rm cum}$, by

$f_{\rm cum}(z)=\frac{\int_{0}^{z}e_{\rm
free}(z^{\prime})dz^{\prime}}{\int_{0}^{L_{z}}e_{\rm
free}(z^{\prime})dz^{\prime}}\times 100\%,$ (14)

which describes how much of the free energy are contained in the volume below
a certain height, and obviously it equals to 100% at $z=L_{z}$. Figure 11(a)
shows a time-height diagram that illustrate the distribution of $e_{\rm free}$
through the vertical domain and evolution through the whole course of the
simulation, with the contours lines indicating $f_{\rm cum}$ levels of 50%,
80% and 95%.

A quantitative comparison between the distribution of free magnetic energy
above the quiet Sun, and active regions (in early and later evolution stages)
is shown in Figure 11(b), (c), & (d). For the quiet Sun stage, we calculate
the temporal average of $e_{\rm free}$ and $e_{\rm free}$ over a time period
of 5.6 hours starting from $t{=}0$ h, and for the active region, the temporal
average is done over a time period of 11.4 (14.4) hours starting from
$t{\approx}10(32)$ h for the early (late) evolution stage, as marked by the
dashed lines in Figure 11(a).

While in the later stage of evolution, the increase of $e_{\rm free}$ occurs
more in the lower half of the vertical domain.

Before the emerging magnetic flux reaches the photosphere, the atmosphere is
still dominated by the small scale magnetic field in the quiet Sun. The
contour lines in Figure 10 suggests that 50% of the free energy is
concentrated in the volume below $z{=}1.2$ Mm and the lowermost 20 Mm above
the photosphere already contains most ($\sim 95\%$) of the free energy. Figure
11(b) shows that the magnetic energy above the quiet Sun drops quickly with
height in the first few Mm, then both $e_{\rm free}$ and $e_{\rm tot}$
continue to decrease in a flatter slope (on logarithmic scale), while $e_{\rm
tot}$ seems to decay slightly faster than $e_{\rm free}$. Above $z{\approx}50$
Mm the two quantities are almost identical to each other, which is also
indicated by their ratio (blue dotted line) that the upper half of the domain
is dominated by free energy (close to 100%). This is because the (upper) quiet
Sun corona is in high $\beta$ regime as mentioned in Section 3.4, and the flow
field gives rise to a turbulent and highly non-potential field. If we
integrate for the energy in all height (i.e., $E_{\rm free}$ and $E_{\rm tot}$
shown in Figure 10), the free energy is about 40% of the total magnetic energy
in the quiet Sun atmosphere.

The “AR early” stage represents a dynamic transition stage from pure quiet Sun
to developed active region through fast flux emergence. Figure 11(c) shows
that, as the emerged active region magnetic field quickly fills the domain,
the magnetic energy in the upper half of the corona is increased by more than
100 times from the quiet Sun values, which makes $e_{\rm free}$ above 50 Mm
about comparable to that above well developed active regions in the late stage
of evolution. This stage exhibits the most extensive distribution of the free
energy with height, as illustrated by the rather flat slope of the $e_{\rm
free}$ ($e_{\rm tot}$ as well) above $z{=}40$ Mm. Meanwhile in the lower half
of the domain, the increase in $e_{\rm tot}$ is clearly more than the increase
in $e_{\rm free}$, which is also demonstrated by the large decrease of their
ratio (e.g., from 60% to 20% at $z{=}20$ Mm and from 30% to less than 10% near
the photosphere). This suggest a significant growth in the potential
component, as strong and stable magnetic flux concentrations form in the
photosphere.

During the “AR late” stage flux emergence is still on going and coronal
magnetic energy continues to grow. Figure 11(d) present that both $e_{\rm
free}$ and $e_{\rm tot}$ in the lower half of the domain increase
significantly by an order of magnitude, and now their vertical profiles decay
exponentially (almost constant slope on logarithmic scale) for the most part
of the atmosphere. The contour lines in Panel (a) shows that about 50% (95%)
of the free energy is contained in the first 10 (40) Mm, which is a much
greater height range compared with the situation of the quiet Sun shown in
Panel (b). The percentage of free energy in the total magnetic energy
increases with height in the most of the vertical domain from about 5% near
the photosphere to about 80% above 100 Mm. It is very interesting to note the
even above strong active regions, the magnetic field in the high corona could
be highly non-potential. When integrated over the whole atmosphere, the free
energy is about 18% of the total magnetic energy.

###  4.3 Force-freeness of the Magnetic Field

Figure 12: Transition from the strongly non-forcefree magnetic field and
high-$\beta$ plasma in the interior to the (close to) forcefree magnetic field
and low-$\beta$ plasma in the corona. (a) The probability density function
(PDF) of the angle between current and magnetic field vectors. The angle is
derived from $r_{f}$ in Equation (15). The PDF is computed for an interval of
$1^{\circ}$. (b) A zoomed view of the same quantity as in Panel (a) in the
domain between $z{=}-5$ and 5 Mm. (c) The PDF of plasma $\beta$, which is
computed for an interval of $\log_{10}\beta=0.1$. Areas filled by the white
background indicate an absent of data point in that range. (d) A zoomed view
of the same quantity in Panel (c) in the domain between $z{=}-5$ and 5 Mm. The
white line in each panel indicates the most probable value of the quantity
considered in that panel as a function of height.

It has been commonly accepted that the magnetic energy dominates over the
plasma energy in the corona and the plasma $\beta$ that measures the ratio of
gas pressure to magnetic pressure is much smaller than unity. In a static
low-$\beta$ corona, if the gas pressure could be ignored, the Lorentz force of
magnetic field must vanish. This leads to the idea of a force-free field,
which is a widely exploited approach to model the corona magnetic field
(Wiegelmann & Sakurai, 2012; Wiegelmann et al., 2014).

In reality, the non-vanishing gas pressure in the corona has to be balanced by
the Lorentz force, which leads to deviations from the force-free state. In
addition, the corona is not static. The corona is heated by dynamical
processes that are a combination of slow braiding of field lines (as in DC
heating) and waves (as in AC heating). As shown in Section 3, the magnetic
energy released through the Lorentz force work is balanced by viscous heating,
suggesting a balance between the Lorentz force and viscous stress. Finally, in
a dynamical corona the evolution of the magnetic field is driven by the
Lorentz force, which can be significant. This raises the question of how
force-free the coronal magnetic field is and how reliable the modeling of
coronal magnetic field is (Peter et al., 2015). Moreover, the first few Mm
above the photosphere are a region where magnetic field and plasma
thermodynamics drastically change, which certainly leads to a great change in
$\beta$ and force-freeness of magnetic field.

The simulation presented in this study is an exceptional test case for
evaluating magnetic forces with non-vanishing pressure gradient force of the
plasma. We summarize a few reasons. First of all, its domain covers the two
extreme regimes, so to say the plasma-dominating interior and magnetic-
dominating corona. Then, by including realistic physical processes in the
energy transport, the simulation produces realistic thermodynamics quantities,
which can provide a more sensible representation of the plasma forces in the
atmosphere. Moreover, the plasma and magnetic field are evolved fully
consistently in an MHD system. Last but not least, the simulation present a
case of strong active regions, which is a more intriguing case in real
observations.

To assess the force-freeness of magnetic field, we use the follow factor

$r_{f}=\frac{\lVert\mathbf{j}\times\mathbf{B}\rVert}{\lVert\mathbf{j}\rVert\cdot\lVert\mathbf{B}\rVert},$
(15)

which is the ratio between the current perpendicular current and the full
current. Obviously, $r_{f}{=}0$ represents a completely force-free field,
where currents are completely field-aligned, and the larger $r_{f}$ is, the
more non-forcefree the magnetic field becomes. The (sharp) angle between
$\mathbf{j}$ and $\mathbf{B}$ is given by $\arcsin(r_{f})$.

In each horizontal layer through the whole vertical domain from $z{=}-9.6$ Mm
to $z{=}113$ Mm, we calculate the histrograms of $\log_{10}\beta$ and
$\arcsin(r_{f})$. The former is done for a range between $10^{-6}$ and
$10^{6}$ with an interval of $\Delta\,\log_{10}\beta{=}0.1$, and the latter
uses an interval of $1^{\circ}$ for angles ranging from 0 to $90^{\circ}$.
Figure 12 displays the probability density functions (PDFs) of
$\log_{10}\beta$ and $\arcsin(r_{f})$. In particular, we also plot detailed
views of the domain between $z{=}-5$ Mm and $z{=}5$ Mm, where $\beta$ and
$r_{f}$ change the most significantly. We note that the histogram at each
height is normalized to unity, respectively.

We first look at the gas-dominating domain below the photosphere. The PDF of
$\arcsin(r_{f})$ shows (most probable) mean angles that are larger than 80∘,
which indicates a highly non-forcefree condition. The PDF has a very wide wing
towards smaller angles, but vanishes as the angle approaches to zero. Panel
(c) shows that almost all grid points below $z{=}0$ Mm fall into the
high-$\beta$ side and the (most probable) mean $\beta$ is over $10^{4}$ for
most of the convective layers. Only in the layer very close to the photosphere
we can see the PDF of $\beta$ extending to the low $\beta$ side, which is
probably due to the strong magnetic flux concentrations, where the magnetic
pressure prevails the gas pressure.

The observation-relevant and the most interesting region is the domain above
the photoshere, where both $r_{f}$ and $\beta$ exhibit the opposite trend
compared with the convective zone. The most noticeable difference in the PDF
shown in Figure 12(a) is that the mean angle jumps to about $5^{\circ}$ and
remains in similar values until $z{\approx}90$ Mm. The shape of the
distribution becomes much narrower, in the sense that the majority of data
points gathering between 0 and $15^{\circ}$. As showed in the zoomed view in
Panel (b), the transition of $r_{f}$ between the two extremely different
regimes occurs in a very thin layer between $z{\approx}0.3$ Mm and
$z{\approx}1$ Mm, where $\arcsin(r_{f})$ values are roughly evenly distributed
between $10^{\circ}$ and $90^{\circ}$, and the mean value changes from over
$80^{\circ}$ to $20^{\circ}$ across this transition layer.

The PDF of $\beta$ also changes drastically from the high $\beta$ side to the
low $\beta$ side. Figure 12(c) shows that the mean $\beta$ quickly drops from
around unity to $10^{-3}$ in the first few Mm above the photosphere, and the
lowest possible $\beta$ obtained in this snapshot (almost $10^{-6}$) is also
found in this region. Above this height, the mean $\beta$ gradually increases
and remains below unity for the rest of the corona. A detailed inspection on
distribution of $\beta$ around the photosphere is presented in Panel (d). The
transition of $\beta$ across unity is found in the layers between
$z{\approx}-1$ Mm and $z{\approx}1.5$ Mm, where the PDF completely shifts from
the right-hand-side of $\beta{=}1$ to the left-hand-side with the mean value
jumping from $10^{3}$ to almost $10^{-3}$.

Figure 13: Horizontally averaged $r_{f}$ as a function of height. The black
line shows the simple arithmetic mean at each height. The blue line plots the
current-weighted mean that uses the amplitude of current
($\lVert\mathbf{j}\rVert$) as weights. The red line is the most probable
values of $r_{f}$ at each height as given by the PDF, and it is equivalent to
the white line shown in Figure 12(a), whereas $r_{f}$ is shown in this plot,
instead of $\arcsin(r_{f})$ in Figure 12(a).

At last, we compare in Figure 13 the most probable value of $r_{f}$ that is
given by the PDF displayed in Figure 12 and the current-weighted mean of
$r_{f}$, which is often used to evaluate the forcefreeness of (extrapolated)
coronal magnetic field (see De Rosa et al., 2009), and we also plot the simple
arithmetic mean of $r_{f}$ as a reference. As we have seen in Figure 12 that
the PDF of $r_{f}$ is highly asymmetric with respect to its peak, the
arithmetic and current weighted means are clearly deviated from the most
probable value. When we compare the simple arithmetic mean and the current-
weighted mean, they are very similar below $z{=}0$ Mm, which suggests that the
current is about evenly distributed in all $r_{f}$. However, above $z{=}0$ and
after a thin transition layer, the arithmetic and current-weighted means
become substantially different. The current-weighted mean is much smaller than
the arithmetic mean (by a factor of 0.5 – 0.8). This means that stronger
currents are more aligned (i.e., making a smaller angle) with the magnetic
field. The current-weighted mean of $r_{f}$ is about 0.15 ($8.6^{\circ}$) in
the lower corona below $z{=}30$ Mm, and gradually increases to around 0.2
($11^{\circ}$) in the middle part and 0.4 ($23^{\circ}$) near $z{=}100$ Mm.
The result given by this simulation seems to be similar to the those given by
extrapolations (De Rosa et al., 2009), although these extrapolations assumed
the magnetic field to be forcefree. However, it is important to note that the
properties of $r_{f}$ obtained in this simulation is a physical result that
magnetic forces need to balance the plasma pressure force and the viscous
stress that dissipates the magnetic energy released through Lorentz force work
and heat the coronal plasma.

As a side note, the layer near the top boundary (above $z{\approx}105$ Mm)
differs from the rest of the coronal volume. This is because we implemented a
closed top boundary for mass flux. Gas pressure from plasma accumulated near
the boundary layer is roughly balanced with the magnetic pressure. This only
occurs near the top boundary and is not representative for the majority of the
corona. We also conducted a simple experiment, which proves that if we open
the top boundary for outgoing mass flux, the trapped plasma can easily
propagate out of the domain.

To summarize, the distribution of plasma $\beta$ in the simulation
successfully reproduces the distribution derived from observations of the real
solar atmosphere (Gary, 2001). The plasma $\beta$ in this simulation is
determined by the magnetic pressure and gas pressure that is established
through fully self-consistent evolution (flux emergence from the convection
zone to the corona) and self-consistently maintained coronal energy balance
(heating by dissipation of energy flux into the corona, cooling through
radiation, and redistribution of thermal energy by conduction). Similar to the
thermodynamics of the quiet Sun corona presented Section 3, this is a clear
hint that the simulated active region resembles the properties of the real
solar corona and reproduce a correct coronal heating that determines the
plasma thermodynamics.

The active region run also allows us to analyze the force-freeness in a real-
looking solar atmosphere, which would be extremely difficult to do for the
real Sun. The first 2 Mm above the photosphere is clearly none-forcefree and
the current vector may be pointed to any angle with respect to the magnetic
field. The current and magnetic vectors become indeed more inline in the low
$\beta$ corona. However, even with $\beta<10^{-3}$, the magnetic force never
vanishes and all over the corona, there are none negligible points where the
current presents a greater derivation from the magnetic field. This suggest
the important role of plasma forces for a more accurate modeling of magnetic
field on the real Sun (e.g., Zhu & Wiegelmann, 2018).

##  5 RESULTS III: Coronal Dynamics Illustrated by Emission Measures

Differential emission measures describe the distribution of plasma density
(along the line-of-sight) with temperatures, as given by the relation

$DEM(T)dT=\int\limits_{LOS}n_{e}^{2}(T)dl.$

The integral of DEM over temperature is known as emission measures, i.e.,

$EM(T)=\int\limits_{T}^{T+\Delta T}DEM(T)dT$

DEMs inverted from multi-channel EUV/X-ray images (e.g., Kashyap & Drake,
1998; Hannah & Kontar, 2012; Cheung et al., 2015) are very useful for
investigating thermal properties of coronal plasma. In MURaM simulations, DEM
integrated along the 3 axes are computed during runtime and available as a
direct output at high cadence.

To evaluate the DEM in the simulation, the logarithmic temperature range
starting from $\log T{=4.5}$ to the highest temperature in the domain (time-
dependent) is discretized into a certain number of temperature bins with an
interval of $\Delta\log T{=}0.1$. Then we calculate along a given line-of-
sight the sum of $n_{e}^{2}$ in corresponding temperature bins (with details
given by Rempel, 2017). The line-of-sight integrals are done along the $x$,
$y$, and $z$ axes, respectively, and provide emission measures (in unit of
cm-5) of the transition region and corona from two side views and one top
view. We note that because this calculation implies an integral over
$\Delta\log T$, we use the term _emission measure_ instead of DEM hereafter.

Although intriguing dynamics spreads over the whole evolution, we select three
short time periods to briefly illustrate features and dynamics in the corona.
The first period covers the time when active region magnetic flux first
appears in the photosphere and active regions start to form. The second period
is in the early stage of active region formation, where newly emerged magnetic
flux drives a drastic reorganization of coronal magnetic fields on the large
scale. The third period is in the later half of the evolution, when active
regions are well developed.

###  5.1 Coronal Dimming above Emerging Active Regions

Figure 14: Coronal dimming created in the quiet Sun by an emerging active
region. (a): $B_{z}$ at the photosphere, before the active region emerges.
(b): $B_{z}$ as in Panel (a), but after the flux bundle (at $y{\approx}50$Ṁm)
has reached the photosphere and the flux concentration (at (x,y)=(120,50))
start to form. (c) & (d): Synthetic AIA 193 channel images from the top view
(i.e., along the $z$-axis). The time stamps are the same as in (a) and (b),
respectively. The dimming region is visible in the lower half of the field-of-
view. (e) & (f): synthetic AIA 193 channel images from a side view along the
$y$-axis. The time stamps are the same as in (a) and (b), respectively. The
dimming region spans between $x{\approx}30$ and 150 Mm and reaches
$z{\approx}80$ Mm.

We investigate how the first active region emerged in the quiet Sun changes
the corona. Figure 14 presents a comparison of the photospheric magnetogram
and coronal EUV images before (i.e., the quiet Sun) and after the first
emerging active region appears in the photosphere.

The magnetic field in the quiet Sun is a mixed-polarity field maintained by a
small-scale dynamo in the convection zone part of the simulation. The flux
imbalance noticeable on spatial scales of about 20 to 30 Mm is a manifestation
of the convection in the deeper convection zone. The corona above the quiet
Sun, as we have seen earlier in this paper, shows very smooth structures on
the large scale, with embedded small-scale bright points.

Panels (b) and (d) in Figure 14 show the magnetogram and synthetic AIA 193
channel image when the first active region appears. A flux bundle that spans
from $x{\approx}50$ to 110 Mm at $y{\approx}50$ Mm is breaking through the
photosphere in the area . The flux concentration at $(x,y){\approx}(120,50)$,
which is still rather small in size, already has a field strength well above 1
kG.

The difference in the corona is more evident. The EUV count rate in a small
area right above the forming sunspots is strongly reduced from a few hundred
to almost zero. More interestingly, we can see around the forming sunspots a
large area that becomes significantly dimmer than the quiet Sun. This dimming
region is over 100 Mm wide along the $x$-axis. The presence of this dimming
region is clearer when comparing Panels (e) and (f) . In particular Panel (f)
demonstrates that the dimming region is like a large bubble that already
reaches 80 Mm above the photosphere at this moment. This implies that even in
the very early stage when the sunspots seem to be very small, the emerged
magnetic field can already impact a much larger domain.

Figure 15: Evolution of the coronal thermodynamics in the dimming region over
an emerging active region. (a) & (b): Emission measure weighted temperature
from a side view along the $y$-axis, as defined in Equation (16). The time
stamps are the same as in Figure 14. (c) & (d): Mean density derived from
emission measures, as in Equation (17). (e): Change of the mean temperature
and density of the corona above the emerging active region (averaged within
the boxes marked in Panels (a)–(d)). The solid and dotted lines show the mean
temperature and density, respectively. (f): Change of the mean emission
measure (EM) in the corona above the emerging active region (averaged within
the boxed marked in Panels (a)–(d)). The blue line shows the mean emission
measure in $\log T{=}[6.0,6.3)$, and the red line is for $\log T{=}[6.3,6.5)$.
The black dotted line shows the unsigned magnetic flux in the photosphere.

The dimming region associated with an expanding magnetic field shows some
similarity to dimming regions associated with CMEs, which are caused by
evacuation of mass (e.g., Cheng et al., 2012; Dissauer et al., 2018). However,
this is NOT the case for the dimming region here. We use emission measure to
investigate the change of plasma properties during this period.

We assign the emission measure as weights and evaluate the mean temperature,
$\langle T\rangle$, by

$\log\left(\langle T\rangle\right)=\frac{\sum\limits_{i=1}^{nbins}EM(\log
T_{i})\,\log T_{i}}{\sum\limits_{i=1}^{nbins}EM(\log T_{i})}.$ (16)

$\langle T\rangle$ represents a mean temperature along a line-of-sight, and
more probably reflects the temperature of plasma that contribute the most
emission.

The sum of the emission measure over all temperatures represents the integral
of $n_{e}^{2}$ along the line-of-sight. Therefore, we can estimate the mean
density, $\langle n_{e}\rangle$, by

$\langle n_{e}\rangle=\sqrt{\sum\limits_{i=1}^{nbins}EM(\log
T_{i})\bigg{/}L},$ (17)

where $L$ denotes the length of the line-of-sight.

The maps of $\langle T\rangle$ and $\langle n_{e}\rangle$ for a side view
along the $y$-axis are presented in Figure 15. Panels (a) and (b) compare
$\langle T\rangle$ before and after the active region appears. The mean
temperature of the quiet Sun corona is around 1.5 MK and there is a large area
with even cooler temperatures. When the active region appears, the most
evident feature in the mean temperature is a dome of about 2 MK. This high
temperature dome extends for over 100 Mm in the horizontal direction and
reaches 80 Mm above the photosphere, which exactly matches with the dimming
region seen in the synthetic EUV images. Meanwhile the mean temperature in the
ambient corona is also slightly increased. This phenomenon is a clear evidence
that emerging active regions provides a significant excessive energy flux that
heats the quiet sun corona to higher temperatures, and this heating is already
present and impacts the large scale corona in the very early stage of active
region formation.

Panels (c) and (d) illustrate the change of $\langle n_{e}\rangle$. The quiet
Sun corona shows a mean density on the order of $10^{8}$ cm-3. The newly
emerged active region pushes plasma upwards from the lower atmosphere and the
mean density at $z{\approx}20$ Mm is increased by about an order of magnitude.
We can also discern a very slight decrease of the density in some areas at
$z{\approx}60$ Mm, which implies that the expanding magnetic field is
reorganizing pre-existing plasma in the corona, and perhaps it is driving the
plasma sideways.

For the region above the forming active region as marked by the boxes in
Figure 15, we plot the mean temperature and density as functions of time in
Figure 15(e), and Panel (f) plots the temporal evolution of coronal emission
measures in two adjacent temperature ranges. The changes of mean temperature
and density are consistent with the expectation that the forming active region
gives rise to more intensive heating and pushes masses into the corona. The
most important information is given by the emission measure. In the first
seven hours, the coronal emission measure is mostly in the temperature range
of $\log T{=}[6.0,6.3)$. However, the emission measure in a higher temperature
range keeps growing with the growth of magnetic flux in the photosphere.
During a period of one hour at $t{\approx}8.5$ h, the hotter emission
drastically increases and the cooler emission significantly drops.

To conclude, the coronal dimming region found during the early stages of
active region formation is a manifestation of intense heating that evolves the
quiet Sun corona into an active region corona. So-called dark moat, a dimmed
area around developed isolated active regions, has been found for decades (see
Singh et al., 2021, and references therein). While recent observations
revealed that emerging active regions give rise to emerging dimming regions
(Zhang et al., 2012a; Payne & Sun, 2021) or a large-scale outwards-propagating
dark ribbon (Zhang et al., 2015). The response of the quiet Sun to the
emerging active region in this simulation is consistent the phenomenon
observed by Payne & Sun (2021), particularly the evolution of the emission
measure of the dimming region, and confirms that the emerging dimming are
caused by enhanced heating in the corona.

###  5.2 EUV Emission above Emerging Active Regions

Figure 16: Coronal emission measure in temperatures from $\log T{=}4.5$ to
$\log T{=}7.0$ in the early stage of the simulation. Each panel displays the
total emission measure in a selected temperature range divided the by the
number of bins it contains (with a bin interval of $\Delta\log T{=}0.1$). The
corresponding temperature range and number of bins are noted above each panel.
This figure is a snapshot of an animation that covers about 2.5 hours of
evolution (available in the electronic version of this paper). This snapshot
shows cool plasma in the lower atmosphere being ejected into the corona after
an eruption occurs. The rest of the animation also presents how the ejected
plasma falls back to the lower atmosphere, as well as the formation of many
long coronal loops, which is described in more detail in Section 5.2.

In Figure 16 and the associated animation, we present a side-view (along the
$y$-axis) of the coronal emission measure in a time period of about 2.5 hours
starting at about $t{=}11\,$h. This period is in the early stage of the full
active region evolution, when fast flux emergence is ongoing and active
regions start to form. The emission measure during this time period highlights
a fast changing corona, where we see a large scale eruption of cool material
from the low atmosphere, and formation of large-scale coronal loops and
enhancement of coronal temperature as consequences of this eruption.

The emission measure calculated at runtime has an interval of $\Delta{\rm
log}T=0.1$ and covers temperatures ranging from $\log T{=}4.5$ to the highest
temperature in the domain. For presentations in this section, we re-bin the
emission measure in neighboring bins and reform the emission measure array to
less bins. Figure 16 is a snapshot of the animation at $t{=}$11 h 44 m, and
presents the emission measure that has been re-binned into six bins. This
snapshot captures a large-scale eruption, which is the most interesting
highlight in this time period.

####  5.2.1 Large Scale Eruption

The animation shows that the eruption starts at $t{=}$10 h 59 m, when the
emission in the hottest bin suddenly increases and a hot loop appears at
$(x,z){\approx}(80,10)$. Right after that time, large amount of cool plasma
($\log T{<}6.0$) is ejected from the low atmosphere and reaches over 80 Mm in
height in about 6 minutes. Large parcels of ejected plasma quickly destruct
into small fragments that exhibit various behaviors. The falling plasma
behaves like coronal rains (e.g., Schrijver, 2001; Antolin & Rouppe van der
Voort, 2012; Antolin et al., 2015) and traces the magnetic field lines.
However, we also see many fragments that seem to “float” in the corona and can
dynamically move up and down. Their movement illustrates the dynamical
evolution of the magnetic field over emerging active regions, which was not
considered by other numerical models of coronal rains (e.g., Müller et al.,
2005; Xia et al., 2017; Kohutova et al., 2020). Moreover, the falling plasma
may undergo a cooling process as in a typical coronal rain scenario, or
heating process that converts the cool plasma into a mass source to hotter
coronal loops.

We also note that the magnetic field in the corona of our simulation does not
form a dip, which was suggested to be the favorable geometry for filament
formation (Mackay et al., 2010) and has been tested by numerical experiments
(e.g., Xia & Keppens, 2016). Otherwise, we could very likely see a
filament/prominence forming by these falling cool plasma elements.

####  5.2.2 Fast Expanding Coronal Loops

In a period of about 20 minutes after $t{=}11$ h 20 m, the falling plasma
gradually fills a large magnetic arcade rooted at $x{\approx}0$ and 130 Mm,
and this forms a large coronal loop mostly visible around $\log T{=}6.0$. If
we compare the panels of $\log T{=}[5.5,6.0)$ and $\log T{=}[6.0,6.2)$,
patches that disappear in the former will show up in the latter. After
$t{=}12$ h, the large loops almost disappear in the bin of $\log
T{=}[5.5,6.0)$, but remain visible in the bin of $\log T{=}[6.2,6.4)$.
Meanwhile, another group of loops emerges in the area between $x{\approx}70$
and 130 Mm. They first become visible in the lower temperature bins ($\log
T{<}6.0$) at about $t{=}11$ h 30 m, and outline magnetic field lines that are
expanding upward. At $t{=}12$ h, these loops completely fade away in the bin
of $\log T{=}[4.5,5.5)$, and are more evident in the bin of $\log
T{=}[6.0,6.2)$. About 35 minutes later, they can only be seen in bins of $\log
T{>}6.0$, and have expanded to a height of about 100 Mm. Finally, these loops
disappear in the bin of $\log T{=}[6.0,6.2)$ and are only visible in the bin
of $\log T{=}[6.2,6.4)$.

The shift of similar structures from bins of lower temperatures to those of
higher temperatures, as illustrated by the two examples here, provides a clear
evidence of heating coronal plasma during fast flux emergence. The heat input
can efficiently heat the cool plasma carried in emerging magnetic field lines
to temperatures well above one MK and maintain the density at loop apexes for
a long time. Another interesting behavior of loops formed in the fast
emergence stage is that they are constantly changing. As we see in the
animation, a loop or thread that is bright in a certain bin may become fade
and disappear in a few minutes, and it very likely to be found a bin of higher
temperatures. Furthermore, if we observe a structure that seems to remain at
one place, this does not necessary imply a static state. However, it is
probably a combined effect of old threads disappearing (moving to higher
temperatures) and new threads forming (heated from lower temperatures).

###  5.3 EUV Emission above Developed Active Regions

Figure 17: Coronal emission measure in temperatures from $\log T{=}4.5$ to
$\log T{\geqslant}8.0$ in a late stage of the simulation. The quantity shown
in each panel is the total emission measure in the temperature range divided
by the number of bins (with a bin interval of $\Delta\log T{=}0.1$), which are
noted above each panel. This figure is a snapshot of an animation that covers
about one hours of evolution (available in the electronic version of this
paper). This snapshot shows the corona right after a high speed jet occurs.
The jet, which is visible at this moment in the middle and bottom rows ($\log
T{\geqslant}6.2$), is propagating from $(x,z){\approx}(80,5)$ to the top-left
corner of the field-of-view. The animation also highlights the evolution of a
set of cooling coronal loops gradually appearing in the $\log T{=}[6.2,6.4)$
panel and transient and stable hot compact loops $\log T{\geqslant}6.8$ in the
active region core, as describe in Section 5.3.

We present a side-view (along the $y$-axis) of the coronal emission measure in
a late stage of the full evolution, when flux emergence in the photosphere
slows down and active regions are well developed. The time period presented
here starts at about $t{=}30.5$ h and lasts for about one hour. The emission
measure features a high speed jet propagating through the vertical domain and
dynamics of typical coronal structures above developed active regions, such as
small-scale surges of cool plasma, warm coronal loops arching over active
regions, and hot compact loops in the active region core.

In this time period, the hottest plasma in the corona is found in temperatures
of several tens MK, which is significantly higher than the highest
temperatures in the early stage. Therefore, we reform the original output of
the emission measure to nine temperature bins and present a snapshot of the
re-binned emission measure at $t{=}$30 h 40 m in Figure 17. This figure is
also a frame of the animation showing the temporal evolution during this time
period.

####  5.3.1 Fast Propagating Hot Jet

Jets are ubiquitously observed on the Sun from the chromosphere to the corona
(Raouafi et al., 2016; Shen, 2021). An intriguing event seen during this
period is a fast propagating jet coming out from the core of the active
region. The jet is first seen in the hottest bin of $\log T{\geqslant}7.3$. It
can be identified in frames around $t{=}30$ h 39 m 40 s as a thin thread along
a line pointing from $(x,z){\approx}(80,0)$ to $(x,z){\approx}(60,60)$. The
jet quickly disappears in the hottest bin in less than 20 s and appears
subsequently in the bin of lower temperatures.

The snapshot of the animation at $t{=}30$ h 40 m 45 s is displayed in Figure
17. At this moment, we can see in the bin of $\log T{=}[7.0,7.3)$ a trace that
extends from its onset location to the upper left corner of the field-of-view.
In the animation the trace shows up almost instantaneously over a large
distance, yielding an apparent propagation speed well over 1000 km/s. Thus it
would be very interesting to clarify by a further analysis (not performed in
this study) whether this is an actual flow or a heating front. Meanwhile, the
lower part of the jet can be observed in bins of intermediate temperatures
(middle row in Figure 17) with a similar shape and propagation direction as it
first appeared in bins of much higher temperatures. The propagation speed of
the structure seen in these temperatures is about 500 km/s, which is a fast
event compared with observed EUV/X-ray jets (Shimojo et al., 1996; Savcheva et
al., 2007).

The animation shows that, in few minutes after the time instance shown in
Figure 17, the jet gradually cools down and becomes visible in lower
temperatures (bins of $\log T{<}6.0$). As it cools down the jet expands in the
transverse direction, thus the structure seen in the cool bins is
significantly wider than that seen in the hot bins. Eventually, after about
$t{=}30$ h 45 m, the ejected plasma gradually falls back down.

Jets triggered by flux emergence have been extensively studied by numerical
simulations in the past decades(e.g., Yokoyama & Shibata, 1995; Archontis et
al., 2005; Moreno-Insertis & Galsgaard, 2013; Fang et al., 2014; Iijima &
Yokoyama, 2017, and references therein). The jet seen in our simulation occurs
in an unprecedentedly complex environment as compared with previous
simulations. It is interesting to study how the configuration of ambient
magnetic field affects the generation and propagation of the jet, as suggested
by Pariat et al. (2015).

####  5.3.2 Warm Coronal Loops

Loops are the most common feature observed in the corona (Reale, 2014). The
middle row of Figure 17 shows prominent coronal loops that are observed in
intermediate temperatures (warm loops) as arches reaching tens of Mm over the
solar surface. At the beginning of the animation, we can see several loops
with apexes at $z{\approx}40$ to 60 Mm in the bin of $\log T{=}[6.2,6.4)$. At
the same time, it is more difficult to isolate a single loop from the bin of
$\log T{=}[6.4,6.6)$, but there are many thin threads above the rather
diffusive background. In the bin of $\log T{=}[6.6,6.8)$, we can identify a
wide loop-like arch that roots at $x{\approx}45$ and 125 Mm and reaches about
$z{\approx}$70 Mm at apex.

In the following 30 minutes, the loop-like arch in the hotter bin in the
middle row slightly decays and several bright loops gradually appear as whole
arches in the cooler bin. At $t{=}31$ h, these loops can be seen in the cooler
bin as nice arches with apexes at $z{\approx}$60 to 80 Mm and feet separations
of about 80 Mm. These geometry properties are consistent with the structure
that was identified earlier in the hotter bin. All these behaviors give a
clear evidence that the new loops seen in the cooler bin result from pre-
existing structures cooling down from higher temperatures. In the next 30
minutes after $t{=}31$ h, the loops vary constantly as existing bright threads
may disappear and new threads also form. But the overall structures in the
bins in the middle row of Figure 17 do not change significantly.

The cooling behaviors are very often observed in real coronal loops (e.g,
Schrijver, 2001; Winebarger et al., 2003; Winebarger & Warren, 2005; Viall &
Klimchuk, 2012; Li et al., 2015; Froment et al., 2015). This have been
suggested as an evidence of impulsive heating in coronal loops and an
explanation on warm loops whose observational properties are inconsistent with
predictions of hydrostatic equilibrium (Aschwanden et al., 2001; Warren et
al., 2002, 2003; Klimchuk et al., 2010). While in our simulation, non-
equilibrium is the most common state of coronal loops over actively evolving
active regions.

####  5.3.3 Hot Loops in Active Region Cores

As the last part of this section, we examine the emission from the hottest
plasma in the corona, which is presented in the bottom row of Figure 17. We
first see a significant amount of emission by hot plasma in the high corona,
as shown in the bin of $\log T{=}[6.8,7.0)$. This is actually caused by the
closed top and (periodic) side boundaries that do not allow for energy flux
leaving the domain. Meanwhile a very high energy input in the corona is
provided by the strong active regions that contain already $10^{23}\,$Mx of
unsigned flux at this moment. The active regions in this simulation are
comparable to the largest active regions seen in the real Sun in terms of
magnetic flux, and real active regions would not be confined in a limited
volume as in this simulation. Therefore the super-hot emission in the high
corona near the top of the simulation domain may not fully represent the real
corona.

In the lower corona, short coronal loops with apexes generally below 20 Mm
give rise to significant emission measure in temperatures around $10^{7}$ K.
These hot loops can be stable or transient structures. For instance, the
stable structure appears as vertical “plume” at $x{=}120$ Mm in the bin of
${\rm log}T{=}[7.0,7.3)$ is a loop that extends mostly along the $y$ direction
(i.e., the line-of-sight for this figure) and connects two strong sunspots at
$y{\approx}60$ and 90 Mm, respectively. This loop is formed within the two
sunspots and remains hotter than 10 MK till the end of the evolution. The
animation of Figure 17 shows many transient hot loops with life times ranging
from a few seconds to tens of minutes.

These hot loops in the active region core resemble hot plasma revealed by
X-ray observations of solar active regions (Schmelz et al., 2009; Reale et
al., 2009; Warren et al., 2012). The core of active regions in the simulation
is highly structured and the surrounding becomes more diffusive, which is
similar to the observation reported by Reale et al. (2007). The presence of
short-lived and steady features is also found in observation Warren et al.
(2010).

Whether a impulsive or steady heating dominates in the active region cores is
not fully clear in observations (see e.g., Tripathi et al., 2011; Warren et
al., 2010; Winebarger et al., 2011, and references therein). However, it is
very likely that in the simulation both mechanisms work. Chitta et al. (2018)
proposed that cancellation of opposite polarity magnetic field can provide
sufficient energy for hot loops. This could be a very viable mechanism for the
transient hot loops in this simulation. In particular, the ongoing flux
emergence continuously brings small patches of magnetic field that can trigger
cancellation with pre-existing magnetic field. However, the long-living and
steady hot loops that last for over 20 hours is more likely to be maintained
by a (statistically) steady heating. This leaves interesting open questions
for more in-depth analysis of the simulation on the energy release in active
region core.

##  6 SUMMARY

We presented the first results of a realistic radiative MHD simulation of the
Sun’s uppermost convection zone and atmosphere. The control equations of this
simulation consider radiative transfer in the lower solar atmosphere, and
account for optically thin radiation and field-aligned thermal conduction in
the corona. This allows for realistic modeling of the thermodynamics of
coronal plasma, hence a direct and quantitative comparison between synthetic
observables and actual observations. The large vertical domain spanning from
the upper most convection zone to over 100 Mm in the corona allows this
simulation to self-consistently account for the coupling between difference
layers of the solar atmosphere though magnetic field, mass and energy flows in
a single simulation domain. We studied the self-consistently maintained quiet
Sun corona and the formation of strong solar active regions though magnetic
flux emergence. The results presented in this paper are summarized as
following.

###  6.1 General Properties of the Quiet Sun and Active Regions

The quiet Sun run reaches a dynamical equilibrium, where the vertical profiles
of temperature and density are consistent with observations. In particular,
the mean coronal temperature, which remains stably at about 1.5 MK, is
maintained by dissipating the Poynting flux that is naturally generated by the
magneto-convection at the solar surface. In addition to the stable background
quiet Sun corona, we also demonstrate various small-scale activities in the
quiet Sun. The most evident ones are the numerous EUV bright point seen in the
synthetic observations and the propagating disturbances as revealed by the
running difference ratio of synthetic AIA images.

The active region run considers the emergence of magnetic flux from the upper
most convection zone to the corona and the formation and evolution of solar
active regions in a time period of about 50 hours. The setup of active-region
scale magnetic flux is a key ingredient of the active region simulation.
Unlike the idealized flux tubes/sheets that were extensively used in previous
studies, our simulation adopts the flux bundles generated in a global-scale
solar convective dynamo (Fan & Fang, 2014). The significantly improved realism
in the magnetic field leads to more realistic sunspots. Moreover, these flux
bundles also give rise to complex active regions that turn out to be very
flare productive. The eruptions occurred in the simulation will be presented a
forthcoming separate paper.

The evolution of the sunspots is consistent with the results presented in
Paper I that focused on the formation of sunspots in the photosphere and their
asymmetric properties. The simulations in Paper I employed the same setup but
with the top boundary placed at about 600 km above the photosphere. Therefore,
we can confirm that the results in that study were not affected by the
boundary.

###  6.2 Coronal Magnetic field and its energy partition

We focused more on the evolution of the coronal plasma and magnetic field
above active regions in the present study. The total unsigned photospheric
magnetic flux grows to more than $10^{23}$ Mx, which makes the simulated
active regions comparable to the largest active regions observed on the real
Sun. The total magnetic energy in the volume of about $197^{2}\times 113$ Mm3
increases from $8.5\times 10^{30}$ erg for the quiet Sun to $4.7\times
10^{33}$ erg at the end of the evolution. About 18% of the total magnetic
energy is free energy. This ratio remain relative constant during the
evolution of the simulated active region. we investigated the distribution of
magnetic energies with height for three evolution stage. The free energy
decays slower than the total magnetic energy. As a results, the magnetic field
in the upper half of the domain (i.e., $z{>}60$ Mm) is highly non-potential,
in the sense that the free energy contributes 80% to almost 100% in the total
magnetic energy.

This active region model gives the first _r_ ealistic presentation of the
force-freeness and plasma $\beta$ in a single domain covering two extreme
regimes. The magnetic field and plasma in this simulation are evolved
consistently from the beginning, and the magnetic force, plasma force, and
viscous stress find their balance spontaneously. We assessed the force-
freeness of the magnetic field by the angle between the current and magnetic
vectors. Beneath the photosphere, the angle is generally larger than 40∘ with
the most probable value found at almost 90∘. This reflects the simple fact
that plasma dominates over the magnetic field in the convection zone. In the
solar atmosphere, the angle quickly switches to the other end of the axis, and
it is generally smaller than 10∘ with the most probable value found at about
5∘ for the most part of the corona.

The implications are twofold. On one hand, the Lorentz force (i.e., the angle
between the current and magnetic field) can not vanish in the real corona.
Even in a purely quiescent area, it is necessary to have a residual Lorentz
force to balance the pressure gradient and to release energy as part of the
coronal heating process, which implies a balance with acceleration and viscous
forces. As shown in our simulation, the current is rarely parallel to the
magnetic field, and the distribution has an enhanced wing on the right side
(larger than) of the most probable value. On the other hand, for most of
corona volume, the angle remains relatively small, which suggests that the
force-free field seems to be an acceptable assumption.

The distribution of the plasma $\beta$ follows the same scenario. Beneath the
photosphere, $\beta$ is in the range between $10^{2}$ and $10^{6}$. While
above the photosphere, $\beta$ rapidly drops to below unity. The lowest (mean)
value of about $10^{-3}$ is seen in the transition region, and the local
minimum can be smaller than $10^{-5}$. A detailed view of the vertical domain
from 5 Mm below to 5 Mm above the photosphere shows that the transition from
the highly non-forcefree and high-$\beta$ regime to the (close to) forcefree
and low-$\beta$ regime occurs in a height range of about just one Mm.

###  6.3 Dynamics of Coronal Plasma

We highlighted three time periods in different evolution stages of active
regions and demonstrate coronal features that are comparable to structures and
activities observed in the real corona.

When the first active region starts to appear in the photosphere, we found in
the synthetic AIA images a bubble-like dimming region that slowly expands from
the active region to almost the whole corona. We investigated the emission
measures and thermodynamics properties in the coronal above the forming active
region. Both the mean temperature and density in the corona increases, and
this trend coincides with the growth of photospheric magnetic flux.
Consequently, the emission in the temperature range of the original quiet Sun
corona drops by an order of magnitude, which leads to the dimming region. As a
comparison, the synthetic images of a filter that is sensitive to higher
temperatures would demonstrate a brightening above the forming active region.

We demonstrated two eruptive activities that occur in the early and late
stages of active region formation, respectively. In the early stage when the
active regions are forming through fast flux emergence, a large-scale eruption
ejects a significant amount of plasma from the lower atmosphere to over 100 Mm
in the corona. This eruption provides an instance of mass circulation on the
large scale. Ejected plasma with temperature around $10^{5}$ K may fall back
to the solar surface along magnetic field lines, or get heated and becomes a
mass source for hotter coronal loops. Some clumps could “float” up and down in
the corona, as they are supported by the magnetic field.

A fast propagating jet is observed in the late stage when the active regions
have developed strong and complex magnetic field. The speed of the bulk of the
jet is over 500 km/s, and the apparent speed of the leading front of the jet
exceeds 1000 km/s. The plasma in the source region of the jet is heated to
extremely high temperature (${>}20$ MK), which is strong evidence of magnetic
reconnection.

We highlighted different types of coronal loops through the evolution. In the
early stage,dynamical and transient large-scale coronal loops quickly form
after the eruption. During the same time, the entire corona undergoes a rapid
increase in temperature and evolves from the quiet Sun to a highly structured
active region corona.

In the late stage when active regions are more developed, a group of long and
warm loops (1–3 MK) forms though cooling of loops are heated to much higher
temperatures (${>}$5 MK). Individual bright thread within the group of loops
appears and disappears in a very short time. However, if we exam the group of
loops as a whole structure, it is sustained for a life-time that is much
longer than the cooling time scale. We found in the active region cores both
stable and transient hot loops in temperatures around 10 MK. These loops
appear to be rather independent structures and implies the presence of
impulsive and stable heating at the same time.

##  7 Conclusion

We present in this paper an introduction on the method and first results of a
comprehensive model of the quiet Sun and solar active regions. This model
provides in a single simulation an unprecedented combination of a
computational domain that is wide enough to host large active region and spans
from the uppermost convection zone to the corona, a billion grid points mesh,
a long temporal evolution including about 11.5 hours of quiet Sun and 48 hours
of active region, strong and complex active region scale magnetic flux, and
sophisticated physics for radiative and conduction energy transport in the
solar atmosphere, each of which is truly a challenge for numerical
simulations.

As summarized above, the simulated quiet Sun and active region can reproduce
the realistic stratification of plasma density and temperature, and many of
the features and dynamics in observations. In particular, this model provides
a fully self-consistent evolution of the plasma and magnetic field and allows
for direct connection between the emission, plasma thermodynamics, and
magnetic field that are difficult to measure in remote sensing observations.
Moreover, the simulation data can also be used to verify results from methods
that assumes the domination of either plasma or magnetic and neglect the
other.

The scope of the current paper is to present overview of the general
properties of the simulated quiet Sun and active regions, with highlighting a
few results that could lead to a wide range of following studies. In addition
to the very wide range of topics covered in this paper, more than 170 flares
occur in 48 hours of the active region stage, with the largest flare reaching
M class and accompanied by a spectacular coronal mass ejection. The overview
of the results related to solar eruptions will be presented in a forthcoming
paper.

The simulation generated in total amount about 400 TB data, including full 3D
snapshots of the whole computational domain, 2D diagnostic slices (plane cuts,
$\tau$ surfaces, and emission measures as presented in this paper) at high
temporal cadence, and the presented results are from the analysis of less than
10% of them. The data will be open to the community. However, due to the
limited manpower and IT resources, these data cannot be made available as a
whole package. For the moment, interest of using the data can be addressed to
authors of this paper.

We would like to acknowledge high-performance computing support from Cheyenne
(doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information
Systems Laboratory, sponsored by the National Science Foundation. The National
Center for Atmospheric Research (NCAR) is sponsored by the National Science
Foundation. F.C. is supported by the Fundamental Research Funds for the
Central Universities under grant 0201-14380041 and had been supported by the
Advanced Study Program postdoctoral fellowship at NCAR and by the George
Ellery Hale postdoctoral fellowship at the University of Colorado Boulder