© ESO 2020

1. Introduction

Recent high-cadence and high-resolution observations have established the presence of waves and oscillations throughout the solar atmosphere (see e.g. reviews by Nakariakov & Verwichte 2005, Banerjee et al. 2007; Zaqarashvili & Erdélyi 2009; Arregui et al. 2012; De Moortel & Nakariakov 2012, Mathioudakis et al. 2013, Arregui 2015, Jess et al. 2015). Most of these perturbations have been interpreted as magnetohydrodynamic (MHD) waves and in many cases are reported to contain a substantial amount of energy, leading to a renewed interest in MHD wave dissipation as a potential coronal heating mechanism (e.g. Parnell & De Moortel 2012; Arregui 2015).

Phase mixing of Alfvén waves (Heyvaerts & Priest 1983) is one of the mechanisms proposed to address the slow dissipation of wave energy with classical transport coefficients; in the presence of a cross-field gradient in the Alfvén speed, waves on neighbouring field lines move out of phase, building up large gradients which, in turn, lead to enhanced dissipation.

As heating occurs in the corona, an enhanced conductive flux is driven downward from the corona to the lower atmosphere which, in turn, leads to an upward enthalpy flux (evaporation), increasing the coronal density. Unless heating continues, the corona will then start to cool whilst the density continues to increase until a balance between thermal conduction and radiation is reached. In the final stage, as the corona continues to cool, mass starts to drain from the loop. (see e.g. Cargill 1994; Klimchuk 2006). Hansteen et al. (2010), for example, studied this mass and energy cycle between the lower atmosphere and the corona and found that both downflows and upflows are present at locations of (strong) magnetic field braiding, leading to redshifts and blueshifts of the order of about 5 km s⁻¹ (see also e.g. Zacharias et al. 2011; Guerreiro et al. 2013). Observationally, a multitude of studies have found evidence for the presence of a complex interplay between upflows and downflows in the Transition Region (TR) and lower corona (e.g. Del Zanna 2008; Feldman et al. 2011; Dadashi et al. 2011, 2012; Tripathi et al. 2012a,b; McIntosh et al. 2012; Winebarger et al. 2013). A comprehensive review of the modelling of impulsive heating in coronal loops can be found in e.g. Reale (2014). Bradshaw & Cargill (2013) presented a detailed discussion on the importance of numerical resolution in the TR to accurately model the mass and energy exchange. Various solutions have been proposed to address this problem such as (controlled) broadening of the TR (e.g. Lionello et al. 2009; Mikić et al. 2013; Johnston & Bradshaw 2019; Johnston et al. 2020) and a jump-condition approach (Johnston et al. 2017a,b, 2019).

Although phase mixing of Alfvén waves has been studied extensively as a potential coronal heating mechanism, it has usually been investigated without taking into account the mass exchange with the lower atmosphere described above. In particular, heating could change the local density profile (e.g. through evaporation or draining) and, hence, the Alfvén speed profile which, in turn, could alter the phase mixing process. This feedback mechanism was investigated for resonant absorption by Ofman et al. (1998). These authors used thermodynamic equilibrium scaling laws to instantaneously adjust the local plasma parameters and found that the feedback caused the heating layers to drift. However, Cargill et al. (2016) show that the “detuning” of the local resonance conditions is a relatively slow process and hence that the instantaneous update to the local plasma conditions used by Ofman et al. (1998) is likely to overestimate the effect of the feedback mechanism.

Cargill et al. (2016) also demonstrated conceptually that heating by phase mixing of Alfvén waves cannot sustain the (assumed) background density profile as heating takes place preferentially where the Alfvén speed (or density) gradient is largest, but cooling and draining occur fastest where the density is higher. In other words, the spatial profile of the heating by phase mixing is not compatible with the imposed density structuring. These authors also found that although the heating feedback mechanism could lead to a highly structured density profile, this only happened on long timescales, comparable to the loop cooling and draining timescales and, hence, the changes in the density profile are dominated by the thermal evolution.

In this paper, we present 2D MHD simulations of phase mixing of Alfvén waves in a coronal loop, including optically thin radiation and thermal conduction, and we study the plasma heating and subsequent evaporation of mass from the lower atmosphere. This allows us to investigate whether the dissipation of phase-mixed Alfvén waves does indeed lead to evaporation from the lower atmosphere and whether such mass flows in turn affect the phase mixing process. This paper is structured as follows. In Sect. 2, we introduce the methods and numerical setup for the coronal loop model used in this study. The results are presented in Sect. 3, followed by the discussion and conclusions in Sect. 4.

2. Model and numerical method

2.1. Equations and initial conditions

We use the 2D Lagrangian Remap code Lare2D (Arber et al. 2001) to solve the (normalised) MHD equations:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}=-\mathrm{\nabla }·\left(\rho \mathit{v}\right),\hfill \end{array}$$(1)

$$\begin{array}{cc}& \rho \frac{D\mathit{v}}{\mathit{Dt}}=-\mathrm{\nabla }P-\rho \mathit{g}+\mathit{j}\times \mathit{B}+\rho \nu \left({\mathrm{\nabla }}^2\mathit{v}\right),\hfill \end{array}$$(2)

$$\begin{array}{cc}& \rho \frac{Dϵ}{\mathit{Dt}}=-P\mathrm{\nabla }·\mathit{v}-\mathrm{\nabla }·{\mathit{F}}_{\mathbf{c}}-{\rho }^2\mathrm{\Lambda }(T)+H_{\mathrm{bg}}(x)\hfill \end{array}$$(3)

$$\begin{array}{cc}& +\eta j^2+Q_{\mathrm{visc}},\hfill \end{array}$$(4)

$$\begin{array}{cc}& P=2\rho T,\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{\partial \mathbf{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{v}\times \mathit{B})-\mathrm{\nabla }\times (\eta \mathrm{\nabla }\times \mathit{B}),\hfill \end{array}$$(6)

$$\begin{array}{cc}& \mathit{j}=\mathrm{\nabla }\times \mathit{B}.\hfill \end{array}$$(7)

Here, D/Dt is the convective derivative, ρ is the mass density, v the velocity, P the gas pressure, B the magnetic field, j the current density, η the resistivity, ρν the dynamic viscosity and $ϵ=\frac{P}{(\gamma -1)\rho }$ the specific internal energy. We include a field-aligned gravitational acceleration g_||, where we use a semi-circular loop profile (e.g. Reale 2010). Equation (4) is the energy equation where we include thermal conduction, optically thin radiation, a background heating function and the viscous heating term, Q_visc = ρν∇v : (∇v)^T. The conductive flux is given by F_c = −κ₀T^5/2(B ⋅ ∇T)B/B² and Λ(T) = χT^α is the optically thin radiative loss function, where χ and α are described in Klimchuk et al. (2008).

To construct a 2D loop model with a cross-field density gradient, we combine a field-aligned hydrostatic (thermodynamic) equilibrium (e.g. Reale 2010; Johnston et al. 2017a) with a background heating function that has a transverse (cross-field) profile. In our 2D Cartesian reference frame, y is the field-aligned direction, x the cross-field direction and z is the invariant direction. The non-uniform profile for the background heating function is given by

$$\begin{array}{c}\hfill H_{\mathrm{bg}}(x)=\frac{H_1+H_2}{2}+\frac{H_2-H_1}{2}tanh(a(-|x|+1)),\end{array}$$(8)

where a = 5, H₂ = 4H₁ and H₁ = 3.24 × 10⁻⁶ J m⁻³ s⁻¹ is the (background) heating required for a field-aligned (1D) equilibrium in the external region. After (numerical) relaxation, this combination of parameters results in a density ratio of ρ_i/ρ_e  =  2.4 at the loop apex, where ρ_i and ρ_e refer to the interior and exterior densities, respectively. The transverse profile of the heating function in Eq. (8) has been chosen as a model-representation of a concentrated background heating. Because of the transverse structuring of the heating function, a transverse density profile is created in the corona (a coronal “loop”) which will lead to phase mixing.

The numerical domain is 4 Mm × 120 Mm and consists of 256 gridpoints in x (cross-field direction) and 4096 gridpoints in y (field-aligned direction). The boundary conditions are periodic in x and zero gradient in y with the velocity set to zero. The setup includes an isothermal model-chromosphere (T = 2 × 10⁴ K), at both ends of the loop, extending over the first and last 8 Mm of the domain. This model-chromosphere acts as a simple mass reservoir and does not include detailed chromospheric physics such as partial ionisation.

In order to fully resolve the TR, we use the approach proposed in Lionello et al. (2009) and Mikić et al. (2013): below a fixed cut-off temperature T_c, the optically thin radiative loss rate is decreased and the parallel thermal conductivity is increased:

$$\begin{array}{c}\hfill {\kappa }_{||}(T)=\{\begin{array}{cc}{\kappa }_0{T}^{5/2}\hfill & T\ge T_{\mathrm{c}}\hfill \\ {\kappa }_0{T}_{\mathrm{c}}^{5/2}\hfill & T<T_{\mathrm{c}},\hfill \end{array}\end{array}$$(9)

and

$$\begin{array}{c}\hfill \mathrm{\Lambda }(T)=\{\begin{array}{cc}\mathrm{\Lambda }(T)\hfill & T\ge T_{\mathrm{c}}\hfill \\ \mathrm{\Lambda }(T){\left(\frac{T}{T_{\mathrm{c}}}\right)}^{5/2}\hfill & T<T_{\mathrm{c}}.\hfill \end{array}\end{array}$$(10)

This approach increases the field-aligned temperature length scale $L_{\mathrm{T}}=T/\left|\frac{\partial T}{\partial y}\right|$, thus broadening the TR. In the simulations presented in this paper, we have set T_c = 5 × 10⁵ K.

Figure 1 shows contour plots of the temperature and density which are obtained after the numerical relaxation. In the coronal part of the loop, there is a clear cross-field gradient in the temperature and the density (as well as in the Alfvén speed, see Fig. 3), with the coronal temperature ranging between 1.0−1.6 MK and the coronal density varying between ρ = 2.5−6 × 10⁻¹³ kg m⁻³ (giving a density ratio of ρ_i/ρ_e = 2.4 at the apex). This density range is lower than values usually measured in Active Region loops but can be seen as representative of the Quiet Sun (see e.g. Brooks et al. 2009, 2012; Brooks 2019). The key aspect for phase mixing is the density contrast which, although modest, is reasonable for Quiet Sun loops.

Fig. 1. Contour plots of the temperature (left) and density (right) after the numerical relaxation.

Overplotted (black lines) on these contour plots are the (initial) locations of the chromosphere-TR and TR-corona boundaries. Here, we have defined the top of the TR to be the location where thermal conduction (∇⋅F_c) changes sign (i.e. changes from an energy loss to an energy gain). The chromospheric boundary is defined as the location where the temperature has decreased to T = 2 × 10⁴ K. The locations of these boundaries are determined in the initial setup and are then assumed fixed in all subsequent calculations. For the particular setup studied in this paper, this assumption is reasonable as the modest heating occurring in the corona does not affect the location of the TR and the chromospheric boundaries. Figure 2 shows a (field-aligned) cross-section of the temperature and density at x = −2 Mm (exterior), x = −1 Mm (shell) and x = 0 Mm (middle of the loop). It is clear from this figure that the temperature and density are well resolved in the transition region. Indeed, the minimum temperature length scale L_T in the loop is of the order of 125 km, occurring in the lower TR. Thus, the field-aligned grid resolution (dy ∼ 30 km) is more than adequate to fully resolve the TR. We note here that the temperature at the top of the TRs remains between 0.51−0.83 MK and, hence, the cutoff temperature T_c = 0.5 MK is always located in the TR.

Fig. 2. Field-aligned temperature (left) and density (right) at x = −2 Mm (black), x = −1 Mm (blue) and x = 0 Mm (green). The small panel inserted on the temperature plot shows a zoomed version in the TR (between y = 7.5 Mm and y = 12 Mm). The symbols represent the numerical gridpoints.

Figure 3 shows contours of the magnetic field components B_x and B_y as well as the Alfvén speed. Before relaxation, the background field is set to be uniform in y with B_y = 10 G. Following the (numerical) relaxation, there is a shallow cross-field gradient in B_y as well as slight expansion of the field in the TR (reflected in the appearance of a small B_x component). The Alfvén speed shows a similar pattern as the contour of the density in Fig. 1, varying between 1000−1800 km s⁻¹ in the corona, where the highest Alfvén speed occurs outside the coronal part of the loop. Finally, we note that our 2D setup results in a plasma beta of the order of 0.01 in the corona.

Fig. 3. Contour plots of the magnetic field components Bx (left) and By (middle) and the Alfvén speed (right) after the numerical relaxation.

2.2. Driver

We implement a sinusoidal driver at the top of the first chromosphere (y = 7.8 Mm) in the invariant z direction. The driver is implemented as an external force, F_z, in the z component of the momentum equation (Eq. (2)) as

$$\begin{array}{c}\hfill F_z=-\rho v_0\omega cos(\omega t).\end{array}$$(11)

Here v₀ = 0.7 km s⁻¹ is 1% of the local Alfvén speed, and $\omega =\frac{2\pi }{P}$ is the angular frequency, where P ≈ 12 s is the period of the driver. This driver generates left and right (up and down) propagating Alfvén waves that propagate along the field (in the y direction). Although most observed waves and oscillations in the solar corona have periods of a few minutes (see e.g. De Moortel & Nakariakov 2012), in this study, we use a high frequency driver to ensure a sufficient number of wavelengths in the corona to lead to significant phase mixing (see Sect. 4 for further discussion).

For the parameters considered in our initial conditions, the thermal conduction and optically thin radiation timescales are of the order of 1000–4000 s. Therefore, to ensure the effects of thermal conduction and optically thin radiation are captured, we continue to drive the simulation for t = 500 P ≈ 6000 s.

3. Results

In this paper, we focus on the analysis of a Lare2D simulation which includes a dynamic viscosity of ρν = 5 × 10⁻⁴ kg m⁻¹s⁻¹ but without resistivity (i.e. η = 0), to avoid diffusion of the background Alfvén speed profile. In addition, we will also present results from the corresponding ideal simulation, as well as a simulation which was run for the same total time but without a driver implemented (i.e. a simulation where the numerical relaxation after imposing the cross-field heating profile continues for the total time interval considered here). In Fig. 4 we show contours of v_z for the viscous simulation at t = 93 s (left), when the first pulse generated by the driver reaches the far TR, and t = 983 s (right), representing a later stage of the simulation. Initially, after a small amount of phase mixing in the lower part of the loop, the Alfvén speed profile reverses in the corona and hence this phase mixing is “undone” before more substantial phase mixing occurs close to the far leg of the loop (i.e. the loop leg furthest away from the driver). Due to the high frequency of the driver, the wavelength of the waves is relatively small in the corona, which allows most of the waves to be transmitted into the transition region. Only a modest amount of wave energy (of the order of 15 %) is reflected back into the coronal part of the domain. Figure 4 (right) shows a contour of v_z at t = 983 s where the oscillation amplitudes are now of the order of 7 km s⁻¹ (compared to 2–3 km s⁻¹ at earlier times). The interference of the driven and reflected waves leads to the partial establishment of a standing wave in the coronal part of the domain. When dissipation and phase mixing are absent (the core part of the loop in the ideal simulation), an almost exact standing wave is established after some time. Either in the shell region (where phase mixing takes place) or when viscosity is present, the coronal part of the domain eventually evolves to a steady-state which is a combination of standing and propagating waves. Phase mixing of the driven and reflected waves in the shell region leads to an in intricate pattern of fine scale structuring. In the viscous simulation, this will lead to heating in the shell regions, as can be seen in the left-hand panel of Fig. 6, which shows contours of the relative change in temperature, (T − T₀)/T₀ at the end of the simulation (t = 5734 s).

Fig. 4. Contour of vz (km s−1) for the viscous simulation at t = 93 s (left panel) and t = 983 s (right panel).

3.1. Long-period field-aligned flows and oscillations

The aim of this paper is to study any changes in the coronal density profile induced by evaporation of TR or chromospheric plasma, following heating in the corona. However, other field-aligned flows are also present in the domain. In order to be able to distinguish these flows from the evaporation, we now look at a detailed comparison between the ideal and the viscous simulations, as well as the non-driven simulation.

Figure 5 shows the field-aligned velocity v_y, integrated over the near TR-coronal boundary for the core of the loop (−1 <  x <  1 Mm) as a function of time. We can see that a long-period (P ∼ 650 s) oscillation is present in all simulations, which is due to the system still evolving under the imposed cross-field background heating profile. The period is related to the slow travel time in the domain (where the time for a slow wave to travel from the top of the first chromosphere to the top of the far chromosphere in the middle of the domain is about 650s). The oscillations are very small amplitude, of the order of |v_y| < 0.1 km s⁻¹. Comparing with the non-driven (dotted line) simulation, we can see that the maxima in the ideal simulations (dashed line) are slightly larger. At this TR-coronal boundary, positive values of the field-aligned velocity v_y correspond to an upflow into the corona. These additional upflows are due to the ponderomotive force associated with the Alfvén waves generated by the forcing term. A more detailed explanation of these ponderomotive upflows is given using a simple 1D model in Appendix A. For the viscous simulation (solid line), we find upflows which are slightly smaller than in the ideal case, due to the effect of viscosity on the upflows. Although the ponderomotive upflows are not the main focus of our study, it is worth remarking here that these upflows are likely over-estimated, as we consider the plasma to be fully ionised everywhere in the numerical domain (see e.g. Laming 2017). However, this will equally affect the ponderomotive upflows in all simulations considered here and, hence, does not affect our estimate of the evaporative upflows, which we obtain by comparing the upflows in the ideal and non-ideal simulations.

Fig. 5. Mean field-aligned velocity vy integrated over the first TR-coronal boundary in the central part of the loop (−1 <  x <  1 Mm) for the viscous (solid line), ideal (dashed line) and non-driven (dotted line) simulations as a function of time.

3.2. Plasma changes

As already mentioned earlier, phase mixing in the shell regions leads to dissipation of the wave energy in the viscous simulations (Fig. 4) and hence to an increase in the local temperature (left-hand panel of Fig. 6). At early times, the dissipation occurs in the far coronal part of the loop (z >  60 Mm) as can be expected from the phase mixing pattern in v_z (LHS of Fig. 4). At later times, when a steady state is established, a small increase in the relative temperature all along the shell regions can be observed, although the increase remains slightly higher in the far leg of the loop at all times. The actual increase in the coronal temperature is very small (of the order of 4000 K or less than 0.5% of the background coronal temperature) due to the limited energy input, which results from the combination of the relatively small amplitude of the wave driver and the limited reflection at the TR-coronal boundary. In other words, any wave energy that has not yet been dissipated by the time the propagating waves reach this boundary is mostly transmitted to the lower atmosphere rather than being contained in the coronal part of the loop. There is no noticeable change in temperature in the core (coronal) part of the loop and the loop environment where phase mixing does not occur.

Fig. 6. Contours of (T − T0)/T0 (left) and (ρ − ρ0)/ρ0 (right) at t = 5340 s for the viscous simulation.

The right-hand panel of Fig. 6 shows the relative change in the density at the end of the simulation. Again, the most significant changes in the coronal part of the loop are limited to the shell region, where an increase of about 1% in the density has occurred (with a corresponding decrease in the TR and chromosphere of the shell region). However, as described above in Sect. 3.1, evaporation due to (coronal) heating is only part of this change in density, as upflows due to the ponderomotive force associated with the Alfvén wave driver and the ongoing evolution of the background due to the enhanced heating profile are also present. Given that the upflows due to the background evolution are oscillatory (dashed line in Fig. 5), it is likely that this process does not contribute significantly to the relative density change in the corona. This is indeed confirmed in the left panel of Fig. 8 where we see that there is essentially no change in mass in the corona for the non-driven simulation (dotted black line).

In Fig. 7, we show time-distance diagrams of the change in relative temperature (top) and density (bottom) along a cross section in the middle of the shell region (x = −1.18 Mm, marked by the vertical dashed lines in Fig. 6). A few features immediately stand out in the evolution of the relative temperature changes. First of all, the time-distance plot of the relative temperature confirms that the heating is indeed asymmetric, with the relative increase in temperature slightly larger in the far leg of the loop (y >  60 Mm). Because of the specific Alfvén speed profile, phase mixing only develops significantly in the loop leg further away from the driver and hence the wave dissipation is more efficient in the far leg. Secondly, there are clear periodic changes in the temperature, propagating upwards from the TR into the corona with a speed of about 120 km s⁻¹ (estimated from the slope of the diagonals in the contour plot). A comparison with Fig. 5 shows that this periodicity must be associated with the ongoing evolution induced by the heating profile, as the same periods are also present in the field-aligned flows in both the non-driven and ideal simulations. As the perturbations associated with both this ongoing relaxation and the ponderomotive force are adiabatic, no net heating (i.e. integrated over time) is expected. Indeed, this is confirmed by analysis of the corresponding time-distance plot for the ideal simulation (not shown), where periodic changes in the temperature are present in the TRs but there is no net increase of temperature over time. For the viscous simulation on the other hand, there is a clear net increase in the temperature over time in both the coronal part of the domain and the lower atmosphere.

Fig. 7. Time distance graphs of the relative changes in temperature (top) and density (bottom) along x = −1.18 Mm (marked by vertical dashed line in Fig. 6).

Focusing on the first 1000 s of the relative temperature change in the viscous simulation, there is a first small heating event occurring roughly between y = 60−80 Mm at t = 500−800 s. Comparing with the contour plots of v_z (the driven Alfvén waves) in Fig. 4, it is clear that this heating is a direct result of the viscous dissipation of the phase-mixed Alfvén waves. Following this heating event, a thermal conduction front can be seen in the form of a downward propagating increase in temperature (i.e. towards higher values of y) with a steep gradient (high velocity). There is some evidence of the temperature change also extending towards the other loop leg (y <  60 Mm) although this is harder to see. Similar heating “events” then appear to re-occur at intervals of about 600–700 s. However, this periodicity is essentially a modulation of the heating from phase mixing by the long-period oscillations discussed above. Although it is not easy to see the downward propagating change in temperature due to conduction, some instances can be identified where the contours are tilted towards the right in the upper half of the domain, which corresponds to propagation from the corona to the TR in the far leg of the loop (y >  60 Mm). Finally, given the background Alfvén speed profile, phase mixing of the driven Alfvén waves will also occur in the TRs and, hence, direct heating in these regions is expected. Therefore the overall, net increase in temperature in the TRs over time is a combination of direct heating and the downward conduction of coronal changes in the temperature.

The bottom panel in Fig. 7 shows the corresponding time-distance diagram for the relative change in density. There are clear upflows present propagating from the lower parts of the loop up into the corona. Again, the periodicity of these density changes is associated with the background long-period oscillations as discussed in Sect. 3.1. The upflows are caused by two effects: the ponderomotive force associated with the driven Alfvén waves and the evaporation from the lower atmosphere into the corona following (coronal) heating by phase mixing of the Alfvén waves. Indeed, comparing with the ideal simulation (where only the ponderomotive mass increase is present) does reveal that for the viscous simulation, there is a larger increase in mass in the shell regions, which is associated with the evaporation from the lower atmosphere.

To confirm this, we look at the change in mass integrated over time in Fig. 8. The right-hand panel shows the change in mass integrated over the coronal part of the shell regions only, for the viscous (solid blue line) and the ideal (dotted-dashed blue line) simulation. An increase in mass is indeed present in both cases but the increase is very clearly larger in the viscous simulation. This excess in mass increase is associated with the evaporative upflows from the lower atmosphere induced by the Alfvén wave heating. The green lines correspond to the time integrated, averaged mass flux (where the averaging is done over the 4 TR-coronal boundaries of the shell regions) and, hence, confirm that the mass increase in the coronal part of the shell is indeed due to mass flux inflowing from the lower atmosphere.

Fig. 8. Left: mass increase (kg m−1) for the three simulations (viscous, solid lines; ideal, dashed lines; non-driven, dotted lines). The black lines correspond to the integrated mass in the coronal part of the domain whereas the green lines represent the mass in the non-coronal part. Right: mass increase (kg m−1; blue lines) but now only integrated over the coronal part of the shell region (viscous, solid lines; ideal, dotted-dashed lines). The green lines represent the time integrated, averaged mass flux (where the averaging is done over the 4 boundaries of the shells).

The left-hand panel of Fig. 8 shows the mass changes for the full loop (i.e. including the loop core and the environment outside the loop) for the three simulations with the coronal mass change plotted in black and the lower atmosphere in green. This figure essentially summarises the findings described above. Looking at the non-driven simulation (dotted line) confirms there is no significant change in mass in the corona or lower part of the domain due to the oscillatory field-aligned perturbations associated with the continued (background) evolution. For the ideal and viscous simulations, there is a clear increase of mass in the coronal part of the loop (with a corresponding decrease in the lower atmosphere). For the full domain (rather than just the shell region as in the right-hand panel), the mass increase is dominated by the ponderomotive effect, given the relatively small difference between the mass increase for the viscous (solid line) and ideal (dashed line) simulations. Saying that, the viscosity also acts on the ponderomotive upflows (see Fig. 5) and, hence, the effect of the evaporative upflows is somewhat larger than indicated by the direct difference between the two simulations.

We can estimate the magnitude of the evaporative upflows by considering the difference in the mean upflows through the TR-coronal boundaries of the shell regions between the viscous and ideal simulations. We find that this evaporative upflow is of the order of 5−20 m s⁻¹ and is approximately equal on both the near (closest to the driver) and far TR-coronal boundaries. Using the continuity equation and the observed mass increase in the shell regions confirms that this evaporative upflow does indeed match the required average evaporation.

Figure 9 shows the time and volume integrated vertical component of the Poynting flux ${\displaystyle {\int }_0^t{\displaystyle \int {\left(E\times B\right)}_y}}\text{d}S\text{d}t({\text{Jm}}^{-1}$, over the TR-coronal boundaries of one of the shell regions of the loop. The green lines represent the near (closest to the driver) TR-coronal boundary and the black lines the far boundary. Firstly, as expected for the ideal simulation, the Poynting flux entering the coronal shell region (through the near boundary - green dashed line) is the same as the Poynting flux leaving this region (black dashed line). For the viscous simulation on the other hand, we can see that the Poynting flux leaving the coronal shell region (solid black line) is substantially smaller than the Poynting flux entering the region (green dashed line). This difference is of the same order as the internal energy increase in the shell region (not shown here). The Poynting flux entering the shell region through the near TR-coronal boundary (green lines) differs in the ideal and viscous simulation due to the different mixture of propagating and standing waves in this region. For the ideal simulation, the standing wave component dominates, which reduces the amount of Poynting flux through the near boundary, as the standing regime creates a node in v_z (see e.g. Prokopyszyn et al. 2019). Comparing with coronal energy requirements, the difference between the average Poynting flux entering and leaving through the TR-coronal boundaries of the shell region for the viscous simulation is of the order of 1.5−2 Jm⁻² s⁻¹, several orders of magnitude too low, even for the Quiet Sun (∼3 × 10² Jm⁻² s⁻¹, Withbroe & Noyes 1977).

Fig. 9. Time integrated Poynting flux ${\displaystyle {\int }_0^t{\displaystyle \int {\left(E\times B\right)}_y}}\text{d}S\text{d}t(\text{J}{\text{m}}^{\text{-1}}$, with S the boundaries of the left shell, for the viscous (solid lines) and ideal simulation (dashed lines). The green lines represent the lower TR-coronal boundary and the black lines the upper boundary.

4. Discussion and conclusions

In this paper, we use 2.5D numerical simulations to investigate whether a pre-existing density profile is modified by the evaporative upflows following (viscous) heating by phase mixing of Alfvén waves. The Alfvén waves are driven through an additional forcing term in the momentum equation. The simulations account for the effects of gravitational stratification, thermal conduction, and optically thin radiation. A fully resolved atmospheric model (including the chromosphere as a mass reservoir) is achieved with relatively modest resolution by artificially broadening the TR using the approach proposed by Lionello et al. (2009) and Mikić et al. (2013). To include a transverse density profile (perpendicular to the background magnetic field) in the initial setup, a heating function which varies in the cross-field direction is imposed.

Throughout the simulations (running for a total time of more than 5000 s), a complex combination of Alfvén waves (v_z) and longitudinal perturbations (v_y) is present in the domain. The longitudinal perturbations are a combination of long-period oscillations resulting from the ongoing evolution due to the imposed cross-field heating profile, the ponderomotive effects associated with the driven Alfvén waves, and the evaporative upflows resulting from heating occurring in the coronal part of the shell regions of the loop due to phase mixing of the Alfvén waves.

By comparing with the results from the ideal simulations, we are able to distinguish the evaporative upflows present in the viscous simulations from the other (ideal) longitudinal perturbations. This allows us to identify the change in mass in the shell regions of the coronal loop which results from the viscous heating by phase mixing of Alfvén waves. For the particular setup studied in this paper, the amount of heating through viscous dissipation of the phase-mixed Alfvén waves in the corona is extremely small and, hence, the evaporative upflows associated with this heating are insignificant. Therefore, in this case, the heating-evaporation cycle has no noticeable effect on the transverse density profile (or Alfvén speed profile) which causes the Alfvén wave phase mixing.

As remarked above, the choice of our high-frequency driver (P ∼ 12 s) implies that a substantial amount of the wave energy is transmitted down to the far TR and chromosphere. Indeed, the short wavelengths associated with the high-frequency driver allow the rapid development of phase mixing in the shell regions of the loop but, at the same time, a substantial amount of energy (about 20–25%) is not dissipated by the time the waves reach the far TR and is lost to the lower atmosphere (see e.g. Hollweg 1984a,b; Berghmans & de Bruyne 1995; De Pontieu et al. 2001). This implies that simply increasing the amplitude of the driven Alfvén wave would only have a limited effect, as a substantial amount of the wave energy would still be lost from the coronal volume, even in the shell regions. In the core of the loop, almost all energy (of the order of 85%) is transmitted to the far TR and chromosphere. Although the TR has been broadened somewhat to ensure sufficient numerical resolution (see Sect. 2), this broadening does not substantially affect the amount of wave energy that is transmitted and reflected. We also note here that the method proposed by Lionello et al. (2009) and Mikić et al. (2013) conserves the total amount of energy that is delivered to the chromosphere and that although there can be small differences with the flows in the modified region (i.e. where T <  T_c), the mass flux into the corona is preserved (Johnston et al. 2020). Therefore, the artificial broadening does not have a significant effect on the evaporative upflows in this study.

In a forthcoming paper, we will extend the study presented here to the low frequency limit, as waves and oscillations observed so far in the solar corona mostly have periods on the order of a few minutes (e.g. De Moortel & Nakariakov 2012). This will lead to longer wavelengths and, hence, the waves would experience more reflection at the boundaries between the coronal part and the lower atmosphere which would, in turn, allow for more energy to be contained in the coronal part of the entire loop. In the initial phase of propagation, there would be less pronounced phase mixing for longer wavelength waves but, over time, phase mixing (in time) would lead to large gradients and, hence, dissipation of the wave energy in the coronal part of the shell regions of the loop. However, longer wavelength waves would take longer to phase mix and therefore, even if more energy is retained in the corona, it remains to be investigated whether this can result in more energy converted into heating and, if so, whether this results in stronger evaporative upflows on relevant timescales.

Even in the long wavelength limit, the simulations would still be highly idealised and mostly a theoretical study of the evaporation. Several other aspects could be included to make the model more representative of actual coronal loops. For example, in the current study we have only considered viscous heating, and neglected the effect of resistivity. As the resistivity is expected to be very low in the solar corona, this might be acceptable, but it would still be useful to examine the effect of including resistivity. In particular, it would be instructive to study how including resistivity balances the effects of changing the background Alfvén speed profile (through diffusing the background magnetic field configuration) and the stronger heating resulting from the additional resistive heating. Further possibilities include a more realistic broadband driver to establish how this affects the energy input into the corona as well as a magnetic field configuration with concentrated sources, where stronger divergence of the magnetic field could enhance the phase mixing process (e.g. De Moortel et al. 2000).

Finally, the study presented here mostly maintains the background initial conditions (i.e. the actual loop profile) through the presence of an imposed artificial background heating function. Without the presence of the background heating, Cargill et al. (2016) argued that the thermal evolution (i.e. the loop cooling) would lead to significant changes in the cross-field density profile (mostly due to draining) on timescales quicker than the heating provided by phase mixing of Alfvén waves. Although heating by phase mixing of Alfvén waves does not appear efficient for the parameters used in this current study, it does not provide a definitive answer to the question whether phase mixing of Alfvén waves can be a self-consistent, stand-alone heating mechanism of solar coronal loops due to the limitations of our model setup described above.

Acknowledgments

This work has received support from the UK Science and Technology Facilities Council (Consolidated Grant ST/K000950/1), the European Union Horizon 2020 research and innovation programme (grant agreement No. 647214) and the Research Council of Norway through its Centres of Excellence scheme, project number 262622.

References

Appendix A: Ponderomotive upflows

In this appendix, we provide a more detailed description of the upflows associated with the Alfvén wave driver. Ponderomotive upflows, resulting from the non-linear magnetic pressure force (e.g. Verwichte et al. 1999) have been discussed in the context of standing Alfvén waves in coronal loops (see e.g. Terradas & Ofman 2004; McLaughlin et al. 2011; Verwichte et al. 2017). In our model setup however, the Alfvén waves are propagating, driven by an additional forcing term in the momentum equation (see Eq. (11)).

To isolate the effects of the ponderomotive force, we consider a simple 1D, uniform model and use the same force driver as in Eq. (11), situated at an arbitrary distance of y = 200 Mm from the bottom boundary. This driver generates both left and right propagating Alfvén waves. After some time, the left-propagating waves reach the lower boundary (at y = 0 Mm) where they reflect. We find that as these reflected waves (now right-propagating) interact with the left-propagating waves generated by the driver (at y = 200 Mm), a standing mode is

generated. After t = 347 s, the right-propagating waves have reached the location of the driver again, setting up a standing mode in the entire lower part of the domain (0 <  y <  200 Mm). Hence, although the driver generates propagating waves, the reflection from the lower boundary indirectly leads to a standing wave in the lower part of the domain as the length of the lower part of the domain (i.e. the distance between the bottom boundary and the driver) was chosen to be a multiple of the wavelength of the driven waves. Figure A.1 (left) shows both the velocity (solid line) and magnetic (dot-dashed line) perturbations at the location of the driver (y = 200 Mm) with time. We see that when the reflected waves reach the driver location (at t = 347 s), this location becomes a node in v_z and an anti-node in b_z from that time onwards. The presence of this standing wave node leads to a rapid increase in the magnetic pressure force, driving an upflow (v_y) as can be seen in Fig. A.1 (right). This upflow then leads to a compression of plasma, which propagates upwards at the slow speed.

Fig. A.1. Left: vz/max(vz) and Bz/max(Bz) with time at the location of the driver (y = 200 Mm). Right: upflows vy (km s−1) as a function of time at the location of the driver.

Figure A.2 shows the relative mass integrated over the region y = 230−260 Mm, as a function of time. We can see a rapid increase in mass after t = 503 s, which accounts for the time required for the mass increase to propagate (at the slow speed) from the location of the driver (at y = 200 Mm) to the region where we are measuring the mass increase (starting at y = 230 Mm). By t = 659 s, the mass increase has reached y = 260 Mm, which marks the end of the region in which we are sampling the mass increase and from then on, the mass stays constant in this interval.

Fig. A.2. Relative mass increase ∫(ρ − ρ0)dV/∫ρ0dV with time in the region y = 230−260 Mm. The red vertical lines mark the times when the slow wave, that is generated at t = 347 s at the location of the driver, arrives and leaves this region.

In the example described above, the lower boundary was set to be a perfectly reflecting boundary. In the model used in Sect. 2, the downward (left) propagating waves will reflect of the density gradient in the chromosphere but this reflection happens gradually rather than at a defined point (such as the lower boundary in the simple model described above). Figure A.3 shows the relative mass increase for a 1D model with the same field-aligned thermodynamic equilibrium as our 2D model described in Sect. 2. As in the 2D model, the driver is now located at y = 7.8 Mm (the top of the first chromosphere). The region where the mass is evaluated is taken as y = 44−50 Mm (in the corona). As the reflection of the downward propagating waves is now happening gradually rather than instantaneously, the mass increase (solid line) is happening more gradually as well.

Fig. A.3. Relative mass increase ∫(ρ − ρ0)dV/∫ρ0dV with time in the region y = 44−50 Mm for the simulation that includes an atmosphere. The dashed line corresponds to the non-driven simulation.

All Figures

	Fig. 1. Contour plots of the temperature (left) and density (right) after the numerical relaxation.
In the text

	Fig. 2. Field-aligned temperature (left) and density (right) at x = −2 Mm (black), x = −1 Mm (blue) and x = 0 Mm (green). The small panel inserted on the temperature plot shows a zoomed version in the TR (between y = 7.5 Mm and y = 12 Mm). The symbols represent the numerical gridpoints.
In the text

	Fig. 3. Contour plots of the magnetic field components Bx (left) and By (middle) and the Alfvén speed (right) after the numerical relaxation.
In the text

	Fig. 4. Contour of vz (km s−1) for the viscous simulation at t = 93 s (left panel) and t = 983 s (right panel).
In the text

	Fig. 5. Mean field-aligned velocity vy integrated over the first TR-coronal boundary in the central part of the loop (−1 <  x <  1 Mm) for the viscous (solid line), ideal (dashed line) and non-driven (dotted line) simulations as a function of time.
In the text

	Fig. 6. Contours of (T − T0)/T0 (left) and (ρ − ρ0)/ρ0 (right) at t = 5340 s for the viscous simulation.
In the text

	Fig. 7. Time distance graphs of the relative changes in temperature (top) and density (bottom) along x = −1.18 Mm (marked by vertical dashed line in Fig. 6).
In the text

Fig. 8. Left: mass increase (kg m−1) for the three simulations (viscous, solid lines; ideal, dashed lines; non-driven, dotted lines). The black lines correspond to the integrated mass in the coronal part of the domain whereas the green lines represent the mass in the non-coronal part. Right: mass increase (kg m−1; blue lines) but now only integrated over the coronal part of the shell region (viscous, solid lines; ideal, dotted-dashed lines). The green lines represent the time integrated, averaged mass flux (where the averaging is done over the 4 boundaries of the shells).

In the text

	Fig. 9. Time integrated Poynting flux ${\displaystyle {\int }_0^t{\displaystyle \int {\left(E\times B\right)}_y}}\text{d}S\text{d}t(\text{J}{\text{m}}^{\text{-1}}$, with S the boundaries of the left shell, for the viscous (solid lines) and ideal simulation (dashed lines). The green lines represent the lower TR-coronal boundary and the black lines the upper boundary.
In the text

	Fig. A.1. Left: vz/max(vz) and Bz/max(Bz) with time at the location of the driver (y = 200 Mm). Right: upflows vy (km s−1) as a function of time at the location of the driver.
In the text

	Fig. A.2. Relative mass increase ∫(ρ − ρ0)dV/∫ρ0dV with time in the region y = 230−260 Mm. The red vertical lines mark the times when the slow wave, that is generated at t = 347 s at the location of the driver, arrives and leaves this region.
In the text

	Fig. A.3. Relative mass increase ∫(ρ − ρ0)dV/∫ρ0dV with time in the region y = 44−50 Mm for the simulation that includes an atmosphere. The dashed line corresponds to the non-driven simulation.
In the text