© ESO 2021

1 Introduction

The Sun has an aura of hot plasma called the corona, which has a temperature of a few million kelvin (Edlén 1943). The image of the corona has been captured by the Atmospheric Imaging Assembly (Lemen et al. 2012) of NASA’s Solar Dynamics Observatory (Pesnell et al. 2012). The corona is not unique to the Sun and has been observed to surround low-mass main-sequence stars in general (Güdel et al. 1997). Due to its high temperature, the corona is the principal source of stellar XUV (X-ray plus EUV: extreme ultraviolet) emissions.

Stellar XUV emissions drive the expansion and thermal escape of planetary atmospheres (Vidal-Madjar et al. 2003; Lecavelier des Etangs et al. 2012; Owen & Wu 2013; Ehrenreich et al. 2015; Airapetian et al. 2017). Therefore, describing the stellar XUV emission in terms of the stellar fundamental parameters (luminosity, mass, radius, etc.) is essential in understanding the evolution of a planet and its habitability. Since stellar XUV emissions decrease over time (Güdel et al. 1997; Ribas et al. 2005; Telleschi et al. 2005; Claire et al. 2012; Guinan et al. 2016), presumably in response to the stellar spin-down (Kraft 1967; Skumanich 1972; Barnes 2003; Irwin & Bouvier 2009; Matt et al. 2015) caused by magnetised stellar wind (Weber & Davis 1967; Kawaler 1988; Shoda et al. 2020), the long-term evolution of the XUV activity of the host star needs to be described as well (Tu et al. 2015; Johnstone et al. 2021).

While the observational characteristics of stellar X-ray emissions have been established, including their correlations with the rotation (Pallavicini et al. 1981; Wright et al. 2011) and magnetic field of the star (Pevtsov et al. 2003; Vidotto et al. 2014; Kochukhov et al. 2020), stellar EUV emissions are poorly characterised because they are difficult to observe. Extreme-ultraviolet photons suffer from strong absorption by the interstellar medium (Rumph et al. 1994), for which stellar EUV spectra are observable only for nearby stars in a limited range of wavelengths (≤360 Å, Ribas et al. 2005; Johnstone et al. 2021). One thus needs to indirectly estimate or reconstruct the EUV spectrum of a star from other observables. Proposed reconstruction methods include the inversion of the differential emission measure from UV and/or X-ray observations (Sanz-Forcada et al. 2011; Duvvuri et al. 2021), the empirical correlation between other observable lines and XUV emission (Linsky et al. 2014; Youngblood et al. 2017; France et al. 2018; Sreejith et al. 2020) and/or a combination of them (Diamond-Lowe et al. 2021). However, the empirical estimations of stellar EUV emission are based on observations with large uncertainty and a limited number of samples, and therefore require further validation from different perspectives.

To circumvent the intrinsic difficulty of stellar EUV observation, this study uses the solar-stellar connection to theoretically estimate the stellar XUV emission. The solar-stellar theoretical connection is a natural strategy for predicting stellar properties. For example, by extending the solar coronal theory, Shibata & Yokoyama (2002) modelled the X-ray characteristics of the stellar coronae, which was later extended by Takasao et al. (2020), who considered the size distribution of active regions. In this study, by extending the solar atmospheric model, the structure of the upper stellar atmosphere and the XUV emission are predicted. To this end, the coronal heating problem must be explicitly solved.

To solvethe coronal heating problem, the following three issues must be addressed (Klimchuk 2006). The first is energy generation and transfer. The source of the coronal thermal energy probably comes from the magneto-convection on the surface (Steiner et al. 1998). The energy generation and transfer to the corona, possibly in the form of Alfvén waves (Alfvén 1947; Osterbrock 1961; Kudoh & Shibata 1999; De Pontieu et al. 2007; McIntosh et al. 2011; Srivastava et al. 2017), needs to be solved. The second is energy dissipation in the corona. The magnetic (and kinetic) energy in the corona needs to dissipate to sustain the high-temperature corona. The field-braiding process must be considered as a promising mechanism of magnetic-energy dissipation (Parker 1972, 1988). The third is thermal response to the heating. The density and temperature of a coronal loop are determined by the energy balance among heating, conduction, and radiation. The Rosner-Tucker-Vaiana (RTV) scaling law (Rosner et al. 1978)originates from the thermal response to coronal heating. Thus, the RTV scaling law or its generalised form (Serio et al. 1981; Zhuleku et al. 2020) should be reproduced by the model (Antolin & Shibata 2010). For these issues, a model of the coronal heating should: (1) include the photosphere (stellar surface) and chromosphere, (2) appropriately consider the field-braiding process in the corona, and (3) implement the thermal conduction and radiative cooling.

Classically, numerical models of a corona have often focused on the thermal responses to heating events by one-dimensional (expanding) flux-tube models (Antiochos & Sturrock 1978; Peres et al. 1982; Antiochos et al. 1999). A zero-dimensional model of the coronal thermal evolution is also proposed (Klimchuk et al. 2008). The advancement of numerical techniques and increase in computational power have made models highly sophisticated. Several models deal with the realistic energy generation by explicitly solving the magneto-convection (Hansteen et al. 2015; Rempel 2017). Other models have focused more on energy dissipation with simpler numerical settings (Moriyasu et al. 2004; Rappazzo et al. 2008; van Ballegooijen et al. 2011; Dahlburg et al. 2016). These models have predicted that the signature of coronal heating could be explained by convection-driven energy injection. By extending these studies, we aimed to construct a stellar coronal model that satisfies the three requirements.

In deriving the X-ray and EUV spectra from simulations, care needs to be taken in the spatial resolution at the transition region (the temperature-jump region between the chromosphere and the corona). The numerical resolution around the transition region is found to significantly affect the coronal density (Bradshaw & Cargill 2013) and the coronal emission measure distribution (EMD). Because the coronal emission is proportional to the EMD, the XUV emission predicted by a simulation significantly depends on the numerical resolution. The required resolution is of the order of kilometres or less, which is impractical in realistic three-dimensional (3D) simulations. Previous studies have attempted to solve this ‘transition-region problem’ by introducing the artificial broadening of the transition region by tuning the magnitude of radiative cooling and thermal conduction (Johnston et al. 2017, 2021; Johnston & Bradshaw 2019; Iijima & Imada 2021). However, this treatment may yield an unrealistic EMD in the transition-region temperature and therefore is inappropriate for the calculation of the XUV emission that includes a significant contribution from the transition region.

Considering the difficulty of the transition-region problem, we perform a series of 1D, high-resolution magnetohydrodynamic (MHD) simulations for a wide range of parameters. This 1D model facilitates a coronal loop simulation with sufficiently high numerical resolution. The field-braiding process (or turbulence), which is essentially 3D, must be appropriately modelled when solving the coronal heating problem by 1D simulation. In this study, an approximated formulation of turbulent dissipation, developed in previous studies (Dmitruk et al. 2002; Shoda et al. 2018), is employed as it is likely to reproduce the average heating rate of the coronal loop (van Ballegooijen et al. 2011). For simplicity, we focus on the Sun-like stars that exhibit solar mass, luminosity, radius, and metallicity.

The rest of the manuscript is structured as follows. In Sect. 2, the coronal model and numerical methods are detailed. Section 3 produces the numerical results, focusing on the dependences of coronal properties on the loop length and magnetic filling factor that are likely to vary with the star (Reale & Micela 1998; Reale et al. 2004; Reiners et al. 2009; See et al. 2019). The XUV spectra obtained from the simulations are also presented. In Sect. 4, analytical arguments on the energy flux injected into the corona are presented. In Sect. 5, the interpretations and limitations of the proposed model are discussed. Section 6 summarises the study. Several details are provided in the appendix, including the resolution dependence of the model (see Appendix B).

2 Model

2.1 Model overview and notation

As mentioned earlier, the aim of this study is to model the XUV (X-ray+EUV) emissions from the Sun and Sun-like stars (including young Sun) with a range of magnetic activity level. To this end, the luminosity, mass, radius, and metallicity of the stars are fixed to the solar value. $$L=L_{\odot },M=M_{\odot },R=R_{\odot },Z=Z_{\odot }.$$(1)

Dependence on metallicity should appear in the radiative loss function (Sect. 2.5) and spectrum calculation (Sect. 3.1). We study the dependence of coronal properties on the coronal loop length and the filling factorof magnetic fields. Because the filling factor tends to increase with the stellar rotation rate (or the Rossby number, Saar 2001; Reiners et al. 2009), the coronal dependence on the stellar rotation rate is implicitly investigated.

We model a single coronal loop rooted in the stellar surface, and self-consistently solve the energetics and dynamics inside the loop. In other words, we solve the time-dependent, 1D MHD equations for an expanding coronal loop. The differential emission measure (DEM) of a single coronal loop is directly obtained from the simulation, which is then converted to the XUV spectrum by prescribing the chemical composition and integrating the continuum and line emissions asa function of wavelength using the CHIANTI atomic database version 10.0 (Dere et al. 1997; Del Zanna et al. 2021).

Notations of the coordinates and constant parameters used in this study are listed in Table 1. X_⊙ denotes the solar value and X_* the value of X measured atthe photosphere.

2.2 Model of the closed flux tube

A single coronal loop is modelled by a 1D expanding flux tube rooted in the photosphere. Figure 1 illustrates a schematic of the model. As shown in Table 1, the coordinate along the axis of the flux tube is denoted by s. The spatial variation only along the loop is assumed to be non-zero, that is, ∂∕∂s≠0. Similar 1D flux tube models were used in the models of solar spicules (Hollweg et al. 1982; Kudoh & Shibata 1999; Matsumoto & Shibata 2010), solar and stellar coronal heating (Moriyasu et al. 2004; Washinoue & Suzuki 2019), and solar wind acceleration (Suzuki & Inutsuka 2005; Shoda et al. 2018).

Hereinafter, the two perpendicular directions shall be denoted by x and y. Thus, the local xy plane is perpendicular to the axis of the loop. The flux-tube expansion is incorporated through the scale factors in the x and y directions: h_x,y. For simplicity, the loop is assumed to expand isotropically in the perpendicular directions. In terms of scale factors, the isotropic expansion isrepresented by $$h_x=h_y\propto \sqrt{A(s)},$$(2)

where A(s) is the cross section of the coronal loop. As h_s = 1 by definition, Eq. (2) results in $$\begin{array}{ll}\nabla \cdot X\hfill & =\frac{1}{A(s)}\frac{\partial }{\partial s}\left(X_{s}A(s)\right),\hfill \\ \nabla \times X\hfill & =\frac{1}{\sqrt{A(s)}}\frac{\partial }{\partial s}\left(X_x\sqrt{A(s)}\right)e_y\hfill \\ \hfill & -\frac{1}{\sqrt{A(s)}}\frac{\partial }{\partial s}\left(X_y\sqrt{A(s)}\right)e_x\hfill \end{array}$$(3)

for any vector field X, where e_x,y represent the unit vectors in the x, y directions. The 1D spherical coordinate system is reproduced as a special case of A(s) = s².

The flux tube expands in the chromosphere as a response to the exponential decrease in the ambient gas pressure (Cranmer & van Ballegooijen 2005; Ishikawa et al. 2021). As a result of this expansion, the filling factor of the magnetic field (flux tube) f should increase nearly exponentially with altitude. Under this assumption, we model the filling factor f as $$f=\min \left[1,f_{\ast }\exp \left(\frac{r-R_{\odot }}{H_{\text{mag}}}\right)\right],H_{\text{mag}}=c_{\text{mag}}{H}_{\ast },$$(4)

where f_* is the magnetic filling factor on the photosphere and $$H_{\ast }=\frac{k_{\text{B}}{T}_{\odot }{R}_{\odot }^2}{G{M}_{\odot }{m}_{\text{H}}}$$(5)

is the pressure scale height at the photosphere. By this formulation, we assume that the loop expands only in the chromosphere and exhibits a uniform cross section in the corona, which is supported by some solar observations (Klimchuk et al. 1992, but see a recent discussion by Malanushenko et al. 2021). Given that the pressure scale height is uniform from the photosphere up to the chromosphere, c_mag = 2 yields a flux-tube expansion with a constant plasma beta in altitude. In this work, we set c_mag = 2.5 that realises a slightly low-beta chromosphere. We have confirmed that the choice of c_mag does not have a significant influence over the simulation results.

When the coronal loop extends to a region far above the surface, the flux tube also undergoes the radial expansion ∝ r², where r is the radial distance from the stellar centre. Considering the chromospheric and radial expansions, the cross section A is expressed as $$A\propto r^{2}f.$$(6)

We define the inclination of the flux tube by prescribing r as a function of s. We considernearly vertical flux tubes, so that the vertical line of sight is nearly aligned with the axis of the flux tube (Fig. 1). It is easier to calculate the DEM along the vertical line of sight for this structure. In particular, we set $$\frac{\text{d}r}{\text{d}s}=\tanh \left[\frac{10\left(l_{\text{loop}}-s\right)}{l_{\text{loop}}}\right],{r|}_{s=0}=R_{\odot }$$(7)

where l_loop is the half-loop length. The actual shape of the flux-tube axis is displayed in Fig. 2. Combining Eqs. (4)–(7), the cross sectionA(s) is well defined as a function of s.

Fig. 1 Schematic picture of the system. One-dimensional dynamics along the axis of the closed flux tube is simulated. The axis isindicated by the dotted line, while the flux tube surface is denoted by green lines. The flux tube is intended to be nearly vertical and super-radially expanding. The geometry of the loop is defined by the filling factor f and radial distance r as functions of field-aligned distance s.

Fig. 2 Shape of the flux-tube axis defined by Eq. (7). x and z denote the horizontal and vertical coordinates, respectively.

2.3 Basic equations

The MHD equations with the equation of state of partially ionised hydrogen, gravity, thermal conduction, radiative cooling, and phenomenology of turbulent heating are selected as the basic equations of the model, which are expressed in the form of conservation law (for derivation, see Appendix A) $$\frac{\partial }{\partial t}U+\frac{1}{r^{2}f}\frac{\partial }{\partial s}\left(F{r}^{2}f\right)=S.$$(8)

The conserved variables U and the corresponding fluxes F are given by $$U=\left(\begin{array}{c}\rho \\ \rho v_s\\ \rho v_x\\ \rho v_y\\ B_x\\ B_y\\ e\end{array}\right),F=\left(\begin{array}{c}\rho v_s\\ \rho v_s^2+p_T\\ \rho v_s{v}_x-B_s{B}_x/(4\pi )\\ \rho v_s{v}_y-B_s{B}_y/(4\pi )\\ v_s{B}_x-v_x{B}_s\\ v_s{B}_y-v_y{B}_s\\ \left(e+p_T\right)v_s-B_s\left(v_{\perp }\cdot B_{\perp }\right)/(4\pi )\end{array}\right),$$(9)

where $$v_{\perp }=v_x{e}_x+v_y{e}_y,B_{\perp }=B_x{e}_x+B_y{e}_y,$$(10) $$p_T=p+\frac{B_{\perp }^2}{8\pi },e=e_{\text{int}}+\frac{1}{2}\rho v^2+\frac{B_{\perp }^2}{8\pi }.$$(11)

Here, e_int denotes the internal energy per unit volume and is defined in Sect. 2.4. The source term S is given by $$S=\left(\begin{array}{c}0\\ \left(p+\frac{1}{2}\rho v_{\perp }^2\right)/L-\rho \frac{G{M}_{\odot }}{r^2}\frac{\text{d}r}{\text{d}s}\\ \frac{1}{2L}\left(-\rho v_s{v}_x+\frac{B_s{B}_x}{4\pi }\right)+\rho D_x^v\\ \frac{1}{2L}\left(-\rho v_s{v}_y+\frac{B_s{B}_y}{4\pi }\right)+\rho D_y^v\\ \frac{1}{2L}\left(v_s{B}_x-v_x{B}_s\right)+\sqrt{4\pi \rho }{D}_x^b\\ [1.5pt]\frac{1}{2L}\left(v_s{B}_y-v_y{B}_s\right)+\sqrt{4\pi \rho }{D}_y^b\\ -\rho v_s\frac{G{M}_{\ast }}{r^2}\frac{\text{d}r}{\text{d}s}+Q_{\text{cnd}}+Q_{\text{rad}}\end{array}\right),$$(12)

where $L^{-1}=\text{d}/\text{d}s\ln \left(r^{2}f\right)$ denotes the length scale of the flux-tube expansion. The conduction term Q_cnd is defined in terms of conductive flux q_cnd as $$Q_{\text{cnd}}=-\frac{1}{r^{2}f}\frac{\partial }{\partial s}\left(q_{\text{cnd}}{r}^{2}f\right).$$(13)

Because the mean free path of an electron is generally smaller than the system size, the Spitzer–Härm flux (Spitzer & Härm 1953) is applied to q_cnd: $$q_{\text{cnd}}=-\frac{\left|B_s\right|}{\left|B\right|}{\kappa }_{\text{SH}}{T}^{5/2}\frac{\partial T}{\partial s},$$(14)

where κ_SH = 10⁻⁶ erg cm⁻¹ s⁻¹ K^−7∕2. The radiative cooling Q_rad and turbulent dissipation $D_{x,y}^{v,b}$ are described in Sects. 2.5 and 2.6, respectively. Because turbulence is considered not as an external force but as a dissipative process, the losses of the kinetic and magnetic energies by turbulence are locally balanced by the gain in the internal energy. Thus, the presence of the turbulence terms does not affect the conservation of the total energy. For the same reason, we do not explicitly consider the numerical dissipation of velocity and magnetic field in the energy equation. We note that the numerical dissipation is unlikely to be the dominant heating mechanism because the coronal Alfvén wave, which has a typical wavelength of ~ 100Mm, is resolved by a sufficiently fine grid in the corona (100km).

2.4 Equation of state

We assume that the plasma consists of neutral hydrogen atoms, protons, and electrons. The internal energy per unit volume e_int is composed of the conventionalthermal energy p∕(γ − 1) and the latent heat of the ionised gas: $$e_{\text{int}}=\frac{p}{\gamma -1}+n_{\text{H}}\chi I_{\text{H}},n_{\text{H}}=\rho /m_{\text{H}},$$(15)

where χ is the ionisation degree and n_H is the number density of hydrogen atoms (proton + neutral hydrogen). I_H is the ionisation energy of the hydrogen atom (I_H = 13.6eV). We assume that the ionisation degree could be determined from the approximated version of the Saha-Boltzmann equation, in which only the ground state is considered as the bound state (low-temperature limit). $$\frac{{\chi }^2}{1-\chi }=\frac{2}{n_{\text{H}}{\lambda }_{\text{e}}^3}\exp \left(-\frac{I_{\text{H}}}{k_{\text{B}}T}\right),$$(16)

where λ_e is the thermal de Broglie wavelength of electron. $${\lambda }_{\text{e}}=\sqrt{\frac{h^2}{2\pi m_{\text{e}}{k}_{\text{B}}T}}.$$(17)

When chromospheric hydrogen is no longer in thermal equilibrium, the ionisation degree will deviate from the Saha–Boltzmann value (Goodman & Judge 2012), which is beyond the scope of this study. Once the ionisation degree χ is obtained, the pressure and temperature are related by $$p=\left(n_{\text{e}}+n_{\text{H}}\right)k_{\text{B}}T=\left(1+\chi \right)n_{\text{H}}{k}_{\text{B}}T.$$(18)

2.5 Radiation

The radiative cooling rate per unit volume Q_rad is given by $$Q_{\text{rad}}={\xi }_{\text{rad}}{Q}_{\text{rad}}^{\text{thck}}+\left(1-{\xi }_{\text{rad}}\right)Q_{\text{rad}}^{\text{thin}},$$(19) $${\xi }_{\text{rad}}=1-\exp \left(-\frac{p}{p_{\text{rad}}}\right),p_{\text{rad}}/p_{\ast }=0.1,$$(20)

where p_* is the pressure at the surface (photosphere) and $Q_{\text{rad}}^{\text{thck}}$ and $Q_{\text{rad}}^{\text{thin}}$ approximatethe optically thick and thin cooling rates, respectively. The optically thick and thin functions are seamlessly connected via ξ_rad. Instead of solving the radiative transfer, we model $Q_{\text{rad}}^{\text{thck}}$ and $Q_{\text{rad}}^{\text{thin}}$ as follows.

The radiative heating and cooling are approximately balanced near the photosphere to maintain a nearly constant surface temperature. The optically thick radiative loss is approximated by an exponential cooling function that forces the local temperatureto approach the reference value (Gudiksen & Nordlund 2005): $$Q_{\text{rad}}^{\text{thck}}=\frac{1}{\tau }\left(e_{\text{int}}^{\text{ref}}-e_{\text{int}}\right),\tau =0.1\text{s}{\left(\frac{\rho }{{\rho }_{\ast }}\right)}^{-1/2},$$(21)

where ρ_* is the photospheric mass density and $e_{\text{int}}^{\text{ref}}$ is the reference internal energy density corresponding to the reference temperature T^ref. We simply assume T^ref = T_⊙, namely, the stellar atmosphere exhibits isothermal behaviour in the absence of other cooling/heating mechanisms.

The optically thin cooling is expressed in terms of the radiative loss function Λ(T) by n_e n_HΛ(T). In this work, we define the loss function over a wide range of temperature (10³ K ≤ T ≤ 10⁷ K) as follows. In the high-temperature range (T ≥ 1.5 × 10⁴ K), the cooling rate is referred from the CHIANTI atomic database with the photospheric abundance (no first ionisation potential (FIP) effect). In the low-temperature range (T ≤ 1.0 × 10⁴ K), the cooling function is deduced by Goodman & Judge (2012), which partially consider the non-LTE effect. In the intermediate-temperature range (1.0 × 10⁴ K < T < 1.5 × 10⁴ K), a bridging law between two loss functions asre employed following the method of Iijima (2016). The chromospheric heating effect by backward coronal radiation is still missing in Λ(T) defined above. To introduce this effect, we quench Λ(T) in the chromospheric temperature range, which gives $Q_{\text{rad}}^{\text{thin}}$ $$Q_{\text{rad}}^{\text{thin}}=n_{\text{e}}{n}_{\text{H}}{\Lambda }_{\text{eff}}(T),{\Lambda }_{\text{eff}}(T)=\Lambda (T)\exp \left(-\frac{T_{\text{chr}}^2}{T^2}\right),$$(22)

where T_chr = 2.0 × 10⁴ K. The effective radiative loss function Λ_eff(T) is displayed in Fig. 3, along with the original radiative loss function Λ(T) and the CHIANTI loss function defined in T ≥ 10⁴ K.

2.6 Phenomenological model of coronal turbulence

Although the mechanism of the coronal heating is debated, it is of no doubt that magnetic field feeds heat to the corona that maintains the million-Kelvin temperature by dissipation. The formation of tangential discontinuities (electric current sheets) in response to the continuous shuffling of the foot points of coronal magnetic fields is a plausible mechanism of magnetic-field dissipation (Parker 1972, 1983; Sturrock & Uchida 1981; van Ballegooijen 1986; Galsgaard & Nordlund 1996). The ubiquitous current-sheet formation should lead to small-scale impulsive energy release, which is likely to feed a sufficient amount of energy to the corona, possibly in the form of micro- and nano-flares (Parker 1988; Shimizu 1995; Aschwanden & Parnell 2002).

The formation of current sheets can be interpreted as turbulent cascading (Rappazzo et al. 2007, 2008; Verdini et al. 2012). Because the coronal loop is threaded by a strong mean magnetic field, the MHD turbulence evolving in the coronal loop can be accurately approximated by the reduced-MHD turbulence, in which the energy cascades preferentially in the perpendicular direction (Shebalin et al. 1983; Cho & Vishniac 2000; Cho & Lazarian 2003). The energy-cascading (or heating) rate of the reduced-MHD turbulence, Q_heat, is precisely approximated by the mean-field quantities as follows (Hossain et al. 1995; Matthaeus et al. 1999; Dmitruk et al. 2002; Verdini & Velli 2007) $$Q_{\text{heat}}\approx c_{\text{d}}\rho \frac{z_{\perp }^{+}{z}_{\perp }^{-}{}^2+z_{\perp }^{-}{z}_{\perp }^{+}{}^2}{4{\lambda }_{\perp }},$$(23)

where $z_{\perp }^{\pm }$ denotes the (root-mean-squared (RMS)) amplitude of the perpendicular Elsässer variables and λ_⊥ is the correlation length of the Elsässer variable (Alfvén wave) perpendicular to the mean field. c_d is a dimensionless parameter. By this approximation, we estimate the averaged heating rate using the coronal turbulence (field braiding).

The approximated heating rate in Eq. (23) is implemented by adding the source terms $D_{x,y}^v$, $D_{x,y}^b$ given by Shoda et al. (2018). $$D_{x,y}^v=-\frac{c_{\text{d}}}{4{\lambda }_{\perp }}\left(\left|z_{x,y}^{+}\right|z_{x,y}^{-}+\left|z_{x,y}^{-}\right|z_{x,y}^{+}\right),$$(24) $$D_{x,y}^b=-\frac{c_{\text{d}}}{4{\lambda }_{\perp }}\left(\left|z_{x,y}^{+}\right|z_{x,y}^{-}-\left|z_{x,y}^{-}\right|z_{x,y}^{+}\right),$$(25)

where $z_{x,y}^{\pm }=v_{x,y}\mp B_{x,y}/\sqrt{4\pi \rho }$. The role of these terms is explained below. Without the conservation part $\propto \partial \left(r^{2}fF\right)/\partial s$, the perpendicular components of the equation of motion and induction equation are expressed as $$\frac{\partial }{\partial t}\left(\rho v_{x,y}\right)=\rho D_{x,y}^v,\frac{\partial }{\partial t}{B}_{x,y}=\sqrt{4\pi \rho }{D}_{x,y}^b.$$(26)

In the limit of the reduced-MHD approximation (time-independent density, ∂ρ∕∂t = 0), Eqs. (24)–(26) are reduced to $$\frac{\partial }{\partial t}{z}_{x,y}^{\pm }=-\frac{c_{\text{d}}}{2{\lambda }_{\perp }}\left|z_{x,y}^{\mp }\right|z_{x,y}^{\pm },$$(27)

which yields the energy conservation law of $$\frac{\partial }{\partial t}{e}_{\perp }=-c_{\text{d}}\rho {\sum }_{i=x,y}\frac{\left|z_i^{-}\right|z_i^{+}{}^2+\left|z_i^{+}\right|z_i^{-}{}^2}{4{\lambda }_{\perp }},$$(28)

where e_⊥ is the sum of the kinetic and magnetic energies emerging from the fluctuations of the perpendicular components: $$e_{\perp }=\frac{1}{4}\rho \left(z_{\perp }^{+}{}^2+z_{\perp }^{-}{}^2\right)=\frac{1}{2}\rho v_{\perp }^2+\frac{B_{\perp }^2}{8\pi }.$$(29)

Comparing Eqs. (23) and (28), one obtain $$\frac{\partial }{\partial t}{e}_{\perp }\approx -Q_{\text{heat}},$$(30)

indicating that the energy dissipation by the reduced-MHD turbulence is considered appropriately.

The perpendicular correlation length is assumed to increase with the flux-tube radius, that is, $${\lambda }_{\perp }={\lambda }_{\perp ,\ast }\sqrt{\frac{A}{A_{\ast }}}={\lambda }_{\perp ,\ast }\frac{r}{R_{\odot }}\sqrt{\frac{f}{f_{\ast }}},$$(31)

where the perpendicular correlation length at the photosphere is set equal to the typical width of the inter-granular lane: λ_⊥,* = 150 km. However, the best possible free parameter c_d is still debated. In this study, we infer c_d = 0.1 from the previous studies of the solar-wind turbulence (van Ballegooijen & Asgari-Targhi 2017; Chandran & Perez 2019; Verdini et al. 2019).

Fig. 3 Effective and original optically thin radiative loss functions, Λeff(T) and Λ(T), shown with solid and dashed lines, respectively. Also shown, by crosses, are the loss function from the CHIANTI atomic database with photospheric abundance.

2.7 Boundary condition and simulation setting

Both boundaries of the simulation domain are located at the photosphere, and the photospheric temperature is fixed to the effective temperature: $$T_{\ast }=T_{\odot }=5.77\times {10}^3\text{K}.$$(32)

The photospheric magnetic field is known to form localised kilo-Gauss patches (Spruit & Zweibel 1979; Tsuneta et al. 2008). These patches are likely to be in thermal equipartition, which equates the gas and magnetic pressures. For the non-magnetised photosphere, the thermal equipartition field is given by $$B_{\text{eq}}=1.34\times {10}^3\text{G}.$$(33)

In the magnetised photosphere, as the deeper region tends to be observed (e.g. Keller et al. 2004), the ambient gas is likely to exhibit larger pressure than the non-magnetised photosphere. Thus, the equipartition magnetic field should be larger than this equipartition value. Therefore, we set the photospheric axial magnetic field equal to $$B_{\text{s},\ast }=1.5B_{\text{eq}}=2.01\times {10}^3\text{G}.$$(34)

We modelled the energy injection from the photosphere by imposing the velocity and magnetic-field fluctuations. Fluctuations were imposed at both ends of the simulation domain. The vertical and horizontal velocity fluctuations were modelled separately. The upward acoustic waves were excited on the photosphere by employing the time-dependent boundary conditions on the density and axial velocity. $${\rho }_{\ast }={\overline{\rho }}_{\ast }\left(1+\frac{v_{\text{s},\ast }}{a_{\ast }}\right),$$(35) $$v_{\text{s},\ast }={\displaystyle {\int }_{{\omega }_{\text{min}}^l}^{{\omega }_{\text{max}}^l}\text{d}}\omega {\tilde{v}}_{\text{s}}\left(\omega \right)\sin \left[\omega t+{\psi }^l(\omega )\right],$$(36)

where ${\overline{\rho }}_{\ast }$ is the time-averaged photospheric density, $a_{\ast }=\sqrt{k_{\text{B}}{T}_{\ast }/m_{\text{H}}}$ is the isothermal speed of sound on the photosphere, and ${\psi }^l\left(\omega \right)$ is a random phase function that ranges between 0 and 2π. In the numerical implementation of the integral in Eq. (35), the frequency range was evenly divided into 21 bins and the corresponding 21 components were summed. The time-averaged photospheric density is given by equipartition on the photosphere: $${\overline{\rho }}_{\ast }{k}_{\text{B}}{T}_{\ast }/m_{\text{H}}=\frac{B_{\text{s},\ast }^2}{8\pi },$$(37)

which yields $${\overline{\rho }}_{\ast }=4.22\times {10}^{-7}\text{g}{\text{cm}}^{\text{-3}}.$$(38)

Here, ${\overline{\rho }}_{\ast }$ is larger than the typical mass density on the solar surface because the magnetised photosphere should be deeper and denser. The (time-averaged) Alfvén speed v_A on the photosphere is then expressed as $$\overline{v_{\text{A},\ast }}\approx \frac{B_{\text{s},\ast }}{\sqrt{4\pi \overline{{\rho }_{\ast }}}}=8.73\text{km}{\text{s}}^{-1}.$$(39)

We arbitrarily set ${\tilde{v}}_{\text{s}}\left(\omega \right)\propto {\omega }^{-1/2}$ with $$2\pi /{\omega }_{\text{min}}^l=300\text{s},2\pi /{\omega }_{\text{max}}^l=100\text{s}.$$(40)

The minimum frequency corresponds to the cut-off frequency of the acoustic wave at the photosphere (e.g. Felipe et al. 2018). The magnitude of ${\tilde{v}}_{\text{s}}\left(\omega \right)$ was tuned such that the RMS amplitude of v_s,* at the photosphere is 0.6 km s⁻¹. $$\sqrt{\overline{v_{\text{s},\ast }^2}}=0.6\text{km}{\text{s}}^{-1},$$(41)

where the overline denotes the time average. Although the longitudinal-wave excitation on the photosphere is explicitly considered, the effect of the longitudinal-wave input is insignificant; the coronal temperature decreases by only 2% when the longitudinal wave injection is terminated.

The horizontal velocity and magnetic field at the bottom boundary are expressed in terms of the Elsässer variables, which are defined as $$z_{x,y}^{\pm }=v_{x,y}\mp \frac{B_{x,y}}{\sqrt{4\pi \rho }}.$$(42)

The free boundary condition was imposed on the downward Elsässer variables. $${\frac{\partial }{\partial s}{z}_{x,y}^{-}|}_{\ast }=0.$$(43)

The upward Elsässer variable was assumed to be non-monochromatic with respect to the frequency. $$z_{x,y,\ast }^{+}={\displaystyle {\int }_{{\omega }_{\text{min}}^t}^{{\omega }_{\text{max}}^t}\text{d}}\omega {\tilde{z}}_{x,y}^{\pm }\left(\omega \right)\sin \left[\omega t+{\psi }^t(\omega )\right],$$(44)

where ψ^t(ω) is a random phase function that ranges between 0 and 2π. As with Eq. (35), the frequency range was discretised into 21 bins when numerically implementing the integral in Eq. (44). We set ${\tilde{z}}_{x,y}^{\pm }\left(\omega \right)\propto {\omega }^{-1/2}$ with $$2\pi /{\omega }_{\text{min}}^t=1000\text{s},2\pi /{\omega }_{\text{max}}^t=100\text{s}.$$(45)

${\tilde{z}}_{x,y}^{\pm }\left(\omega \right)\propto {\omega }^{-1/2}$ corresponds to the 1∕ω energy spectrum discovered in previous solar simulations and observations (Van Kooten & Cranmer 2017). Given that the typical size of a granule is 1000 km and the typical speed of surface convection is 1 km s⁻¹ (Chitta et al. 2012), the maximum wave period corresponds to the turn-over time of granular motion. Similarly, given that the typical size of an inter-granular lane is 100 km, the minimum wave period corresponds to the turn-over time of inter-granular motion. The magnitude of ${\tilde{z}}_{x,y}^{\pm }\left(\omega \right)$ was tuned such that the root-mean-squared amplitude of $z_{x,y,\ast }^{+}$ at the photosphere is 1.2 km s⁻¹: $$\sqrt{\overline{z_{x,\ast }^{+}{}^2}}=\sqrt{\overline{z_{y,\ast }^{+}{}^2}}=1.2\text{km}{\text{s}}^{-1},$$(46)

where the overline denotes the time average. Although some solar observations have observed the suppression of convective velocity in large-filling-factor regions (e.g. Katsukawa & Tsuneta 2005), we dismiss this effect for simplicity.

By this formulation, the energy flux of the upward Alfvén wave at the footpoint of the flux tube is given by $$\begin{array}{ll}{F}_{\text{A},\ast }\hfill & =\frac{1}{4}\overline{\rho \left(z_{x,\ast }^2+z_{y,\ast }^2\right)v_{\text{A},\ast }}\approx \frac{1}{4}\overline{{\rho }_{\ast }}\left(\overline{z_{x,\ast }^2}+\overline{z_{y,\ast }^2}\right)\overline{v_{\text{A},\ast }}\hfill \\ \hfill & =2.65\times {10}^9\text{erg}{\text{cm}}^{-2}{s}^{-1},\hfill \end{array}$$(47)

which is sufficiently larger than the energy flux required to sustain the solar corona (Withbroe & Noyes 1977).

2.8 Numerical method

A non-uniform grid system was used to resolve the computational domain 0 ≤ s ≤ 2l_loop. A uniform cell size of Δs_min was used below the critical height s < s_ge, above which the cell size expands until it reaches the maximum Δs_max. In particular, in 0 ≤ s ≤ l_loop, the size of the i-th cell, Δs_i, was iteratively defined as $$\begin{array}{l}\Delta s_i=\max \left[\Delta s_{\text{min}},\min \left[\Delta s_{\text{max}},\Delta s_{\text{min}}+\frac{2{\epsilon }_{\text{ge}}}{2+{\epsilon }_{\text{ge}}}\left(s_{i-1}-s_{\text{ge}}\right)\right]\right],\\ s_i=s_{i-1}+\frac{1}{2}\left(\Delta s_{i-1}+\Delta s_i\right).\end{array}$$(48)

Letting N be the total number of cells, we express the cell size in the latter half of the domain l_loop ≤ s ≤ 2l_loop by $$\Delta s_i=\Delta s_{N+1-i},s_i=s_{i-1}+\frac{1}{2}\left(\Delta s_{i-1}+\Delta s_i\right).$$(49)

The maximum cell size was fixed to Δs_max = 100 km. The relation between the minimum cell size and f_* is $$\Delta s_{\text{min}}=5\text{km}(f_{\ast }<0.05),\Delta s_{\text{min}}=2\text{km}(f_{\ast }\ge 0.05).$$(50)

This relation is used because a higher resolution is required at the transition region in the large-f_* runs (for details, see Appendix B). The grid expansion rate also depends on f_* as ε_ge = 2.13  (f_* < 0.05) and ε_ge = 1.89  (f_*≥ 0.05). The grid expansion height is fixed to s_ge = 10 000 km, which is greater than the typical height of the transition region that needs to be resolved with the minimum cell size.

In numerically solving Eq. (8), we rewrite the basic equations in terms of the cross-section-weighted conserved variables $\tilde{U}$ and the corresponding fluxes $\tilde{F}$ defined by $$\tilde{U}=\left(\begin{array}{c}\tilde{\rho }\\ \tilde{\rho }{\tilde{v}}_s\\ \tilde{\rho }{\tilde{v}}_x\\ \tilde{\rho }{\tilde{v}}_y\\ {\tilde{B}}_x\\ {\tilde{B}}_y\\ \tilde{e}\end{array}\right)=\left(\begin{array}{c}\rho r^{2}f\\ \rho v_s{r}^{2}f\\ \rho v_x{r}^{2}f\\ \rho v_y{r}^{2}f\\ B_{x}r\sqrt{f}\\ B_{y}r\sqrt{f}\\ e{r}^{2}f\end{array}\right),$$(51) $$\tilde{F}=\left(\begin{array}{c}\tilde{\rho }{\tilde{v}}_s\\ \tilde{\rho }{\tilde{v}}_s^2+{\tilde{p}}_T\\ \tilde{\rho }{\tilde{v}}_s{\tilde{v}}_x-{\tilde{B}}_r{\tilde{B}}_x/(4\pi )\\ \tilde{\rho }{\tilde{v}}_s{\tilde{v}}_y-{\tilde{B}}_r{\tilde{B}}_y/(4\pi )\\ {\tilde{v}}_s{\tilde{B}}_x-{\tilde{v}}_x{\tilde{B}}_r\\ {\tilde{v}}_s{\tilde{B}}_y-{\tilde{v}}_y{\tilde{B}}_r\\ \left(\tilde{e}+{\tilde{p}}_T\right){\tilde{v}}_s-{\tilde{B}}_r\left({\tilde{v}}_{\perp }\cdot {\tilde{B}}_{\perp }\right)/(4\pi )\end{array}\right),$$(52)

where $${\tilde{p}}_T=\tilde{p}+\frac{{\tilde{B}}_{\perp }^2}{8\pi }=p_T{r}^{2}f.$$(53)

Using $\tilde{U}$ and $\tilde{F}$, the basic equation is given by $$\frac{\partial }{\partial t}\tilde{U}+\frac{\partial }{\partial s}\tilde{F}=\tilde{S},$$(54)

where $$\tilde{S}=\left(\begin{array}{c}0\\ \left(\tilde{p}+\frac{1}{2}\tilde{\rho }{\tilde{v}}_{\perp }^2\right)/L-\tilde{\rho }\frac{G{M}_{\ast }}{r^2}\frac{\text{d}r}{\text{d}s}\\ \frac{1}{2L}\left(-\tilde{\rho }{\tilde{v}}_s{\tilde{v}}_x+\frac{{\tilde{B}}_r{\tilde{B}}_x}{4\pi }\right)+\tilde{\rho }{D}_x^v\\ \frac{1}{2L}\left(-\tilde{\rho }{\tilde{v}}_s{\tilde{v}}_y+\frac{{\tilde{B}}_r{\tilde{B}}_y}{4\pi }\right)+\tilde{\rho }{D}_y^v\\ \sqrt{4\pi \tilde{\rho }}{D}_x^b\\ \sqrt{4\pi \tilde{\rho }}{D}_y^b\\ -\tilde{\rho }{\tilde{v}}_s\frac{G{M}_{\ast }}{r^2}\frac{\text{d}r}{\text{d}s}+Q_{\text{cnd}}{r}^{2}f+Q_{\text{rad}}{r}^{2}f\end{array}\right).$$(55)

With this variable conversion, any MHD solver designed for the Cartesian coordinate system can be directly applied to Eq. (54). In this study, the Harten–Lax–van Leer discontinuities (HLLD) approximated Riemann solver (Miyoshi & Kusano 2005) is used to calculate $\tilde{F}$ at the cell boundary. For spatial reconstruction, the fifth-order accurate monotonicity-preserving method (Suresh & Huynh 1997) is used to reconstruct the cross-section-weighted conserved variables $\tilde{U}$ in s ≤ s_ge and 2l_loop − s ≤ s_ge, whereas the monotonic upstream-centred scheme for the law of conservation (van Leer 1979) with a minmod flux limiter is used in s_ge < s < 2l_loop − s_ge.

Thermal conduction and the other parts are solved independently by the second-order operator-splitting procedure as follows. First, the thermal conduction is solved for a half step Δt∕2. $${\tilde{U}}^n\stackrel{\Delta t/2}{\to }{\tilde{U}}^{\ast }(\text{thermal}\text{conduction}\text{only}),$$(56)

where ${\tilde{U}}^n$ is the n-th step value of $\tilde{U}$. Then, the rest of the basic equations are solved for a full step Δt. $${\tilde{U}}_i^{\ast }\stackrel{\Delta t}{\to }{\tilde{U}}_i^{\ast \ast }(\text{without}\text{thermal}\text{conduction}).$$(57)

Finally, the thermal conduction is solved again for a half step Δt∕2. $${\tilde{U}}_i^{\ast \ast }\stackrel{\Delta t/2}{\to }{\tilde{U}}_i^{n+1}(\text{thermal}\text{conduction}\text{only}).$$(58)

With this procedure, we avoid the severe constraints on the Δt from thermal conduction when updating the MHD equations.

The third-order strong-stability-preserving (SSP) Runge–Kutta method is used in the time integration of the MHD equations (Shu & Osher 1988; Gottlieb et al. 2001). The super-time-stepping method (Meyer et al. 2012, 2014) is used to solve the thermal conduction, which reduce the numerical cost and time with minimum loss of accuracy.

3 Simulation result

3.1 Fiducial (solar) case: atmosphere and spectrum

First, we discuss the simulation run with l_loop = 20 Mm and f_* = 1.0 × 10⁻² as the fiducial case. In this case, the coronal field strength is ≈20 G, which is within the range of the solar coronal magnetic field strength measured by the coronal seismology technique (Nakariakov & Ofman 2001; Verwichte et al. 2004; Jess et al. 2016), and thus the fiducial case is regarded as the solar case.

Figure 4 illustrates the time-averaged properties of a quasi-steady coronal loop. Panels show the mass density (top) and temperature (middle) along the loop axis and the differential emission measure (DEM, bottom), defined by $$\text{DEM}(T)=n_{\text{e}}^2\left(T\right)\frac{\text{d}{l}_{\text{los}}}{\text{d}T},$$(59)

where n_e(T) is the electron density withtemperature T and l_los is the lengthalong the line of sight. Practically, dividing the temperature range into bins and considering the vertical line of sight (l_los = r), the DEM and associated emission measure distribution (EMD) are numerically obtained as follows $$\text{DEM}(T_i)\Delta T_i\equiv \text{EMD}(T_i)=n_{\text{e}}^2\left(T_i\right)\Delta r(T_i),$$(60)

where n_e(T_i) is the total number density of an electron that exhibits a temperature in [T_i − ΔT_i∕2, T_i + ΔT_i∕2] and Δr(T_i) is the total radial extension of where T_i − ΔT_i∕2 ≤ T < T_i + ΔT_i∕2. The DEM is calculated in the temperature range of 10⁴ K ≤ T ≤ 10⁷ K with an equal spacing in the logarithmic scale of T. The DEM in the range of T < 10⁴ K is not calculated because, in the low-temperature range, the atmosphere is not optically thin and the DEM loses its meaning. We note that the unit of EMD is cm⁻⁵, while the ‘volume’ EMD, which has often been used in the literature (Güdel et al. 2003; Scelsi et al. 2005), represents the distribution of an emission measure over the whole coronal volume and has a different unit (cm⁻³).

The high-temperature (T > 10⁶ K) corona is successfully reproduced in the fiducial case. Since we are imposing the phenomenological, mean-field formulation of the coronal heating (field braiding/turbulence), the heating tends to be more constant in time and more uniform in space than the actual three-dimensional case in which the heating is intermittent in time and space. Nevertheless, the time-averaged heating rateshould be similar between the 1D approximation and the 3D simulation, because the previous 3D simulation of solar wind yielded a similar mean field to the 1D simulation with turbulence phenomenology (Shoda et al. 2019).

Under the assumption that the upper atmosphere (T ≥ 10⁴ K) is optically thin in the wavelength of interest, the XUV spectrum is obtained from the EMD using the open-source package ChiantiPy based on CHIANTI database ver 10.0 (Del Zanna et al. 2021). In particular, the specific intensity I_λ was calculated from the EMD using the ChiantiPy.core.Spectrum module with the coronal abundance given by Schmelz et al. (2012). Although the radiative loss functionΛ(T) is constructed with photospheric abundance, because the loss function is nearly independent of the FIP effect in the radiation-dominated temperaturerange T ≤ 3 × 10⁵ K, the inconsistency in abundance do not violate the simulation results. Assuming that the stellar corona is a uniformly bright sphere, the spectral flux density at the heliocentric distance r is deduced by (see, for example, Sect. 1.3 of Rybicki & Lightman 1979) $$F_{\lambda }=\pi I_{\lambda }{\left(\frac{R_{\odot }}{r}\right)}^2,$$(61)

which yieldsthe X-ray and EUV luminosities as $$L_{\text{X}}=4\pi r^2{\displaystyle {\int }_{5\AA }^{100\AA }\text{d}}\lambda F_{\lambda },L_{\text{EUV}}=4\pi r^2{\displaystyle {\int }_{100\AA }^{912\AA }\text{d}}\lambda F_{\lambda }.$$(62)

In terms of energy, X-ray photons are in the range of 0.12−2.48 keV. Caution must be exercised when calibrating the X-rays because a subtle difference in the bandpass of the instrument can result in large differences in the derived response function and X-ray luminosity (Zhuleku et al. 2020). The procedure for the translation between different instruments is available in Judge et al. (2003).

To test the capability of our model in the prediction of XUV spectrum, we compare in Fig. 5 the spectral flux density obtained from our simulation (red) and that from the observation in a solar activity minimum (blue). The observed spectrum is obtained from the coordinated observation in the Whole Heliosphere Interval (WHI, from March 20, 2008, to April 16, 2008, Woods et al. 2009; Chamberlin et al. 2009). Figure 5 shows that, below the Lyman edge (≤ 912 Å), the simulated spectrum is in a good agreement with the observed spectrum. Above the Lyman edge (≥ 912 Å), the continuum is underestimated, and the emission lines are overestimated in the simulated spectrum, possibly because the optically thin approximation is inadequate in this wavelength range. In this study, because the focus is on the spectrum below the Lyman edge, the simulation is validated with respect to spectrum prediction.

Fig. 4 Simulated loop properties of the fiducial case. Top: time-averaged (solid line) and snapshot (dashed line) profiles of mass density along the loop axis. Middle: time-averaged (solid line) and snapshot (dashed line) profiles of temperature along the loop axis. Bottom: time-averaged differential emission measure (solid line). Annotations ‘TR’ in the middle panel indicate the locations of the transition region in the snapshot profile.

Fig. 5 Observed (blue) and simulated (red) spectral flux density of the Sun measured at 1 au. The observed spectrum is retrieved in the solar activity minimum. The two spectra are in a good agreement below the Lyman edge (≤912 Å).

3.2 Loop-length dependence: density and temperature

The coronal loop length is a fundamental parameter that affects the coronal density and temperature. Here, we show the relation between the time-averaged coronal properties and the coronal loop length. For simplicity, we fix the magnetic filling factor to the fiducial (solar) value: f_* = f_⊙ = 0.01. Hereinafter, the time-averaged value of X will be denoted by X_ave.

In the discussion on the behaviour of the coronal properties, the results must be compared with the analytical RTV scaling law (Rosner et al. 1978). We note that the half-loop length in the RTV scaling law denotes the length from the transition region to the apex of the loop, whereas l_loop stands for the length from the stellar surface to the apex of the loop. For better comparison, instead of l_loop, we use the effective half-loop length $l_{\text{loop}}^{\text{eff}}$, which denotes the coronal length, and is defined as $$l_{\text{loop}}^{\text{eff}}={\displaystyle {\int }_{s_{\text{TR}}}^{l_{\text{loop}}}\text{d}}s,$$(63)

where s_TR (< l_loop) is where T_ave = 10⁵ K.

Another factor that should be considered is the gravitational stratification. The effective half-loop length often exceeds the coronal pressure scale height. In such cases, the loop-top pressure in the original RTV scaling law is replaced by the coronal-base pressure (Serio et al. 1981). Given that the pressure is continuous across the transition region, we define the coronal-base pressure p_base as the pressure measured at T_ave = 10⁵ K. For the coronal electron density, both the loop-top value n_e,top and coronal-base value n_e,base are measured. In contrast to pressure, the density is discontinuous across the transition region, and therefore the definition of the coronal-base density is not trivial. Here, we define n_e,base as the value measured in the coronal base: $$n_{\text{e},\text{base}}\equiv \frac{1}{s_{\text{cor},2}-s_{\text{cor},1}}{\displaystyle {\int }_{s_{\text{cor},1}}^{s_{\text{cor},2}}{n}_{e,\text{ave}}}\text{d}s,$$(64)

where we set s_cor,1 = 10 Mm and s_cor,2 = 15 Mm. The loop-top density and temperature are given by $$n_{e,\text{top}}=n_e\left(s=l_{\text{loop}}\right),T_{\text{top}}=T\left(s=l_{\text{loop}}\right).$$(65)

Figure 6 shows the $l_{\text{loop}}^{\text{eff}}$ dependencesof the time-averaged coronal properties (density, temperature and pressure). The loop-top temperature and coronal-basepressure obey a power–law relation with respect to the effective half-loop length, which is formulated as $$T_{\text{top}}\propto l_{\text{loop}}^{\text{eff}}{}^{0.39},$$(66) $$p_{\text{base}}\propto l_{\text{loop}}^{\text{eff}}{}^{0.25}.$$(67)

The (generalised) RTV scaling law predicts that the loop-top temperature obeys the following relation $$T_{\text{top}}^{\text{RTV}}\propto {\left(p_{\text{base}}{l}_{\text{loop}}^{\text{eff}}\right)}^{1/3},$$(68)

where we dismiss the exponential correction term as it is negligible (Serio et al. 1981). A comparison of Eqs. (66)–(68) reveals that the simulation results are consistent with the RTV scaling law.

An alternative form of the RTV scaling law predicts a relation among the coronal energy flux F_cor, loop length $l_{\text{loop}}^{\text{eff}}$, and loop-top temperature$T_{\text{top}}^{\text{RTV}}$, and is expressed as $$T_{\text{top}}^{\text{RTV}}\propto {\left(F_{\text{cor}}{l}_{\text{loop}}^{\text{eff}}\right)}^{2/7}.$$(69)

A comparison of Eqs. (66) and (69) indicates that the energy flux injected into the corona is larger for longer loops. In terms of the heating rate per unit volume Q, the RTV predictions, $$T_{\text{top}}^{\text{RTV}}\propto Q^{2/7}{l}_{\text{loop}}^{\text{eff}}{}^{4/7},$$(70) $$p_{\text{base}}^{\text{RTV}}\propto Q^{6/7}{l}_{\text{loop}}^{\text{eff}}{}^{5/7},$$(71)

and simulation results from Eqs. (66) and (67) indicate that $Q=F_{\text{cor}}/l_{\text{loop}}^{\text{eff}}$ is a decreasing function of $l_{\text{loop}}^{\text{eff}}$. These conclusions shall be directly validated in the following section.

The bottom panel of Fig. 6 shows the variations in loop-top and coronal-base electron densities over $l_{\text{loop}}^{\text{eff}}$. Given that the RTV scaling law predicts a larger loop-top density at a higher loop-top temperature for a uniform corona, the decrease in the loop-topdensity is attributed to gravitational stratification. We note that the coronal-base density exhibits a weaker dependence on $l_{\text{loop}}^{\text{eff}}$ than the coronal-base pressure, which is contradictory if the coronal-base temperature is constant in $l_{\text{loop}}^{\text{eff}}$. It may be interpreted that the coronal-base temperature increases with $l_{\text{loop}}^{\text{eff}}$ in response to the increasing T_top with $l_{\text{loop}}^{\text{eff}}$.

Fig. 6 Coronal-property dependences on the effective loop length. Top: relation between the effective half loop length $l_{\text{loop}}^{\text{eff}}$ (see Eq. (63) for definition) and the time-averaged loop-top temperature (Ttop, red circles) and coronal-base pressure (pbase, blue diamonds). Bottom: relation between the effective half loop length $l_{\text{loop}}^{\text{eff}}$ and the time-averaged loop-top electron density (ne,top, red circles) and coronal-base electron density (ne,base, blue diamonds). In both panels, lines represent the power-law fittings to the symbols.

3.3 Loop-length dependence: energy flux and XUV emission

For a better interpretation of the behaviour of coronal properties, the variation in the energy flux entering the corona F_cor with the effective half-loop length $l_{\text{loop}}^{\text{eff}}$ must be revealed.

Directly measuring F_cor from the simulation data is, however, not a trivial task. To obtain F_cor, one needs to measure the energy flux at the base of the corona, which moves in time and is broadened after time averaging. Because a significant amount of energy is reflected back at the transition region, a slight difference in the position of the coronal base yields a significant error or uncertainty in F_cor. Therefore, instead of directly measuring F_cor, we measure the backward conductive flux F_cnd, which should be balanced with F_cor by energy conservation.

The top and bottom panels in Fig. 7 show the loop-length (l_loop and $l_{\text{loop}}^{\text{eff}}$) dependence of the coronal conductive flux. The panels display the simulation runs with a fixed magnetic filling factor of f_* = 1.00 × 10⁻². The red circles and blue diamonds represent the conductive flux measured at T_ave = 1.0 × 10⁶ K and averaged over 10 Mm ≤ s ≤ 15 Mm, respectively.The black solid lines represent the power-law fittings to the blue diamonds. The coronal conductive flux increases with the loop length. However, the trend deviates from a simple power law. When the loop length is sufficiently small ($l_{\text{loop}}^{\text{eff}}\lesssim {10}^{1.5}\text{Mm}$) or large ($l_{\text{loop}}^{\text{eff}}\gtrsim {10}^{2.5}\text{Mm}$), the conductive flux weakly depends on the loop length. The energy flux increases around $l_{\text{loop}}^{\text{eff}}~{10}^{2.0}\text{Mm}$, with the minimum and maximum values differing by a factor of four.

An approximate power-law fit to the blue diamonds yields the following scaling law $$F_{\text{cnd}}\propto l_{\text{loop}}^{\text{eff}}{}^{0.48},$$(72)

which, in combination with the RTV scaling law, predicts $$T_{\text{top}}^{\text{RTV}}\propto {\left(F_{\text{cor}}{l}_{\text{loop}}^{\text{eff}}\right)}^{2/7}\propto l_{\text{loop}}^{\text{eff}}{}^{0.42}.$$(73)

The simulated dependence of T_top on $l_{\text{loop}}^{\text{eff}}$, Eq. (66), is reproduced by the semi-analytical arguments. Thus, the results shall be explained semi-analytically once the theoretical behaviour ofF_cor (or equivalently F_cnd) has been derived. In Sect. 4, we propose a simple model to produce F_cor.

The enhanced energy injection to the corona produces enhanced XUV emissions. Figure 8 depicts the loop-length dependence of the predicted XUV luminosity. For each coronal loop calculation, the XUV spectrum is calculated as in Fig. 5 and converted to XUV luminosity through Eq. (62). Both L_X and L_EUV exhibit increasing trends as those found in the coronal energy flux F_cnd. In particular, L_X exhibits thistrend significantly. The inferred power–law relations are $$L_{\text{EUV}}\propto l_{\text{loop}}^{\text{eff}}{}^{0.29},$$(74) $$L_{\text{X}}\propto l_{\text{loop}}^{\text{eff}}{}^{0.65}.$$(75)

The dependence of L_X on $l_{\text{loop}}^{\text{eff}}$ is further explained semi-analytically. The X-ray luminosity should be proportional to the emission measure of the corona, that is, $$L_{\text{X}}\propto n_{\text{cor}}^2{l}_{\text{loop}}^{\text{eff}},$$(76)

where n_cor is the typical number density of the coronal loop that lies between n_base and n_top. Assuming that n_cor follows the prediction of the RTV scaling law, $$n_{\text{cor}}\propto T_{\text{top}}^{\text{RTV}}{}^2/l_{\text{loop}}^{\text{eff}}\propto F_{\text{cor}}{}^{4/7}{l}_{\text{loop}}^{\text{eff}}{}^{-3/7},$$(77)

the X-ray luminosity is connected to the coronal energy flux F_cor and loop length $l_{\text{loop}}^{\text{eff}}$ by $$L_{\text{X}}\propto F_{\text{cor}}{}^{8/7}{l}_{\text{loop}}^{\text{eff}}{}^{1/7}\propto l_{\text{loop}}^{\text{eff}}{}^{0.69},$$(78)

where the simulated scaling law Eq. (72) is used in the second proportional relation. The actual scaling relation Eq. (75) is in good agreement with the semi-analytical prediction Eq. (78), validating the RTV scaling law in describing the stellar coronal properties.

We note that EUV luminosity L_EUV is weakly dependent on the loop length compared with L_X because a portion of the EUV photons originates from the upper chromosphere and transition region, which are barely affected by the variation in the coronal loop length. Within the investigated loop-length range, the ratio L_EUV∕L_X decreases from ~10 to ~ 3 as the effective loop length increases.

Fig. 7 Relation between the coronal conductive flux and loop length. Top: half loop length versus coronal conductive flux measured at Tave = 106 K (red circles) and averaged in 10 Mm ≤ s ≤ 15 Mm (blue diamonds). Also shown, by the solid line, is the power-law fitting to the blue diamonds. Bottom: same as the top panel but with respect to the effective half loop length $l_{\text{loop}}^{\text{eff}}$.

Fig. 8 (Effective) loop length versus X-ray luminosity LX (red circles) and EUV luminosity LEUV (blue diamonds). The solid and dashed lines show the power-law fittings to the simulation results. The top and bottom panels show the relation with respect to the half loop length lloop and the effective hal loop length $l_{\text{loop}}^{\text{eff}}$ (see Eq. (63)), respectively.

Fig. 9 Magnetic filling factor on the photosphere f* versus time-averaged loop-top temperature (Ttop, top panel) and loop-top electron density (ne,top, bottom panel). The half loop length is fixed to lloop = 20 Mm. Red circles show the simulation results and the solid lines show the power-law fittings to them.

3.4 Filling-factor dependence: density and temperature

The magnetic filling factor of the photosphere appears to have scaled with the Rossby number, or equivalently the magnetic activity level (Saar 1996, 2001; Reiners et al. 2009). Hence, the dependence of the coronal properties on the filling factoris worth investigating as a proxy of the dependence of the stellar corona on the activity level. The simulation results with various magnetic filling factors shall be discussed below, with the half-loop length fixed to l_loop = 20 Mm.

Figure 9 shows the filling-factor dependence of the time-averaged loop-top temperature T_top and electron density n_e,top. The coronal density and temperature increase with an increase in the magnetic filling factor following a power–law relationship $$T_{\text{top}}\propto f_{\ast }^{0.29},$$(79) $$n_{\text{e},\text{top}}\propto f_{\ast }^{0.49}.$$(80)

Although not explicitly demonstrated, the effective half-loop length also weakly depends on f_* as $$l_{\text{loop}}^{\text{eff}}\propto f_{\ast }^{0.04}.$$(81)

We note that the RTV scaling law semi-analytically predicts a relation of $$T_{\text{top}}^{\text{RTV}}\propto \sqrt{n_{\text{e},\text{top}}{l}_{\text{loop}}^{\text{eff}}}\propto f_{\ast }^{0.27},$$(82)

which is close to the simulation result Eq. (79). Alternatively, in terms of the heating rate per unit volume Q, the RTV scaling law coupled with Eqs. (79) and (81) yields $$Q\propto f^{0.935},$$(83)

while the RTV scaling law Eqs. (80) and (81) yield $$Q\propto f^{0.868}.$$(84)

The two trends are mostly consistent with each other.

The RTV scaling law predicts that the loop-top temperature T_top depends on the coronal energy flux F_cor and effective half-loop length $l_{\text{loop}}^{\text{eff}}$ as Eq. (69). A comparison between Eqs. (69) and (79) reveals the mostly proportional relationship between F_cor and f_*: $$F_{\text{cor}}{l}_{\text{loop}}^{\text{eff}}\propto f_{\ast },$$(85)

which shall be further investigated in the following section.

Fig. 10 Relation between filling factor and coronal conductive flux. Top: magnetic filling factor on the photosphere f* versus backward conductive flux measured in the corona. The half loop length is fixed to lloop =20 Mm. The red circles and blue diamonds represent the flux measured at Tave = 106 K and averaged in 10 Mm ≤ s ≤ 15 Mm. The black solid line is the power-law fitting to the blue diamonds. Bottom: same as the top panel but the vertical axis denotes $F_{\text{cnd}}{l}_{\text{loop}}^{\text{eff}}$.

3.5 Filling-factor dependence: energy flux and XUV emission

The energetics of the corona with various magnetic filling factors shall be discussed here. The top panel of Fig. 10 shows the relation between the filling factor f_* and the coronal backward conductive flux measured at T_ave = 10⁶ K (red circles) and averaged over 10 Mm ≤ s ≤ 15 Mm (blue diamonds). The two conductive fluxes exhibit different behaviours, with the red circles appearing to saturate within the large f_* limit. This difference may arise from the different physical meanings of T_ave = 10⁶ K, where the values of the red circles are measured. T_ave = 10⁶ K is located near the loop top for the small f_* cases and near the transition region for the large ones. Therefore, the blue diamonds may be interpreted as the coronal energy flux.

Given the dependence of the loop-top temperature and RTV scaling law on the filling factor, the relation among the coronal energy flux F_cor, effective loop length $l_{\text{loop}}^{\text{eff}}$, and filling factor f_* is predicted by Eq. (85). This prediction is confirmed in the bottom panel of Fig. 10, which shows the following scaling relation: $$F_{\text{cor}}{l}_{\text{loop}}^{\text{eff}}\propto f_{\ast }^{0.95}.$$(86)

This consistency indicates that the RTV scaling law must be applicable to stellar coronae with varying activity levels. The simulation results validate the previous studies on stellar coronal energetics based on the RTV scaling law (Shibata & Yokoyama 2002; Takasao et al. 2020).

Figure 11 shows the spectra of XUV emission for different magnetic filling factors. For all wavelengths below the Lyman edge (912 Å), the XUV emissions increase as the magnetic filling factor increases. This is a natural consequence of the increased coronal density shown in Fig. 9. In particular, as seen in the central panels of Fig. 11, the enhancement of the energy flux from f_* = 0.01 to f_* = 1 is centred on ~ 10 with a variance of factor ~5 in the EUV range (100 Å < λ ≤ 912 Å), whereas a significantly larger enhancement occurs in the X-ray range (5 Å < λ ≤ 100 Å). This is because of the sensitivity of the X-ray emissivity to the temperature. Also shown in the right panels of Fig. 11 are the spectra of the total photon number emitted from the star per unit time (photon number luminosity) given by 4πr²λF_λ∕(hc). For all f_*, the spectrum of the photon number luminosity is flatter than the energy spectrum, implying that the photon population with respect to wavelength is approximately uniform in the XUV range.

The magnetic-flux dependences of the X-ray luminosity L_X and EUV luminosity L_EUV are also investigated. Figure 12 reveals the variation in L_X (blue diamonds) and L_EUV (red circles) with the unsigned magnetic flux ${\Phi }_{\text{unsign}}^{\text{mag}}$ defined by $${\Phi }_{\text{unsign}}^{\text{mag}}=4\pi R_{\odot }^2{B}_{\text{s},\ast }{f}_{\ast }.$$(87)

Because the photospheric magnetic field is fixed, the filling factor f_* is the only variable in the definition of ${\Phi }_{\text{unsign}}^{\text{mag}}$. We find thatboth L_X and L_EUV follow power laws with respect to ${\Phi }_{\text{unsign}}^{\text{mag}}$ and are expressed as follows $$L_{\text{EUV}}\propto {\Phi }_{\text{unsign}}^{\text{mag}}{}^{0.81},$$(88) $$L_{\text{X}}\propto {\Phi }_{\text{unsign}}^{\text{mag}}{}^{1.19},$$(89)

which are represented by the solid and dashed lines in Fig. 12, respectively. A scaling law is observationally discovered between the unsigned magnetic flux and the X-ray luminosity (Pevtsov et al. 2003) $$L_{\text{X}}^{\text{obs}}\propto {\Phi }_{\text{unsign}}^{\text{mag}}{}^{1.15},$$(90)

which is shown by the dotted line in Fig. 12. Although the magnitudes of the X-ray luminosity from the simulation are smaller than the observed values, possibly because the proposed model reproduces the coronal properties in the activity minimum (a quiet, quasi-steady corona without flaring activities), the power-law index from the simulation (1.19) is nearly identical to that from the observation (1.15). The reproduction of the observational ${\Phi }_{\text{unsign}}^{\text{mag}}$–L_X relation provides another support to our model. The proposed ${\Phi }_{\text{unsign}}^{\text{mag}}$–L_EUV relationship, Eq. (88), can be tested by the Sun-as-a-star observations in the EUV range (e.g. Toriumi et al. 2020).

The filling-factor dependence of the EUV photon number luminosity ${\Phi }_{\text{photon}}^{\text{EUV}}$ is also investigated. ${\Phi }_{\text{photon}}^{\text{EUV}}$ is defined as the total number of EUV photons emitted from the star, that is, $${\Phi }_{\text{photon}}^{\text{EUV}}=4\pi r^2{\displaystyle {\int }_{100\AA }^{912\AA }\text{d}}\lambda \frac{\lambda F_{\lambda }}{hc},$$(91)

which acts as a proxy of the photo-ionised hydrogen number per unit time. Figure 13 shows the EUV photon number luminosity with varying magnetic filling factors. As with L_EUV, ${\Phi }_{\text{photon}}^{\text{EUV}}$ also follows a power law with respect to f_*, which is given by $${\Phi }_{\text{photon}}^{\text{EUV}}\propto f_{\ast }^{0.78}\propto {\Phi }_{\text{unsign}}^{\text{mag}}{}^{0.78}.$$(92)

A comparison of Eqs. (88) and (92), reveals that L_EUV and ${\Phi }_{\text{photon}}^{\text{EUV}}$ follow marginally different power laws with respect to f_* (or equivalently ${\Phi }_{\text{unsign}}^{\text{mag}}$).

Fig. 11 Filling-factor dependence of the XUV spectrum. Left: spectral energy flux at 1 au for f* = 0.01 (orange, solid), f* = 0.1 (green, dashed), and f* = 1 (blue, dotted). Centre: energy-flux ratio of f* = 1 to f* =0.01. Right: photon-number luminosity spectrum for f* = 0.01 (orange, solid), f* = 0.1 (green, dashed), and f* = 1 (blue, dotted). Bottom panels are the same as top panels but with logarithmic x axis thatemphasise the short-wavelength region. For better visualisation, spectra are averaged over 20 Å-bins in the top panels.

3.6 L_X –L_EUV and L_X –${\Phi }_{photon}^{EUV}$ relations

One of the objectives of this study is to propose a method to estimate the stellar EUV luminosity L_EUV from other observable quantities. Because the X-ray luminosity L_X is measurable and should be correlated with L_EUV (see e.g. Sanz-Forcada et al. 2011), the relationship between L_X and L_EUV over a wide range of filling factors and loop lengths warrants investigation.

The top panel of Fig. 14 highlights the correlation between X-ray luminosity L_X and EUV luminosity L_EUV. The black circles correspond to the simulation runs listed in Table 2. Although the simulation runs are spread over wide ranges in the filling factor and loop length, a single power-law fit adequately applies to the L_X –L_EUV relation, which is given by $$\log L_{\text{EUV}}=9.93+0.67\log L_{\text{X}},$$(93)

where L_EUV and L_X are measured in the unit of erg s⁻¹. The red points with the vertical error bars indicate the observational estimation by Sanz-Forcada et al. (2011). The pure theoretical estimation of L_EUV in this study conforms with the pure observational estimation of the same by Sanz-Forcada et al. (2011). Although the power-law index of L_X –L_EUV relation is 0.86, according to Sanz-Forcada et al. (2011), the coefficient 0.67 in Eq. (93) is, in fact, close to the observational value (0.681) derived from the Extreme Ultraviolet Explorer (EUVE) observations (Johnstone et al. 2021). The fact that the independent methods produces a similar L_X –L_EUV relation validates the proposed model of stellar coronae and their XUV emissions.

The correlation between the X-ray luminosity L_X and the EUV photon-number luminosity ${\Phi }_{\text{photon}}^{\text{EUV}}$ is also highlighted in the bottom panel of Fig. 14. They, too, follow a power–law relationship given by $$\log {\Phi }_{\text{photon}}^{\text{EUV}}=20.40+0.66\log L_{\text{X}},$$(94)

where ${\Phi }_{\text{photon}}^{\text{EUV}}$ is measured in unit of s⁻¹.

Fig. 12 Unsigned magnetic flux ${\Phi }_{\text{unsign}}^{\text{mag}}$ versus X-rayluminosity LX (blue diamonds) and EUV luminosity LEUV (red circles). Results for the fixed half loop length lloop = 20 Mm are shown. Solid and dashed lines are the power-law fittings to the simulation results. Also shown, by the dotted line, is the observational relation by Pevtsov et al. (2003).

Fig. 13 Magnetic filling factor on the photosphere f* versus photon number luminosity in the EUV range. Each simulation run corresponds to each red circle. The black solid line shows the power-law fitting to the simulation results.

4 Analytical arguments on coronal energy flux

Several scaling relations with respect to the loop length and magnetic filling factor are provided in Sect. 3. In this section, we analytically describe the coronal energy injection process governing the scaling relations. Because energy is mainly transported as transverse fluctuations of the velocity and magnetic field (Alfvén waves), the efficiency of Alfvén-wave transmission into the corona is key to the analysis. This section is devoted to reviewing the linear theory of Alfvén-wave propagation in a stellar atmosphere and to provide analytical arguments for the origin of the scaling laws discovered through the simulation.

Fig. 14 Inferred relations between X-ray luminosity and EUV emission. Top: correlation between the X-ray luminosity LX and EUV luminosity LEUV. X-ray and EUV are defined in terms of wavelength λ as 5 Å < λ ≤ 100 Å and 100 Å < λ ≤ 912 Å, respectively. Each black circle corresponds to the simulation run listed in Table 2, and the black solid line shows the power-law fitting to them. Also shown, by the red points with vertical error bars, are the empirical estimation by Sanz-Forcada et al. (2011). Bottom: same as the top panel but for the X-ray luminosity LX and photon number luminosity ${\Phi }_{\text{photon}}^{\text{EUV}}$.

4.1 Linear theory of Alfvén-wave propagation and reflection

A significant amount of Alfvén waves excited at the photosphere is reflected before reaching the solar corona (Cranmer & van Ballegooijen 2005; Verdini & Velli 2007). The relationship between the Alfvén-wave transmission efficiency and the loop properties in the linear regime remains unclear and shall be discussed in this section.

Figure 15 illustrates the schematic of the concerned system. The convection–magnetic field interaction excites the Alfvén waves propagating upwards along the expanding flux tube. In our simulation setting, the Alfvén speed is maintained at a nearly constant value during the tube expansion (Eq. (4)). Because the Alfvén-wave reflection is triggered by a gradient in the Alfvén speed (Stein 1971; An et al. 1989; Vasquez 1990), the Alfvén waves reflected in the flux-tube expansion region are feeble. Therefore, the bulk of the reflection occurs in a layer between the flux-tube merging height (where the flux-tube expansion ends) and the coronal base, which is in between the two vertical dashed lines in Fig. 15.

The Alfvén waves travelling through the reflection layer is considered in the following analysis. For simplicity, the thickness of the layer is disregarded. This approximation is validated by the fact that the reflection-layer thickness (~ 10⁰ Mm) is notably smaller than the typical wavelength of Alfvén waves therein (~ 10¹⁻² Mm). Then, the problem is simplified as a linear-wave transmission through a discontinuity, as described in the lower part of Fig. 15. Hereinafter, we denote the variables in the region on the left/right of the discontinuity as region I/II, respectively.

The equation describing the Alfvén-wave propagation is $$\left[\frac{{\partial }^2}{\partial t^2}-v_{A,\text{I/II}}^2\frac{{\partial }^2}{\partial s^2}\right]\delta v_{\text{I/II}}=0,$$(95)

the solution of which is given by $$\delta v_{\text{I/II}}=a_{\text{I/II}}{e}^{i\left(-k_{\text{I/II}}s+\omega t\right)}+b_{\text{I/II}}{e}^{i\left(k_{\text{I/II}}s+\omega t\right)}.$$(96)

We note that a_I/II (b_I/II) denotes the amplitudes of the upward (downward) propagating Alfvén wave in region I/II, respectively. Given the frequency ω, the dispersion relation yields the wavenumber k_I/II in each region. $$k_{\text{I/II}}=\omega /v_{\text{A},\text{I/II}}.$$(97)

The corresponding magnetic-field fluctuation is deduced from the linearised induction equation, that is, $$\frac{\partial }{\partial t}\delta B_{\text{I/II}}=-B_0\frac{\partial }{\partial s}\delta v_{\text{I/II}}.$$(98)

Because the fluctuations of the velocity and magnetic field are required to be continuous across the boundary, the amplitudes must satisfythe following relations. $$a_{\text{I}}+b_{\text{I}}=a_{\text{II}}+b_{\text{II}},$$(99) $$-k_{\text{I}}{a}_{\text{I}}+k_{\text{I}}{b}_{\text{I}}=-k_{\text{II}}{a}_{\text{II}}+k_{\text{II}}{b}_{\text{II}}.$$(100)

Solving these two equations yields $$a_{\text{I}}=\frac{k_{\text{I}}+k_{\text{II}}}{2k_{\text{I}}}{a}_{\text{II}}+\frac{k_{\text{I}}-k_{\text{II}}}{2k_{\text{I}}}{b}_{\text{II}}.$$(101)

The energy transmission rate T_AW is defined asthe ratio of the forward-propagating Alfvén-wave energy in regions I and II. Namely, $$T_{\text{AW}}=\frac{{\rho }_{\text{II}}{v}_{\text{A},\text{II}}{\left|a_{\text{II}}\right|}^2}{{\rho }_{\text{I}}{v}_{\text{A},\text{I}}{\left|a_{\text{I}}\right|}^2}.$$(102)

We obtain a general form of the energy transmission rate within the limit of ρ_II ∕ρ_I ≪ 1 by combining Eqs. (101) and (102): $$T_{\text{AW}}=\frac{4v_{\text{A},\text{I}}{\left|a_{\text{II}}\right|}^2}{v_{\text{A},\text{II}}\left({\left|a_{\text{II}}\right|}^2+a_{\text{II}}{b}_{\text{II}}^{\ast }+a_{\text{II}}^{\ast }{b}_{\text{II}}+{\left|b_{\text{II}}\right|}^2\right)},$$(103)