Issue 
A&A
Volume 652, August 2021



Article Number  A32  
Number of page(s)  13  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202040192  
Published online  04 August 2021 
Stellar Xrays and magnetic activity in 3D MHD coronal models
Max Planck Institute for Solar System Research,
37077
Göttingen,
Germany
email: zhuleku@mps.mpg.de
Received:
21
December
2020
Accepted:
28
May
2021
Context. Observations suggest a powerlaw relation between the coronal emission in Xrays, L_{X}, and the total (unsigned) magnetic flux at the stellar surface, Φ. The physics basis for this relation is poorly understood.
Aims. We use threedimensional (3D) magnetohydrodynamics (MHD) numerical models of the coronae above active regions, that is, strong concentrations of magnetic field, to investigate the L_{X} versus Φ relation and illustrate this relation with an analytical model based on simple wellestablished scaling relations.
Methods. In the 3D MHD model horizontal (convective) motions near the surface induce currents in the coronal magnetic field that are dissipated and heat the plasma. This selfconsistently creates a corona with a temperature of 1 MK. We run a series of models that differ in terms of the (unsigned) magnetic flux at the surface by changing the (peak) magnetic field strength while keeping all other parameters fixed.
Results. In the 3D MHD models we find that the energy input into the corona, characterized by either the Poynting flux or the total volumetric heating, scales roughly quadratically with the unsigned surface flux Φ. This is expected from heating through fieldline braiding. Our central result is the nonlinear scaling of the Xray emission as L_{X} ∝ Φ^{3.44}. This scaling is slightly steeper than found in recent observations that give powerlaw indices of up to only 2 or 3. Assuming that on a real star, not only the peak magnetic field strength in the active regions changes but also their number (or surface filling factor), our results are consistent with observations.
Conclusions. Our model provides indications of what causes the steep increase in Xray luminosity by four orders of magnitude from solartype activity to fast rotating active stars.
Key words: Sun: corona / stars: coronae / magnetohydrodynamics (MHD) / methods: numerical / Xrays: stars
© J. Zhuleku et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
1 Introduction
Stellar coronal Xray emission is observed to increase with stellar rotation rate (e.g., Pizzolato et al. 2003; Wright et al. 2011; Reiners et al. 2014; Magaudda et al. 2020). It is widely assumed that an increase in rotation could be responsible for stronger dynamo action leading to larger surface magnetic field. Some active stars (e.g., M dwarfs) that rotate rapidly (typical periods of 1–2 days) show high average photospheric magnetic field strengths which can reach up to 8 kG or even more (Reiners 2012). Because of this increased photospheric magnetic field, we can expect that a stronger upwarddirected Poynting flux is generated that can heat the corona to higher temperatures, leading to stronger Xray emission in the corona. An indication of such behavior has been found in stellar observations (e.g., Vidotto et al. 2014) revealing a close relation between coronal Xray emission and surface magnetic flux.
The scaling relationship between the coronal Xray emission L_{X} and the surface magnetic flux Φ has been extensively studied using solar and stellar observations, and follows a power law, L_{X} ∝ Φ^{m}. In early studies, the powerlaw index m was found to be close to unity (Fisher et al. 1998; Pevtsov et al. 2003), that is, the Xray radiation scales almost linearly with magnetic flux. However, more recent studies suggest a much steeper power law with m = 1.80 (Vidotto et al. 2014) or even steeper with m = 2.68 (Kochukhov et al. 2020). The physical mechanism relating the observed Xray emission to the surface magnetic flux is still under debate.
To study the impact of the surface magnetic field on the coronal Xray emission in the environment of a realistic setup, the use of 3D magnetohydrodynamic (MHD) models is required. In addition, the 3D numerical simulations will provide a useful tool with which to further test the validity of a simplified analytical model. The main advantage of 3D models is the selfconsistent treatment of the corona. The heating originates from the Ohmic dissipation of currents induced by photospheric magnetoconvective motions. This drives the magnetic field similar to Parker’s field line braiding (or nanoflare) model (Parker 1972, 1983). The Parker fieldline braiding model has been extensively studied in numerical models. The energy was shown to cascade in current sheets from large scales to dissipative scales before converting to heating (Rappazzo et al. 2008). In addition, 3D MHD simulations of footpoint motions have proven successful in forming a selfconsistent corona (see e.g., Gudiksen & Nordlund 2002, 2005a,b; Bingert & Peter 2011, 2013; Hansteen et al. 2015; Dahlburg et al. 2016, 2018; Matsumoto 2021).
The numerical models are able to provide the necessary energy flux in the corona to heat it to temperatures beyond 1 MK and this energy flux is consistent with observations (Bingert & Peter 2011; Hansteen et al. 2015). Furthermore, extreme ultraviolet (EUV) synthetic spectra from these 3D simulations can explain some aspects of the actual observations (Peter et al. 2004; Dahlburg et al. 2016; Warnecke & Peter 2019a). This confirms the validity and efficiency of Parker’s field line braiding model to create a hot corona. These models can also be used to study the effects of magnetic helicity injection in the photosphere of active stars on the resulting coronal Xray emission (Warnecke & Peter 2019b). This showed that an increase of photospheric magnetic helicity without changing the vertical magnetic field increases the coronal Xray emission following simple powerlaw relations. However, the effect of the surface magnetic activity on the coronal Xray emission in 3D MHD models of solar and stellar coronae has not yet been studied.
In our study, we focus on the effect that the photospheric magnetic field strength has on the coronal Xray emission. This is motivated by the observation that stars more active than the Sun host a stronger surface magnetic field. For that reason, we choose to increase the strength of the vertical surface magnetic field at the bottom boundary of our computational domain; that is, we treat the peak (or average) magnetic field strength as a free parameter. All the other parameters remain the same in all numerical experiments. By varying only one parameter (i.e., the surface magnetic field) we can study the exact relation between magnetic flux and coronal emission. Other parameters important for the coronal energy input, such as the photospheric velocity distribution,are poorly constrained by observations for other stars and are therefore not changed (or varied) in the present work. Our main objective is to relate the synthetic Xray emission from the numerical models with the surface magnetic flux therein and compare this to the observed relationships. Furthermore, we compare our numerical results to our earlier analytical model (Zhuleku et al. 2020), which is briefly summarized in Sect. 2.
2 Analytical scaling relations
The dependence of the coronal Xray emission on the surface magnetic field has been extensively studied for the Sun as well as for other stars (Fisher et al. 1998; Pevtsov et al. 2003; Vidotto et al. 2014; Kochukhov et al. 2020). In Zhuleku et al. (2020), we developed an analytical model to describe the L_{X} ∝ Φ^{m} relation, where L_{X} is the coronal Xray emission and Φ the total surface unsigned magnetic flux.
Our model is based on the wellknown Rosner, Tucker & Vaiana (RTV; Rosner et al. 1978) scaling laws. These authors derived scalings relating the volumetric heating rate and the loop length to the coronal temperature and pressure. Alternatively, we can express the RTV scaling laws in a way that allows us to relate temperature T and number density n with the volumetric heating rate H and loop length L, (1) (2)
Using these scaling laws together with other relations, we derive an analytical expression for the Xray emission L_{X}, (3)
The powerlaw index m depends only on four parameters α, β, γ, and δ. The first, α, is related to the temperature sensitivity of the instrument used for the Xray observations. In general, the Xray radiation is proportional to the density squared and to a temperaturedependent function R(T). This function R(T) is the sum of all the contribution functions of emission lines and continua and is also known as the temperature response function, and differs from one instrument to the next. We found that the temperature response function for temperatures below log _{10}T[K] = 7 can be expressed as a power law (cf. Zhuleku et al. 2020, their Fig. 1), (4)
The parameters β and γ characterize the relation between the energy flux, or the vertical Poynting flux S_{z}, injected in the corona and the vertical surface magnetic field B and the volumetric heating rate H, again through power laws, (5) (6)
The case β = 1 represents Alfvén wave heating and β = 2 represents nanoflare heating (see also Sect. 6.1). The parameter γ was considered to be unity, γ = 1. Finally, δ relates the surface area covered by a magnetic structure (e.g. a whole active region) to the total magnetic flux, (7)
Solar studies suggest a value of δ = 0.819 (Fisher et al. 1998).
In our analytical study, we found the powerlaw index m to be in the range from roughly one to almost two (Zhuleku et al. 2020). This result agrees relatively well with most observations, but at least one more recent observation finds an even steeper powerlaw connection of L_{X} ∝ Φ^{m} with m of 2.68 (Kochukhov et al. 2020).
In the numerical study presented in the following sections, we assume that the area covered by the magnetic field remains the same, that is, the change of the surface magnetic flux is solely due to the (average) vertical magnetic field strength. This is equivalent to choosing δ = 0 in Eq. (7). In this case, we get a much steeper power law in the analytical model, with m up to 4 (Zhuleku et al. 2020, Sect. 5.2, Eq. (15)). The numerical study in this paper provides a more detailed comparison to observations that closely match those of Kochukhov et al. (2020).
3 Numerical model setup
3.1 Basic equations
Our model is based on the work of Bingert & Peter (2011, 2013). We use the PENCIL CODE (Pencil Code Collaboration 2021), where we numerically solve the MHD equations from the photosphere up to the corona. The PENCIL CODE is a finite difference code with a sixthorder numerical scheme used for spatial derivatives and thirdorder RungeKutta numerical scheme used for time derivatives. The size of the computational box is 128 × 128 × 128 grid points in Cartesian coordinates (x, y, z), representing a 50 × 50 × 50 Mm^{3} volume witha gridscale of 390 km in all directions. This comparably small grid size allows us to run a sufficiently large number of numerical experiments to study the relationship between the magnetic activity and coronal emission.
In principle, a higher resolution could affect the results, in particular because we are only barely resolving the granulation scale (on the Sun). However, earlier studies showed that reducing the grid spacing from around 550 km by about a factor of four gives a similar overall appearance of the corona, although the finer details become visible (Chen et al. 2014; Chen & Peter 2015). Based on these findings we consider the results we get at the grid spacing of 390 km to give a good representation of the injected photospheric Poynting flux then powering the hot corona. This is also supported by earlier studies which used a grid spacing comparable to that of our model. At comparable resolution, Gudiksen & Nordlund (2005b,a) found a loopdominated active region. In a more recent datadriven model, Warnecke & Peter (2019a) used a grid spacing very close to ours (350 km), albeit for a much larger computational domain, and found very good agreement between model and observations.
The MHD equations are the continuity, momentum, energy, and induction equation connecting density ρ, velocity u, and temperature T with the magnetic field B, and pressure p: (8) (9) (10)
Here the Lagrangian derivative D∕Dt is defined as D∕Dt = ∂∕∂t + u ⋅∇. The current density is given by j = ∇×B∕μ_{0}. The adiabatic index of an ideal gas is given by γ = 5∕3. We use a constant gravitational acceleration, g = (0, 0, −g), where g = 274 m s^{−2}, c_{V} is the specific heat capacity at constant volume, and is the rate of the strain tensor. The resistivity η and viscosity ν are constants through the whole computational box with values η = ν = 5 × 10^{10} m^{2} s^{−1}. At the bottom boundary we reduce η to prevent excessive diffusion of the photospheric magnetic field.
The magnetic field is calculated through the vector potential as, (11)
which ensures that the ∇⋅B = 0 is satisfied at all times. For our simulations we choose the resistive gauge ϕ = η∇⋅A so that the induction equation reads, (12)
To close the system of equations we need the equation of state for an ideal gas which connects the gas pressure with temperature, (13)
where ρ is the mass density, k_{B} is the Boltzmann constant, μ is the mean atomic weight and m_{p} is the proton mass.
The radiative losses L_{rad}, are modeled with an optically thin approximation using the radiative loss function P(T). The function P(T) depends on the abundances as well as the radiation processes that take place, such as the ionization rate and recombination rate (Meyer 1985; Murphy 1985; Cook et al. 1989). For a detailed discussion on the implementation of the radiative loss function in the code, see Bingert (2009).
The (Spitzer) heat flux along the magnetic field lines q reads, (14)
where K_{0} = 10^{−11} W(mK)^{−1} (Spitzer 1962) and is the unit vector of the magnetic field. To speed up the simulations, we replace Eq. (14) by a nonFourier heatflux scheme. See Warnecke & Bingert (2020) for a detailed description and discussion. Additionally, we use a semirelativistic correction to the Lorentz force similar to the work of Boris (1970), Gombosi et al. (2002) and Rempel (2017). For details of the implementation, see Chatterjee (2020) and Warnecke & Bingert (2020).
One concern for all existing 3D models including heat conduction is the limited resolution in the transition region. Because the heat conductivity gets very inefficient for low temperatures (∝ T^{5∕2}), the temperaturegradient has to become very steep in the transition region. This can necessitate a grid spacing along the magnetic field down to a km or less, which can be achieved in 1D models, in particular when using adaptive grids (e.g., Hansteen 1993) that then also allow tiny condensations to be resolved (e.g., Müller et al. 2003). When not having full resolution, this mainly leads to underestimation of the coronal density (e.g., Bradshaw & Cargill 2013). This can be easily understood, because at low resolution the transition region gets smeared out and can radiate the downward conducted heat at a lower density, which also leads to lower density in the corona. In 3D models, such resolution is (currently) impossible, but there are ways to counteract these shortcomings (e.g., Lionello et al. 2009). Despite their limited resolution, 3D MHD models can capture the essential response of the upper atmosphere to heat input and are consistent with the RTV scaling relations (e.g., Bourdin et al. 2016). In our models we are mostly concerned with the change in Xray luminosity with magnetic activity, which depends on the change of density (or pressure) in the corona, but not on its absolute value. Indeed, we are likely to find that our numerical model underestimates the coronal density. Still, even if the absolute values of the density are not fully correct, we can still retrieve the proper scaling of density, and hence of the Xray emission, which is the aim of our study.
In order to avoid instabilities because of steep gradients in temperature and density, we include in Eq. (10) a (numerical) isotropic heat conduction term and a mass diffusion term (see Bingert 2009). As is commonly used in many codes for modeling the solar corona (see e.g., Gudiksen et al. 2011; Rempel 2017) the mass diffusion term is used to smooth out density fluctuations on the grid scale and has no effect on the coronal dynamics. We also include a shock viscosity term for numerical reasons.
3.2 Initial and boundary conditions
The simulations are driven by (horizontal) motions on the solar surface that drive the magnetic field anchored in the photosphere. For the spatial distribution of the vertical magnetic field we employ an observed solar magnetogram (see Fig. 1). This is a snapshot of the active region AR 11102 as observed with the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) on August 30, 2010. A detailed description of a similar active region and how it is implemented in the model is discussed in Warnecke & Bingert (2020). These kinds of magnetograms represent a typical situation for a solar active region and similar magnetograms have been used as input for datadriven simulations in earlier studies (see e.g., Gudiksen & Nordlund 2005a,b; Bingert & Peter 2011).
For the initial condition, we use a potentialfield extrapolation to fill the box with magnetic field. The initial temperature stratification follows a vertical profile mimicking the temperature increase into the solar corona. The initial density is calculated from hydrostatic equilibrium and the system is initially at rest with all velocity components set to zero. Both initial temperature and density stratifications are motivated by the model of Vernazza et al. (1981).
In the horizontal xy plane, all variables are periodic. At the bottom boundary, the temperature T and density ρ have fixed values whereas the horizontal velocities, u_{x}, u_{y} have zero vertical gradients. The vertical velocity u_{z} is set to satisfy the divergencefree condition ∇⋅u = 0. We also prescribe a photospheric velocity driver which generates random photospheric motions; these are timedependent and create a velocity distribution and spectrum similar to the one observed on the solar photosphere. As a result, we shuffle the footpoints of the magnetic loops mimicking the solar photospheric flows in a similar way to, e.g., Gudiksen & Nordlund (2002, 2005a,b). At the top boundary, all velocity components are zero. To prevent any heat flux going in or out of the computational domain, the gradients of temperature and density are set to zero. The magnetic field is potential both for the bottom and the top boundary.
Fig. 1 Initial vertical magnetic field at the bottom boundary of the simulation. This is based on a magnetogram from active region AR 11102 observed with HMI. See Sect. 3.2. 
4 Numerical experiments
4.1 Setup
The main idea of this work is to start with an active region hosting only a small amount of total (unsigned) magnetic flux and then run models with increasing magnetic flux. Our aim is to study how the increase of the photospheric magnetic field strength contributes to the heating, and thus the Xray emission of the solar or stellar corona.
Here we report the results from a set of five numerical experiments. The original total unsigned surface magnetic flux of AR11102 is around 7 × 10^{20} Mx, which is a typical value of a small solar active region. We increase the flux by multiplying the B_{z} component of the surface magnetic field by a constant ranging from 1 to 20 (see Table 1). The spatial structure of the magnetic field remains unchanged. In all cases, the total magnetic flux at the bottom boundary is zero, that is, the surface magnetogram is balanced.
The heating of the coronal loops originates from the dissipation of Poynting flux into heat. More specifically, the conversion of the photospheric magnetic energy to thermal energy is due to the dissipation of the currents created by the random photospheric motions. The stronger the magnetic field, the more Poynting flux reaches the corona, resulting in a higher temperature and Xray emission.
In our setup, the different total unsigned magnetic fluxes correspond to peak values of the magnetic field ranging from 1 to 20 kG inside the spots; hence our naming convention in Table 1. The value of 20 kG is very high for a solar active region. Recent observations of solar active regions measured maximum values of the magnetic field strength on the order of 8 kG (Castellanos Durán et al. 2020). However, this high value was observed in only very small area on the active region light bridge. Typically the peak magnetic field strengths on solar active regions are on the order of 2–3 kG. Therefore, of the numerical experiments listed in Table 1, the 5B run can be considered to represent a typical solar active region, with a typical value for the total unsigned flux. However, higher values of surface magnetic field could be a common feature of very active stars (Reiners et al. 2014), and so the runs 10B and 20B could be considered to be representative of more active stars.
The increase in the surface magnetic field is not the only way to increase the surface magnetic flux. Alternatively, we could increase the horizontal extent of an active region while keeping the peak magnetic field strength the same. This will be the main focus of a future study.
Summary of numerical experiments.
4.2 Synthesized emission: Xrays and EUV
Coronal loops are mostly observed in extreme ultraviolet (EUV) wavelengths and Xrays because these wavelengths give access to the 1 MK plasma in the corona, and in these wavelengths, the dilute corona is also visible in front of the solar (or stellar) disk, which is very bright in the visible range. Therefore, we synthesize the emission in one EUV and one Xray band. For the EUV band we choose the widely used 171 Å band as seen by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the Solar Dynamic Observatory (SDO; Pesnell et al. 2012). Throughout the present paper we only consider the 171 Å channel of the EUV spectra and not the whole EUV range. We use the 171 Å channel of AIA as a reference for plasma at 1 MK and how it compares with plasma at higher temperature emitting in Xrays. For the Xray emission, we use the Alpoly filter of the solar Xray telescope (XRT; Golub et al. 2007) on board the Hinode observatory (Kosugi et al. 2007).
As mentioned in Sect. 2, the optically thin radiative losses through lines and continua in the corona observed in a given wavelength band are given through (15)
where n_{e} is the electron density and R(T) the temperature response (or contribution) function. The response needs to be calculated for each instrument (or filter) using the effective area depending on wavelength and the spectral lines forming in the wavelength region covered by the instruments. To calculate R(T) for the 171 Å channel of AIA and the Alpoly filter of XRT, we use the routines in the CHIANTI data base v9 (Dere et al. 1997, 2019) as they are available in the SOLARSOFT package^{1}. The response functions R(T) for AIA 171 Å and XRT Alpoly are shown in Fig. 2. To calculate the synthetic images as they would be observed by AIA or XRT, ε from Eq. (15) has to be integrated through the computational domain along the chosen line of sight. The samples for both channels for the 10B and 20B runs are shown in Fig. 3. These are integrated along the y direction which would correspond to an observation near the limb.
Fig. 2 Temperature response function for the AIA instrument on board SDO and the XRT on board Hinode. The black line shows the 171 Å channel of the AIA, the red line shows the Alpoly filter of the XRT. For illustration purposes, the black dashed line indicates a powerlaw approximation to XRT at temperatures below 10^{7} K. See Sect. 4.2. 
4.3 Horizontal averages
The increase of the surface magnetic field results in an increase in temperature and density in the coronal part of the domain. We calculate the horizontal averages of temperature T, density ρ, and the vertical component of the Poynting flux S_{z} and then average these in time for an interval of 1 h. During this time interval the computation reached a relaxed state, that is, the respective (spatially averaged) quantities show only rather small changes around a mean value (see Sect. 4.4 and Fig. 4). The average vertical stratification of T, ρ, and S_{z} is shown in Fig. 5.
4.3.1 Average Poynting flux deposited in the corona
The photospheric horizontal motions lead to an upwarddirected flux of magnetic energy, the Poynting flux. Here we concentrate on its vertical component, (16)
which is shown in Fig. 5c for the different runs. In the main part of the computational domain, the u × B × B term dominates and (on average) is positive, that is, upwards directed. The first term including the current, j, is significant only near the bottom where boundary effects of the driving cause high currents. Energetically, this is not relevant, because there the density is high enough that the heating through the currents has virtually no effect.
As we increase the total unsigned magnetic flux from one experiment to the next, the energy stored in the corona increases. The magnetic energy in excess of that of a potential field will be (partly) dissipated and converted into heat. The higher amount of dissipated (free) magnetic energy in the runs with higher magnetic flux leads to higher coronal temperatures and density (see Figs. 5a,b). This is just as expected from the RTV scaling laws (Eqs. (1) and (2)) as we discuss later in Sect. 6.2 (see also Fig. 6).
The Poynting flux for the runs 1B and 2B (i.e., black and blue solid lines in Fig. 5c) coincide in the lower part of the atmosphere (below 6 Mm). We find this to be the case only for the specific time frame used for the time averaging (i.e., from 3.5 to 4.5 h). For the time period before 3.5 h the Poynting fluxes for the runs 1B and 2B differ by a factor of two, as expected. This peculiar behavior is only observed for the run 1B for which the magnetic field is too low to even produce a proper corona. We observe this particular behavior only at the lower atmosphere while for the coronal part shown in Fig. 4, we see a clear distinction between runs 1B and 2B. However, as we only consider values on the coronal partfor our results, this peculiar behavior will not have any effect on the results presented in the following section.
Furthermore, the Poynting flux of the 1B run at the base of the corona (i.e. at an average temperature of around 0.1 MK) barely reaches 50 W m^{−2}. Typical estimations based on observations suggest an energy requirement of around 100 W m^{−2} for the quiet Sun and 10^{4} W m^{−2} above activeregions (Withbroe & Noyes 1977). Therefore, we cannot expect this run with the lowest magnetic activity to produce a megakelvin corona.
4.3.2 Average temperature and density
All of our simulations selfconsistently form a hot upper atmosphere, where the temperature is about two orders of magnitude higher than at the surface (cf. Fig. 5a). A higher total unsigned flux in the photosphere (cased 1B through 20B) corresponds not only to higher Poynting fluxes, but also to higher temperatures and density. The values shown in Figs. 5a,b are averages only, so the peak values are significantly higher, up to 5 MK and more.
The experiment with the lowest magnetic activity (run 1B) fails to create a megakelvin hot corona, as expected. Still, we consider it in the analysis of the powerlaw relation in Sect. 5. The main focus of this work is to relate the coronal emission to the surface magnetic activity through a number of numerical experiments, and in this sense a model that is not active enough to produce a megakelvin corona also provides valuable insight.
Besides the increased temperature and density, the models with higher magnetic activity also have the transition region located at lower heights. From this, it is clear that the height where the average temperature reaches 0.1 MK is lower for the runs with more magnetic flux. The higher energy input leads to a higher heat flux back to the Sun. Because the radiation is most efficient at lower temperatures (at 0.1 MK and below), in equilibrium the transition region will be found at lower temperatures and thus higher densities where it can radiate the energy. Consequently, the density (and the pressure) throughout the corona will be higher, as is indeed seen in our simulations.
The average density profile ρ displays similar (qualitative) behavior to the temperature. In the coronal part, the density is high for the runs with higher magnetic flux (Fig. 5b). Following a steep drop over many orders of magnitude in the low atmosphere, the density remains almost constant in the coronal part. This is simply because of the large barometric (pressure) scale height at high temperatures. At 1 MK, this scale height is about 50 Mm and thus comparable to the vertical extent of our computational domain; hence the horizontally averaged pressure and density are roughly constant in the coronal part of our box.
Fig. 3 Side view of the computational domain showing coronal emission. Left two panels: emission as it would be seen by the 171 Å channel of AIA starting from around 1 MK. Two righthand panels: Xray emission as seen by XRT sampling higher temperatures. Here we show snapshots of the two more active models, runs 10B and 20B. The emission here is integrated along the y direction which corresponds to an observation near the limb (of the Sun or a star). The snapshot is taken at t = 230 min; i.e., in the relaxed state. See Sect. 4.2. 
Fig. 4 Total coronal heating H_{tot} as a function of time. The vertical dashed lines indicate the time span used for the temporal averages. The colors represent the different runs as indicated in the legend (cf. Table 1). See Sect. 4.4. 
4.4 Temporal evolution
The quantities we consider in our model vary significantly with time, especially during the early phase of the simulations. To investigate the average behavior of our model, we have to consider a timeframe of the numerical model where the system reaches a relaxed (or evolved) state. During that state, the quantities show (comparably small) variations around an average value.
To illustrate this, we first consider the total heating in the coronal part of the computational domain. We define the coronal part as the volume above the height where the horizontally averaged temperature is 10^{5} K. Because the temperature gradient in the transition region around 10^{5} K is rather steep, the exact choice of this temperature is inconsequential. We therefore define the total coronal heating H_{tot} as the volume integral over this coronal part, here symbolized by the subscript ‘cor’, (17)
The temporal variation of the heating is shown in Fig. 4 for the five models with different magnetic activity. We see a clear ordering of the heating with the surface magnetic flux (increasing for run 1B through 20B), which we discuss in more detail in Fig. 7 and Sect. 6.1. In terms of the temporal evolution, we find that the heat input reaches a relaxed state rather quickly, probably within less than an hour. This is expected, because the stresses applied on the magnetic field in the photosphere will propagate with the Alfvén speed. The corresponding Alfvén crossing time for perturbations to cross the whole box is on the order of minutes.
The situation for the relaxation time is different when considering the coronal emission in Xrays and the 171 Å channel of AIA. Here it turns out that we have to wait for about three hours before the models reach a relaxed state. Mainly, this is because of the radiative cooling time under typical coronal conditions which is on the order of 1 h (Aschwanden 2005). Therefore, we examine the temporal evolution of the coronal radiation in some more detail below.
We now turn to the variability of the Xray and EUV emission represented by the 171 Å channel of AIA (see Fig. 8). Because the lower cool part of the atmosphere does not produce any significant amount of Xrays or EUV, we simply integrate the coronal emission over the whole computational domain. This is equivalent to the luminosity originating from the domain, L_{X} and L_{171}. Because we use the temperature response functions for XRT and AIA (Sect. 4.2, Fig. 2), we get the counts per second as expected for the respective instrument from the whole loop in the box (cf. Fig. 3). For the comparison between the different model runs, it is important that we use the same scale for the different models. Based on this, we see a clear scaling of the coronal emission with magnetic activity, in a similar way to what we see for the heat input. We discuss this in more detail in Fig. 9 and Sect. 6.3.
Overall, the coronal emission shows an initial drop on the timescale of almost an hour, in particular for EUV channel of AIA 171 Å in the runs with low magnetic activity. This is because the atmosphere of the initial condition is rather hot, and the plasma is cooling down until it reaches a new equilibrium after a few coronal radiative cooling times (Fig. 8). Here we see that all model runs reach a relaxed state after about 2.5–3 h.
The relative variability for the 1B and 2B runs is larger than for the runs with higher magnetic activity. Still, the absolute variability seen in the more active runs is much higher (the logarithmic plot in Fig. 8b shows the relative variation). Also, when following run 1B further in time, the emission in the 171 Å channel is not dropping further. The situation is different when considering other channels like the Xray regime which considers hotter plasma. There, the variability is relatively low even for the low magnetic activity runs.
Thus, with some safety margin we can assume that from 3.5 to 4.5 h the models have reached a relaxed state (see vertical dashed lines in Figs. 4 and 8). All the time averages discussed in our study are taken over this time frame.
Fig. 5 Horizontal averaged quantities as a function of height. We show temperature T (panel a), density ρ (panel b), and the vertical component of the Poynting flux S_{z} (panel c). The colors represent the different runs as indicated in the legend (cf. Table 1). The quantities are averaged horizontally for each snapshot and then in time for 1 h (between the timepoint at 3.5 and that at 4.5 h, as indicatedin Fig. 4 by the vertical dashed lines. For the Poynting flux we omitted the first three grid points that show boundary effects. See Sect. 4.3. 
5 Scaling relations in numerical experiments
To characterize the model runs with different magnetic activity, that is, unsigned magnetic surface flux, we investigate the scaling relations in the form of power laws between different parameters. In Sect. 4 we show that the enhancement of the photospheric magnetic flux leads to a substantial increase in temperature, density, Poynting flux, and coronal emission. In this section, we discuss the powerlaw relations of various quantities averaged in space and time. We discuss these results in Sect. 6 including the study of Zhuleku et al. (2020).
We first concentrate on the relation between the vertical component of the average Poynting flux ⟨ S_{z} ⟩ and both the unsigned surface magnetic flux Φ and the averaged total coronal heating ⟨H_{tot}⟩ (see Fig. 7). Each cross in the two figures represents the average value of the respective quantities in each individual numerical model (cf. Table 1). Here we take the horizontal average of the Poynting flux ⟨ S_{z} ⟩ at the height where the horizontally averaged temperature is 10^{5} K, that is, the base of the corona and average it in time from 3.5 to 4.5 h (see Sect. 4.4). This represents the energy flux (per unit area) into the corona. The total averaged coronal volumetric heating ⟨ H_{tot}⟩ is calculated according toEq. (17) and then averaged between times 3.5 and 4.5 h as discussed in Sect. 4.4. The unsigned surface flux Φ is the integral over the bottom boundary, that is, the stellar surface which is constant in time. The bars in Fig. 7 indicate the standard deviation of the respective quantity in time. For relations displayed in Fig. 7 we perform powerlaw fits (indicated by the red line) that result in (18) (19)
Here Γ corresponds to 1∕γ from the analytical model in Eq. (6) and Zhuleku et al. (2020). The two scalings Eqs. (18) and (19) imply that the heating increases roughly quadratically with magnetic flux, ⟨ H_{tot}⟩∝ Φ^{1.94}.
As a next step, we relate the coronal temperature T and density ρ to the total coronal heating ⟨H_{tot}⟩. Here we test to what extent the coronal temperature and density in our numerical model deviate from the analytic RTV scaling laws as given in Eqs. (1) and (2). To this end, we calculate the average temperature ⟨ T⟩ and density ⟨ ρ⟩ in the corona in a height range from z = 10 Mm to 20 Mm for each of the models. We choose this particular height range because this is where the bright structures appear (see Fig. 3). If we were to also include higher regions of the box, the averages would no longer represent the visible parts of the corona. In addition, we average ⟨ T⟩ and ⟨ ρ⟩ in time from 3.5 to 4.5 h as discussed in Sect. 4.4. We show the corresponding plots including the powerlaw fits of ⟨ T⟩ and ⟨ ρ⟩ as a function of the total coronal heating H_{tot} in Fig. 6. The powerlaw fits yield (20) (21)
Finally, we address the relation of the Xray emission to the heating rate and the unsigned surface magnetic flux. As the photospheric magnetic flux is larger for models with higher magnetic activity, the total dissipated energy is also larger. Because the Xray emission is expected to increase with the heating rate, it should also be larger for higher unsigned magnetic flux. For the analysis, we consider the averaged emission, ⟨ L_{X} ⟩, that is integrated over the whole computational domain and averaged in the same way as the other quantities from 3.5 to 4.5 h. The corresponding relations are plotted in Figs. 9 and 10a, and the powerlaw fits give (22) (23)
We apply the same integration and time averaging to the EUV emission as seen by AIA in its 171 Å channel. Its relation to the heat input is displayed in Fig. 10b and the powerlaw fit reveals an almost linear relation, (24)
Fig. 6 Scaling of average temperature ⟨T⟩ and density ⟨ρ⟩ with averagecoronal heat input ⟨H_{tot}⟩. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux listed in Table 1. The bars represent the standard deviation of the spatial averages of T and ρ in time. The red lines are powerlaw fits to the data. The blue lines indicate what is expected from the RTV scaling laws. See Sects. 5 and 6.2. 
Fig. 7 Scaling of average Poynting flux ⟨S_{z}⟩ with unsigned surface flux Φ and average coronal heating ⟨H_{tot}⟩. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux as listed in Table 1. As a reference for the magnetic flux, Φ_{ref}, we choose the magnetic flux of the least active setup, run 1B (cf. Table 1). The bars represent the standard deviation of S_{z} in time. The red lines are powerlaw fits to the data. See Sects. 5 and 6.1. 
6 Discussion
6.1 Energy input into the corona
The solar coronal heating problem and the underlying physical mechanism has been extensively discussed over the last 70 yr. Two of the main mechanisms that are being considered are the Alfvénwave model (e.g., van Ballegooijen et al. 2011) and the fieldline braiding or nanoflare model (Parker 1972, 1983). In these two processes, the Poynting flux, and hence the heat input, have a different dependence on the magnetic field. A simple estimate for these dependencies was given by Fisher et al. (1998), and here we follow their arguments. For an Alfvén wave with constant amplitude, the Poynting flux will be proportional to the propagation speed, that is, the Alfvén velocity, and hence to the magnetic field B, S_{z} ∝ B. For the fieldline braiding, the Poynting flux will be set by the driving motions with speed u and the magnetic field, S ∝u × B × B, and hence the scaling will be S_{z} ∝ B^{2} (assuming that the driving motions do not change with B).
In our numerical experiments we find a powerlaw scaling between Poynting flux and (unsigned) magnetic flux with a powerlaw index of about 1.7 ± 0.4; see Eq. (18) and Fig. 7a. Within the uncertainties this is consistent with the above (analytical) estimate of 2. This is not surprising, because in our model we drive the coronal magnetic field through footpoint motions consistent with the fieldline braiding scenario. Still, it is reassuring to recover this scaling by our numerical model.
Observationally, it is clear that Alfvén waves are present in the corona (e.g., Tomczyk et al. 2007). Still, it remains unclear as to whether or not the energy they carry is sufficient to energize the corona. Waves might play a role in the quiet Sun corona, but they appear to be unable to heat active regions (McIntosh et al. 2011). Our particular model can contribute little to thisdiscussion, because we do not fully resolve Alfvén waves, which is because of the comparably large dissipation.
Finally, we expect that the total energy dissipated in the coronal volume matches the Poynting flux at the base of the corona. In this case, the Poynting flux should scale linearly with the total amount of energy dissipated in the corona. In our numerical experiments, we find a powerlaw relation with a powerlaw index of about 0.9 ± 0.2, which is consistent with linear; see Eq. (19) and Fig. 7b.
Fig. 8 Temporal evolution of the coronal emission from integration over the whole computational domain. Panel a: Xray emission as seen by XRT in the Alpoly filter. Panel b: EUV emission as it would be seen by AIA in the 171 Å channel. The vertical dashed lines indicate the timespan used for time averaging. The colors represent the different runs as indicated in the legend (cf. Table 1). See Sect. 4.4. 
Fig. 9 Scaling of average Xray emission ⟨L_{X}⟩ with unsigned surface flux Φ. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux listed in Table 1. As a reference for the magnetic flux, Φ_{ref}, we choose the magnetic flux of the least active setup, run 1B (cf. Table 1). The bars represent the standard deviation of L_{X} in time. The red lines are powerlaw fits to the data. The blue lines show the analytic powerlaw relations based on the RTV scalinglaws; see Sects. 5, 6.3 and 6.4. 
6.2 RTV scaling laws compared to numerical experiments
One obvious check for the numerical experiments is to what extent the average quantities will follow the RTV scaling laws of Eqs. (1) and (2). Because of the spatial and temporal variability we cannot expect a perfect match, but the average quantities should roughly follow these scalings, as was found in an earlier model for one single (small) solar active region (Bourdin et al. 2016).
The original RTV scalings in Eqs. (1) and (2) include a dependence on the loop length L. However, in our numerical models the physical size of the computational box is kept constant. Therefore, the coronal loops can be considered to have similar lengths. Consequently, in this study we only have to consider the dependence of temperature and density on (total) heating in the coronal volume, H_{tot}.
For a comparison to the RTV scalings one should not only compare the powerlaw indices of the scaling, but also the absolute values of the temperature T and density ρ as given in the original paper by Rosner et al. (1978). To calculate the predictions from the RTV scalings, that is, T_{RTV} and ρ_{RTV}, we use a loop length of L = 30 Mm, which is similar to the average loop length we find in the emission patterns synthesized from the model (see Fig. 3). For the volumetric heating rate we use the (total) heating in the coronal volume, ⟨ H_{tot}⟩, divided by the coronal volume, which gives the (average) volumetric heating. The resulting variation of T_{RTV} and ρ_{RTV} with ⟨ H_{tot}⟩ is shown in Fig. 6 (blue lines). The power law indices of 0.29 and 0.57 for these correspond to 2/7 and 4/7 in Eqs. (1) and (2).
The powerlaw relation of the average temperatures ⟨T⟩ and density ⟨ρ⟩ in the numerical models with ⟨H_{tot}⟩ are close to what is expected from RTV. However, both the temperature and the density are significantly lower, by about a factor two and three; see Fig. 6 and Eqs. (20) and (21).
The reason for this underestimation is mainly due to the averaging process of the temperature and density. The bright loops are hotter and denser than the ambient corona, meaning that the averages underestimate the values corresponding to the RTV scaling relations. Still, we conclude that the numerical models and the analytical scaling relations are consistent (as an orderofmagnitude estimation).
6.3 Relation of Xray emission to surface magnetic flux
The central result of our study is the powerlaw relation between the synthetic Xray emission L_{X} and the unsigned surface magnetic flux Φ. We show this relation L_{X} ∝ Φ^{m} with a powerlaw index m of about 3.4 ± 0.3 in Fig. 9; see also Eq. (22).
This L_{X} versus Φ^{m} relationship has been extensively discussed in numerous observational studies of the Sun and of stars of different spectral types (see e.g., Fisher et al. 1998; Pevtsov et al. 2003; Vidotto et al. 2014). However, the physical reasons behind the observed powerlaw relation remain unclear. Observations of various solar features and solarlike stars suggest a relation close to but slightly steeper than linear, i.e., m ≳ 1 (e.g., Fisher et al. 1998; Pevtsov et al. 2003).
We can estimate this powerlaw index m from observations by comparing the (average) Xray emission and magnetic field of the Sun to those of other stars. The Sun has an average magnetic field of roughly 10 G (Cranmer 2017), active stars show average fields of about 1 kG (Saar 1996; Saar & Brandenburg 2001). This reveals a difference of a factor of 100 in magnetic field strength. At the same time the Xray emission increases from the Sun to active stars by a factor of 1000 (Xrayrotation activity relation) (see e.g., Pizzolato et al. 2003). From this one can conclude that the L_{X} versus Φ increases as a power law with an of index m = 1.5. However, more recent studies find the power law to be close to quadratic (m = 1.8; Vidotto et al. 2014) or even steeper (m = 2.68; Kochukhov et al. 2020). While Kochukhov et al. (2020) considered stars similar to the Sun in terms of spectral type, all these were considerablymore active than the Sun in terms of Xray luminosity. The large sample of Vidotto et al. (2014) also mostly contained stars significantly more active than the Sun.
An almost linear relation of L_{X} versus Φ could be simply understood by increasing the number of active regions on a (solarlike) star. If the filling factor of active regions is low and one simply doubles the number of active regions, both the total unsigned magnetic flux on the stellar surface and the Xray emission would double; hence, the linear relation between L_{X} and Φ. This, in principle, could explain the findings of Fisher et al. (1998) and Pevtsov et al. (2003).
Active stars show Xray luminosities (compared to their bolometric luminosity) that can be three or more orders of magnitude larger than that of the Sun (e.g., Vidotto et al. 2014). In this case there would simply not be enough space on the star to cover it with enough(solarlike) active regions. This means that the Xray emission per active region has to increase. This is exactly what we find in our model when increasing the magnetic flux of an active region while keeping its size (i.e. area) the same. This leads to a steep increase in Xray luminosity with surface magnetic flux, an increase that is even steeper than observed: the powerlaw index m for L_{X} ∝ Φ^{m} is about 3.4 in our model, while the largest value found in observations is below 2.7 (Kochukhov et al. 2020). This overestimation of the powerlaw index by our model indicates that on real stars we might find a mixture of an increase in the numbers of active regions together with an increase in the peak magnetic field strength (or magnetic flux per active region) for more active stars.
Furthermore, the powerlaw index m might be overestimated by our model. Strong surface magnetic fields could quench the horizontal motions and consequently reduce the resulting Poynting flux (e.g., Gudiksen & Nordlund 2005a). This is not included in our model. In Appendix A, we discuss the impact that this quenching could have on our estimations of the powerlaw relations. We find that a realistic quenching of the horizontal velocities by about a factor of three would reduce the index m of the powerlaw relation between Xray emission and magnetic field, L_{X} ∝ Φ^{m} by about 25% to around 2.6. This would bring the value for m derived from our model closer to the index derived by recent observations.
Fig. 10 Scaling of average Xray and 171 Å emission as observed by AIA, ⟨L_{X}⟩, and ⟨ L_{171}⟩ with averagecoronal heating ⟨H_{tot}⟩ for model runs listed in Table 1. The bars represent the standard deviation of L_{X} and L_{171} in time. The red lines are powerlaw fits to the data. The blue line in panel a indicates what is expected from the RTV scaling laws. See Sects. 5 and 6.5. 
6.4 Analytical model for scaling of Xray emission
We now compare the scaling of L_{X} ∝ Φ^{m} to some basic analytical considerations. In an earlier study we derived an analytical scaling of the Xray emission with surface magnetic flux that is based on the RTV scaling relations (Zhuleku et al. 2020). As discussed in the summary of that model in Sect. 2, there we used a scaling of active region size with (unsigned) magnetic flux based on solar observations, see Eq. (7). In contrast, in our numerical model we keep the size or area covered by the active region constant; hence, here the powerlaw index in Eq. (7) is δ = 0. This simplifies our analytical scaling from Eq. (3) to (25)
According to Eq. (4), α is the powerlaw index relating the temperature response (or contribution) function to temperature. Here we approximate this for Xrays as seen by the XRT on Hinode (see Fig. 2), where a powerlaw fit yields (26)
We provide a more extensive discussion on the temperature responses for different instruments in Zhuleku et al. (2020).
We take the other two parameters β and γ in Eq. (25) from the relations of the Poynting flux to the unsigned surface magnetic flux and the total heating in Eqs. (18) and (19), with γ = 1∕Γ in Eq. (19). This results in the analytic scaling (based on RTV) of (27)
This is overplotted on the results from the numerical experiments in Fig. 9 as a blue line. We conclude that this powerlaw index is in good agreement with the value we obtained from the numerical models, see Eq. (22) and Fig. 9.
Just as for the comparison to the RTV scalings in Sect. 6.2, also here it is not sufficient to find a match of the powerlaw index; the absolute values of the derived Xray emission also have to be of the same order for a good match. To get the constant of proportionality in Eq. (27), we have to assign a volume of the emitting structure described in the analytical scaling. For the comparison to our numerical model we therefore assign the volume of the loop(s) dominating the coronal emission. Using the loops in Fig. 3 as a guideline, these have a length and a radius of about L ≈ 30 Mm and r ≈ 2.5 Mm. With the volume of V = πR^{2}L we then find the total Xray radiation from the loops based on the (RTV) scaling relations of . Within an order of magnitude these match what we find in the numerical study. The deviation is mainly because of the factor of about three difference in density (cf. Sect. 6.2) which enters the emission quadratically.
6.5 Xray and EUV emission versus coronal energy input
Our current understanding of solar physics would lead us to believe that a higher energy input into the corona should result in a higher Xray emission. Indeed, this is what we find in Fig. 10a. However, we find an almost quadratic relationship, ⟨ L_{X} ⟩∝⟨H_{tot}^{⟩1.8}, cf. Eq. (24). This nonlinearity is counterintuitive: Assuming that the plasma in our numerical model has reached a relaxed state, then whatever energy reaches the corona should be radiated away as Xray emission (i.e., naively, we would expect a linear dependence).
This (roughly) quadratic relation can be understood when going back to the RTV scalings. The Xray emission is given by L_{X} ∝ n^{2} R_{X}. Combining this with the fit to the Xray response function in Eq. (26) and the RTV scaling relations in Eqs. (1) and (2) provides the analytical scaling for the Xray emission with (average) coronal heat input, (28)
This is shown in comparison to the results from the numerical model in Fig. 10a as a blue line. As before, we assume a constant loop length, because the active region size is the same for all numerical experiments.
This consideration provides an understanding of the nonlinear relation between Xray output and heat input. In the corona the energy input will be only partly balanced by the Xray output (heat conduction, wind outflow, and radiation at other wavelengths also play their respective roles). Because the Xray output shows a nonlinear temperature dependence, naturally it will also be connected in a nonlinear fashion to the heat input. The simple RTV scaling relations, together with the temperature response of Xray emission give rise to the roughly quadratic dependence of Xrays on heat input.
The situation changes when investigating another wavelength region, which essentially implies probing a different temperature range. For this we use the channel as seen by AIA at 171 Å. The emission from the 171 Å channel of AIA in this band mostly originates from temperatures around 1 MK (cf. Fig. 2). If we apply the same analysis as for the Xray emission, the numerical experiments show an almost linear scaling with heat input, ⟨ L_{171}⟩∝⟨H_{tot}^{⟩1.1}; see Fig. 10b. The reason this relation is less steep than for the Xrays is because of the different temperature response in the EUV.
The response of the 171 Å channel of AIA cannot be simply approximated by a power law (see Fig. 2). However, over the temperature range of most our models, as a very rough zerothorder approximation we could assume the response in the 171 Å channel to be constant. This would be consistent with the powerlaw index for the response function, that is, α in Eq. (4), being equal to zero. Following the above procedure for the Xrays, in analogy to Eq. (28) we would find for the 171 Å channel emission. This result compares well with the numerical experiments (cf. Fig. 10b), but certainly it has to be taken cum grano salis, because of our very rough assumption of α = 0. This is why we do not overplot an analytical approximation in Fig. 10b as we did for the other scaling relations. Still, this consideration underlines the importance of the temperature response in establishing a relationship between coronal emission and heat input.
Fig. 11 Scaling of averages of Xray emission ⟨L_{X}⟩ with 171 Å emission as observed by AIA ⟨L_{171}⟩ for model runs listed in Table 1. The bars represent the standard deviation of L_{X} in time. The red line is a powerlaw fit to the data; see Sect. 6.6. 
6.6 Relating Xrays to EUV emission
Finally, we briefly investigate the scaling relation between the Xray and the emission of the 171 Å channel of AIA. Solar and stellar studies haveaddressed the relation between the radiative fluxes in different wavelength bands, also called fluxflux relations. One prominent example is the relation between the coronal Xray radiative flux and the radiance in the chromospheric Ca II line with (e.g., Schrijver 1983). Powerlaw relations have also been derived with other chromospheric and transition region lines, such as Mg II, Si II, C II, and C IV (Schrijver 1987). The higher the temperature in the source region of the emission, the steeper the powerlaw relation of the respective line to the surface magnetic field (e.g., Testa et al. 2015). This implies that the relation of ⟨ L_{X} ⟩ versus ⟨ L_{171}⟩ is less steep than the relation of Xrays to Ca II.
In the numerical models we find that the relation between Xrays and the 1 MK EUV emission of the 171 Å channel of AIA follows a power law with about ⟨L_{X}⟩∝⟨L_{171}^{⟩1.3}; see Fig. 11. This is also consistent with other numerical studies (see e.g., Warnecke & Peter 2019b, their Fig. 10b). Therefore, as expected, it is less steep than Xrays to Ca II with a powerlaw index of about 1.7.
Our numerical models are not realistic enough when it comes to the cooler parts of the atmosphere, in the transition region around 10^{5} K and in particular in the cooler chromosphere. Still we see a consistent trend for the fluxflux relations, and their investigation is left for future more realistic models.
7 Conclusions
We performed a series of 3D MHD numerical simulations of active regions. With these, we address the question of how the (Xray) emission from stellar coronae scales with the surface magnetic flux. Observationally this powerlaw scaling L_{X} ∝ Φ^{m} is well established with the most recent studies giving powerlaw indices m in the range of below 2 to almost 3 (Vidotto et al. 2014; Kochukhov et al. 2020). So far, the physical basis for this relation is poorly understood in terms of (numerical) models.
For our model series, we assumed that the area covered by the active region remains the same, and we changed the peak (or average) vertical field strength B to change the (unsigned) magnetic flux at the surface by a factor of 20. This resulted in a change of the Xray emission by more than four orders of magnitude. The scaling we found in our numerical experiments is a power law with an index m ≈ 3.4, which is slightly steeper than found in observations. As we discussed in Sect. 6.3, this difference can be understood if on the real star there are not only active regions with higher magnetic flux (but the same area) but also a larger filling factor of active regions, or in other words a larger number of active regions.
The results of our numerical experiments are consistent with an analytical scaling model (Zhuleku et al. 2020). This is based on the RTV scaling laws connecting the temperature and density to the heating rate and size of the coronal structure and the temperature response of the wavelength band in which the observations are performed. With this, we have a clear understanding of how and why the radiative Xray output changes in a nonlinear fashion in response to the heat input, and hence the surface magnetic flux.
In our numerical model we choose the specific approach of increasing the surface magnetic flux by increasing the field strength. This is motivated by the very high field strengths seen on active M dwarfs of up to 8 kG on average (e.g., Reiners 2012). Such a large average suggests that the highest field strengths on these stars could be even higher, perhaps consistent with our model. Still, to investigate further possibilities we will study the response of the coronal emission to an increase in the magnetic flux by increasing the area of the active region(s) in a future project.
Our models for individual active regions cannot be expected to fully account for all aspects of the scaling of coronal emission with (unsigned) surface magnetic flux. Still, our approach indicates a way in which to understand the excessive increase in the observed Xray emission by four orders of magnitude from solartype activity to fast rotating active stars (by e.g., Pizzolato et al. 2003; Vidotto et al. 2014).
Acknowledgements
The simulations have been carried out at GWDG in Göttingen and the Max Planck Computing and Data Facility (MPCDF) in Garching. This work was supported by the International Max Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. J.W. acknowledges funding by the MaxPlanckPrinceton Center for Fusion and Astro Plasma Physics. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 818665 “UniSDyn”).
Appendix A Quenching of convective motions
Fig. A.1 Effect of quenching of horizontal motions by strong magnetic field. We show the relation of the quenched Poynting flux as defined in Eq. (A.5) to the unsigned surface magnetic flux Φ. Each data point represents an average in time and space for each run. The values are normalized to model run 1B represented here by the subscript ‘ref’. The black line shows the powerlaw fit to the case where no quenching is included (i.e. q_{sat} = 1). The red and blue lines show powerlaw fits for q_{sat} = 2∕3 and 1/3. The corresponding powerlaw indices β for are given with the legend. 
The current numerical experiments presented in Sect. 4 might overestimate the amount of Poynting flux injected from the bottom boundary. This is because, at the surface, (strong) magnetic fields would quench the convective motions. To properly account for this effect, the nearsurface magnetoconvection would have to be included selfconsistently into the model. Still, we can estimate what effect this quenching would have on the powerlaw indices we derive.
The horizontal velocities at the surface, u_{h}, that drive the magnetic field will be reduced in the strongfield regions according to (A.1)
The quenching function q can be expressed as a function of the vertical magnetic field at the surface, B_{z,0}, (A.2)
Essentially, if the magnetic field is well below the equipartition field strength, B_{eq}, there is no quenching, while for large B_{z,0} the surface velocities are reduced by the factor q_{sat}. For our experiments, we assume values of q_{sat} of 1, 2/3, and 1/3.
For a smooth transition of the quenching q from 1 to q_{sat} we apply a cubic step function that depends on the logarithm of the magnetic field. We assume that the inflection point of the function q is at the equipartition field strength, B_{eq}. In our models this is (on average) 2.3 kG. The transition from q = 1 to q_{sat} happens over a range of 0.8 in log _{10}B_{z,0}. This quenching as a function of magnetic field (for the case q_{sat} = 1∕3) is very close to the quenching used by Gudiksen & Nordlund (2005b) as defined in their Eq. (12).
In general, the Poynting flux is defined as (A.3)
The first term j × B is not affected by the quenching of horizontal motions. The second term can be further decomposed into two terms, (A.4)
where we now use the horizontal velocities quenched by the magnetic field, . The subscript h refers to horizontal component combining the x and y direction. From the two terms in Eq. (A.4) the first term is less relevant in our context, because the vertical velocities are smaller than the horizontal ones at the surface, and because the horizontal magnetic field is smaller than the vertical one. Only the second term, (A.5)
will be affected by quenching of the horizontal velocities and is of foremost interest in the discussion in this appendix.
To test the effect of the quenching of the horizontal motions on the Poynting flux we therefore evaluate from Eq. (A.5). To this end, we calculate the quenching term q as defined in Eq. (A.2) at each location at the bottom boundary, i.e., the surface, calculate the quenched velocity according to Eq. (A.1), and then derive the vertical component of the Poynting flux due to the quenched horizontal motions, , at the bottom boundary of our model runs. As done for the other (average) quantities in our study, the quenched Poynting flux at the bottom is (horizontally) averaged in space and in time during the relaxed state of the simulations. For better comparison between the different models, we normalize the (average) magnetic flux and Poynting flux at the bottom by the respective values from model run 1B.
To quantify the effect of the quenching on the scaling between magnetic flux Φ and Poynting flux at the surface , we conduct experiments for three cases with q_{sat} = 1, 2/3, and 1/3. These experiments are illustrated in Fig. A.1. For the case with no quenching (q_{sat}=1), we find a powerlaw scaling with a powerlaw index β = 1.5. This is slightly smaller than β in the full model in Fig. 7 and Eq. (18), but still within the range of uncertainties given there. The small difference is because of neglecting the other terms of the Poynting flux in Eq. (A.3) and because in the main text we consider the Poynting flux at the base of the corona (see Sect. 6.1) and not at the bottom boundary as done here. More importantly, the powerlaw relation gets more shallow for stronger quenching, with the powerlaw index β going down to about 1.4 and 1.15 for values of q_{sat} of 2/3 and 1/3. Consequently, for a more realistic setup with a quenching similar to Gudiksen & Nordlund (2005b), i.e., q_{sat} = 1∕3, we would find a reduction of the powerlaw index β by about 25% compared to the case of no quenching.
The reduction of the powerlaw index β by the quenching of horizontal motions will also affect how the Xray luminosity L_{X} scales with the surface magnetic flux Φ. In the framework of the analytical model in Sect. 6.4, the powerlaw index m of L_{X} ∝ Φ^{m} is expressed in Eq. (25) as a function of three parameters, α, β, γ, of which only the parameter β will be affected by the quenching. α and γ will retain the values discussed in Sect. 6.4. Following Eq. (25), m is linear in β. Hence its reduction through the quenching for a more realistic setup would be the same as for β, that is, about 25%. Therefore we would expect the quenching to reduce m as found in our numerical experiments as displayed in Fig. 9 and discussed in Sect. 6.3 to about m = 2.6. This would be very close to the value of 2.68 found by Kochukhov et al. (2020) in their recent observational study, but still considerably higher than the traditional value that is assumed to be not considerably higher than unity (Fisher et al. 1998; Pevtsov et al. 2003).
References
 Aschwanden, M. J. 2005, Physics of the Solar Corona, An Introduction with Problems and Solutions, 2nd edn. (SpringerVerlag Berlin Heidelberg New York) [Google Scholar]
 Bingert, S. 2009, Heating of the Corona in a 3D MHD Forward Model Approach (Breisgau: University of Freiburg) [Google Scholar]
 Bingert, S., & Peter, H. 2011, A&A, 530, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bingert, S., & Peter, H. 2013, A&A, 550, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boris, J. P. 1970, in A Physically Motivated Solution of the Alfven Problem, NRL memorandum report [CrossRef] [Google Scholar]
 Bourdin, P.A., Bingert, S., & Peter, H. 2016, A&A, 589, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bradshaw, S. J., & Cargill, P. J. 2013, ApJ, 770, 12 [Google Scholar]
 Castellanos Durán, J. S., Lagg, A., Solanki, S. K., & van Noort, M. 2020, ApJ, 895, 129 [CrossRef] [Google Scholar]
 Chatterjee, P. 2020, Geophys. Astrophys. Fluid Dyn., 114, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, F., & Peter, H. 2015, A&A, 581, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chen, F., Peter, H., Bingert, S., & Cheung, M. C. M. 2014, A&A, 564, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cook, J. W., Cheng, C.C., Jacobs, V. L., & Antiochos, S. K. 1989, ApJ, 338, 1176 [Google Scholar]
 Cranmer, S. R. 2017, ApJ, 840, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Dahlburg, R. B., Einaudi, G., Taylor, B. D., et al. 2016, ApJ, 817, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Dahlburg, R. B., Einaudi, G., UgarteUrra, I., Rappazzo, A. F., & Velli, M. 2018, ApJ, 868, 116 [CrossRef] [Google Scholar]
 Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dere, K. P., Del Zanna, G., Young, P. R., Landi, E., & Sutherland, R. S. 2019, ApJS, 241, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, G. H., Longcope, D. W., Metcalf, T. R., & Pevtsov, A. A. 1998, ApJ, 508, 885 [NASA ADS] [CrossRef] [Google Scholar]
 Golub, L., Deluca, E., Austin, G., et al. 2007, Sol. Phys., 243, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Gombosi, T., Tóth, G., Zeeuw, D., et al. 2002, J. Comput. Phys., 177, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Gudiksen, B. V., & Nordlund, Å. 2002, ApJ, 572, L113 [Google Scholar]
 Gudiksen, B. V., & Nordlund, Å. 2005a, ApJ, 618, 1020 [NASA ADS] [CrossRef] [Google Scholar]
 Gudiksen, B. V., & Nordlund, Å. 2005b, ApJ, 618, 1031 [NASA ADS] [CrossRef] [Google Scholar]
 Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hansteen, V. 1993, ApJ, 402, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Hansteen, V., Guerreiro, N., De Pontieu, B., & Carlsson, M. 2015, ApJ, 811, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Kochukhov, O., Hackman, T., Lehtinen, J. J., & Wehrhahn, A. 2020, A&A, 635, A142 [CrossRef] [EDP Sciences] [Google Scholar]
 Kosugi, T., Matsuzaki, K., Sakao, T., et al. 2007, Sol. Phys., 243, 3 [Google Scholar]
 Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17 [Google Scholar]
 Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902 [Google Scholar]
 Magaudda, E., Stelzer, B., Covey, K. R., et al. 2020, A&A, 638, A20 [CrossRef] [EDP Sciences] [Google Scholar]
 Matsumoto, T. 2021, MNRAS, 500, 4779 [CrossRef] [Google Scholar]
 McIntosh, S. W., de Pontieu, B., Carlsson, M., et al. 2011, Nature, 475, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Meyer, J. P. 1985, ApJS, 57, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, D. A. N., Hansteen, V. H., & Peter, H. 2003, A&A, 411, 605 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murphy, R. J. 1985, PhD thesis, University of Maryland, USA [Google Scholar]
 Parker, E. N. 1972, ApJ, 174, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1983, ApJ, 264, 642 [Google Scholar]
 Pencil Code Collaboration (Brandenburg, A., et al.) 2021, J. Open Source Softw., 6, 2807 [Google Scholar]
 Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Sol. Phys., 275, 3 [Google Scholar]
 Peter, H., Gudiksen, B. V., & Nordlund, Å. 2004, ApJ, 617, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Pevtsov, A. A., Fisher, G. H., Acton, L. W., et al. 2003, ApJ, 598, 1387 [NASA ADS] [CrossRef] [Google Scholar]
 Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2008, ApJ, 677, 1348 [Google Scholar]
 Reiners, A. 2012, Liv. Rev. Sol. Phys., 9, 1 [Google Scholar]
 Reiners, A., Schüssler, M., & Passegger, V. M. 2014, ApJ, 794, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2017, ApJ, 834, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643 [Google Scholar]
 Saar, S., & Brandenburg, A. 2001, ASP Conf. Ser., 248, 231 [Google Scholar]
 Saar, S. H. 1996, in Steller Surface Structure, eds. K. G. Strassmeier, & J. L. Linsky (Dordrecht: Kluwer Academic Publishers), 176, 237 [NASA ADS] [Google Scholar]
 Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [Google Scholar]
 Schrijver, C. J. 1983, A&A, 127, 289 [NASA ADS] [Google Scholar]
 Schrijver, C. J. 1987, A&A, 172, 111 [NASA ADS] [Google Scholar]
 Spitzer, L. 1962, Physics of Fully Ionized Gases, 2nd edn. (New York: Interscience) [Google Scholar]
 Testa, P., Saar, S. H., & Drake, J. J. 2015, Phil. Trans. R. Soc. London, Ser. A, 373, 20140259 [NASA ADS] [CrossRef] [Google Scholar]
 Tomczyk, S., McIntosh, S. W., Keil, S. L., et al. 2007, Science, 317, 1192 [Google Scholar]
 van Ballegooijen, A. A., AsgariTarghi, M., Cranmer, S. R., & DeLuca, E. E. 2011, ApJ, 736, 3 [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [Google Scholar]
 Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., & Bingert, S. 2020, Geophys. Astrophys. Fluid Dyn., 114, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., & Peter, H. 2019a, A&A, 624, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., & Peter, H. 2019b, A&A, submitted [arXiv:1910.06896] [Google Scholar]
 Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, N. J., Drake, J. J., Mamajek, E. E., & Henry, G. W. 2011, ApJ, 743, 48 [Google Scholar]
 Zhuleku, J., Warnecke, J., & Peter, H. 2020, A&A, 640, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1 Initial vertical magnetic field at the bottom boundary of the simulation. This is based on a magnetogram from active region AR 11102 observed with HMI. See Sect. 3.2. 

In the text 
Fig. 2 Temperature response function for the AIA instrument on board SDO and the XRT on board Hinode. The black line shows the 171 Å channel of the AIA, the red line shows the Alpoly filter of the XRT. For illustration purposes, the black dashed line indicates a powerlaw approximation to XRT at temperatures below 10^{7} K. See Sect. 4.2. 

In the text 
Fig. 3 Side view of the computational domain showing coronal emission. Left two panels: emission as it would be seen by the 171 Å channel of AIA starting from around 1 MK. Two righthand panels: Xray emission as seen by XRT sampling higher temperatures. Here we show snapshots of the two more active models, runs 10B and 20B. The emission here is integrated along the y direction which corresponds to an observation near the limb (of the Sun or a star). The snapshot is taken at t = 230 min; i.e., in the relaxed state. See Sect. 4.2. 

In the text 
Fig. 4 Total coronal heating H_{tot} as a function of time. The vertical dashed lines indicate the time span used for the temporal averages. The colors represent the different runs as indicated in the legend (cf. Table 1). See Sect. 4.4. 

In the text 
Fig. 5 Horizontal averaged quantities as a function of height. We show temperature T (panel a), density ρ (panel b), and the vertical component of the Poynting flux S_{z} (panel c). The colors represent the different runs as indicated in the legend (cf. Table 1). The quantities are averaged horizontally for each snapshot and then in time for 1 h (between the timepoint at 3.5 and that at 4.5 h, as indicatedin Fig. 4 by the vertical dashed lines. For the Poynting flux we omitted the first three grid points that show boundary effects. See Sect. 4.3. 

In the text 
Fig. 6 Scaling of average temperature ⟨T⟩ and density ⟨ρ⟩ with averagecoronal heat input ⟨H_{tot}⟩. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux listed in Table 1. The bars represent the standard deviation of the spatial averages of T and ρ in time. The red lines are powerlaw fits to the data. The blue lines indicate what is expected from the RTV scaling laws. See Sects. 5 and 6.2. 

In the text 
Fig. 7 Scaling of average Poynting flux ⟨S_{z}⟩ with unsigned surface flux Φ and average coronal heating ⟨H_{tot}⟩. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux as listed in Table 1. As a reference for the magnetic flux, Φ_{ref}, we choose the magnetic flux of the least active setup, run 1B (cf. Table 1). The bars represent the standard deviation of S_{z} in time. The red lines are powerlaw fits to the data. See Sects. 5 and 6.1. 

In the text 
Fig. 8 Temporal evolution of the coronal emission from integration over the whole computational domain. Panel a: Xray emission as seen by XRT in the Alpoly filter. Panel b: EUV emission as it would be seen by AIA in the 171 Å channel. The vertical dashed lines indicate the timespan used for time averaging. The colors represent the different runs as indicated in the legend (cf. Table 1). See Sect. 4.4. 

In the text 
Fig. 9 Scaling of average Xray emission ⟨L_{X}⟩ with unsigned surface flux Φ. Each data point represents an average for one of the model runs with different unsigned surface magnetic flux listed in Table 1. As a reference for the magnetic flux, Φ_{ref}, we choose the magnetic flux of the least active setup, run 1B (cf. Table 1). The bars represent the standard deviation of L_{X} in time. The red lines are powerlaw fits to the data. The blue lines show the analytic powerlaw relations based on the RTV scalinglaws; see Sects. 5, 6.3 and 6.4. 

In the text 
Fig. 10 Scaling of average Xray and 171 Å emission as observed by AIA, ⟨L_{X}⟩, and ⟨ L_{171}⟩ with averagecoronal heating ⟨H_{tot}⟩ for model runs listed in Table 1. The bars represent the standard deviation of L_{X} and L_{171} in time. The red lines are powerlaw fits to the data. The blue line in panel a indicates what is expected from the RTV scaling laws. See Sects. 5 and 6.5. 

In the text 
Fig. 11 Scaling of averages of Xray emission ⟨L_{X}⟩ with 171 Å emission as observed by AIA ⟨L_{171}⟩ for model runs listed in Table 1. The bars represent the standard deviation of L_{X} in time. The red line is a powerlaw fit to the data; see Sect. 6.6. 

In the text 
Fig. A.1 Effect of quenching of horizontal motions by strong magnetic field. We show the relation of the quenched Poynting flux as defined in Eq. (A.5) to the unsigned surface magnetic flux Φ. Each data point represents an average in time and space for each run. The values are normalized to model run 1B represented here by the subscript ‘ref’. The black line shows the powerlaw fit to the case where no quenching is included (i.e. q_{sat} = 1). The red and blue lines show powerlaw fits for q_{sat} = 2∕3 and 1/3. The corresponding powerlaw indices β for are given with the legend. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.