HEATING OF THE PARTIALLY IONIZED SOLAR CHROMOSPHERE BY WAVES IN MAGNETIC STRUCTURES

Understanding the physical mechanisms leading to heating of the solar atmosphere is one of the primary questions in solar physics and a long-standing puzzle. In recent years, the importance of partial ionization for the processes in solar plasma is becoming clear. It has been shown that the current dissipation enhanced in the presence of neutrals in the plasma not entirely coupled by collisions can be essential for the energy balance of the chromosphere (De Pontieu & Haerendel 1998; Judge 2008; Krasnoselskikh et al. 2010; Khomenko & Collados 2012a; Martínez-Sykora et al. 2012). In the latter paper it was shown how the dissipation of static currents in non-current-free solar magnetic flux tubes allows to release a large amount of energy into the chromosphere, potentially able to compensate its radiative losses, (Khomenko & Collados 2012b).

Solar atmosphere is far from stationary. It is filled with waves and rapidly varying flows, and pierced by dynamic magnetic structures, originating in the solar interior. As it has been demonstrated (Osterbrock 1961; Goodman 2005; Steiner et al. 2008; Goodman & Kazeminezhad 2010; Krasnoselskikh et al. 2010; Goodman 2011; Abbett & Fisher 2012; Shelyag et al. 2012, 2013), the magnetized plasma motions in the solar photosphere create an electromagnetic (Poynting) energy flux sufficient to heat the upper solar atmosphere. A significant part of this flux is produced in the form of Alfvén waves (Shelyag & Przybylski 2014). These waves are notoriously difficult to dissipate if an ideal plasma assumption is employed.

In the present study we extend the analysis by Khomenko & Collados (2012a) including time-dependent oscillatory perturbations. Magnetic waves (fast magneto-acoustic and Alfvén waves) produce perturbations in the magnetic field. These perturbations create currents at all times and locations in the magnetized solar atmosphere. Similarly to static currents, the currents produced by waves in a partially ionized magnetized chromosphere can be efficiently dissipated due to the presence of ambipolar diffusion. This would provide a constant and efficient influx of energy into the chromospheric layers.

It was long known that the presence of neutral atoms in partially ionized plasmas significantly affects excitation and propagation of waves (Kumar & Roberts 2003; Khodachenko et al. 2004; Forteza et al. 2007; Pandey et al. 2008; Vranjes et al. 2008; Soler et al. 2009; Zaqarashvili et al. 2011). The relative motion between the ionized and charged species causes an increase of collisional damping of MHD waves in the photosphere, chromosphere and prominence plasmas, in a similar way that any electro-magnetic wave is dissipated in a collisional medium. These mechanisms were investigated for high-frequency waves and were shown to be important for frequencies close to the collisional frequency between the different species (Khodachenko et al. 2004; Forteza et al. 2007; Pandey et al. 2008). In the cited works the damping was analyzed in terms of the damping time of waves. Here we study the energy input due to this damping.

In the solar chromosphere the collisional coupling of the plasma is still strong (Zaqarashvili et al. 2011; Khomenko et al. 2014). In this case, it is convenient to use a single-fluid quasi-MHD approach i.e., including non-ideal terms not taken into account by the ideal MHD approximation instead of solving full multi-fluid equations. The conservation equations for the multi-species solar plasma with scalar pressure and no heat flux can be written as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{u}})}{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{u}}\otimes {\boldsymbol{u}})+{\rm{\nabla }}p={\boldsymbol{J}}\times {\boldsymbol{B}}+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\boldsymbol{\nabla }}\cdot ({\boldsymbol{u}}e)+p{\rm{\nabla }}\cdot {\boldsymbol{u}}={\boldsymbol{J}}\cdot [{\boldsymbol{E}}+{\boldsymbol{u}}\times {\boldsymbol{B}}],\end{eqnarray} \tag{ 3 }$$

where the internal energy e is computed from the scalar pressure according to the ideal equation of state:

$$\begin{eqnarray}&&e=p/(\gamma -1).\end{eqnarray} \tag{ 4 }$$

The rest of the notation is standard (for details see Khomenko et al. 2014). This system of equations is closed by the generalized Ohm's law, where we only include Ohmic and ambipolar terms:

$$\begin{eqnarray}&&[{\boldsymbol{E}}+{\boldsymbol{u}}\times {\boldsymbol{B}}]=\eta {\boldsymbol{J}}-{\eta }_{A}\displaystyle \frac{[({\boldsymbol{J}}\times {\boldsymbol{B}})\times {\boldsymbol{B}}]}{| B{| }^{2}}\end{eqnarray} \tag{ 5 }$$

with resistivity coefficients (in units of ${{ML}}^{3}/{{TQ}}^{2}$) defined as

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\rho }_{e}}{{({{en}}_{e})}^{2}}\left(\displaystyle \sum _{a}{\nu }_{e;{ai}}+\displaystyle \sum _{b}{\nu }_{e;{bn}}\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\eta }_{A}=\displaystyle \frac{{\xi }_{n}^{2}| B{| }^{2}}{{\alpha }_{n}},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{\alpha }_{n}=\displaystyle \sum _{b}{\rho }_{e}{\nu }_{e;{bn}}+\displaystyle \sum _{a}\displaystyle \sum _{b}{\rho }_{a}{\nu }_{{ai};{bn}}.\end{eqnarray} \tag{ 8 }$$

In these equations the inter-species collision frequencies ν are those between electrons and ions of the specie a (${\nu }_{e;{ai}}$), between electrons and neutrals of the specie b (${\nu }_{e;{bn}}$) and those between ions and neutrals of different species (${\nu }_{{ai};{bn}}$). These frequencies, together with the densities ρ, the electron number density n_e, and the neutral fraction ${\xi }_{n}$ are calculated according to Khomenko & Collados (2012a), Khomenko et al. (2014). We do not include the rest of the terms in the generalized Ohm's law as these terms either do not contribute to the energy balance, as Hall term, or can be considered small for the solar atmospheric conditions (Khomenko et al. 2014). Nevertheless, it should be noted that the battery term can potentially increase the current and, therefore, the dissipation, as it was demonstrated by Díaz et al. (2014).

Finally, the induction equation reads as:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times \left[({\boldsymbol{u}}\times {\boldsymbol{B}})-\eta {\boldsymbol{J}}+{\eta }_{A}\displaystyle \frac{[({\boldsymbol{J}}\times {\boldsymbol{B}})\times {\boldsymbol{B}}]}{| B{| }^{2}}\right].\end{eqnarray} \tag{ 9 }$$

The physical magnitude of η is too small to produce a measurable effect. We used a numerical analog of the η, denoted as ${\eta }_{{\rm{hyp}}}$ (hyperdiffusion). The value of ${\eta }_{{\rm{hyp}}}$ is allowed to vary in space and time according to the structures developed in the simulations, so that it reaches larger values at places where the variations are small-scale (Vögler et al. 2005; Felipe et al. 2010). ${\eta }_{{\rm{hyp}}}$, together with the grid resolution, defines the maximum Reynolds number of the simulations. The effective magnetic Reynolds number due to ${\eta }_{{\rm{hyp}}}$ in the simulation presented in this paper (as computed directly from the derivatives of physical quantities entering its definition) varies between 10⁰ and 10³, with the smallest values localized in narrow layers at the regions where the gradients are strongest, i.e., at the shocks. ${\eta }_{A}$, however, can reach values significantly larger than ${\eta }_{{\rm{hyp}}}$ decreasing the effective physical magnetic Reynolds number to 10⁻¹–10⁰ in the layers above 1200 km, where ${\eta }_{A}$ is large.

The equations above are solved by means of the code mancha (Khomenko et al. 2008; Felipe et al. 2010; Khomenko & Collados 2012a). We evolve n_e in time via the Saha equation. The simulations are done in full 3D. We use 10 grid points Perfectly Matched Layer (PML, Berenger 1994) boundary conditions on the side and top boundaries. At the bottom boundary a time-dependent condition for all variables is provided to drive waves into the domain.

A single 3D Schlüter–Temesvary-like axisymmetric self-similar magnetic flux tube, embedded into the VALIIIC solar atmospheric model (Vernazza et al. 1981), similar to the one in Fedun et al. (2011a, 2011b), Gent et al. (2013), is used as the initial equilibrium model in the simulation (see Figure 1). The components of the magnetic field are calculated as follows:

$$\begin{eqnarray}&&{B}_{x}=-{{xB}}_{0}(z)G\displaystyle \frac{{{dB}}_{0}(z)}{{dz}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{B}_{y}=-{{yB}}_{0}(z)G\displaystyle \frac{{{dB}}_{0}(z)}{{dz}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{B}_{z}={B}_{0}^{2}(z)G,\end{eqnarray} \tag{ 12 }$$

where G is a piecewise-parabolic function with the dimension of ${B}^{-1}$, describing the opening of the magnetic flux tube. ${B}_{0}(z)$ is chosen so that it decreases with height, similar to the gas pressure. This choice allows to avoid negative values of gas pressure within the magnetic flux tube, while keeping plasma-β low. The magnetic field components were used to recalculate (using the magnetohydrostatic equation) the pressure and density to obtain a 3D magnetohydrostatic configuration. The simulation box is set to have physical dimensions of 4 × 4 × 1.84 ${\mathrm{Mm}}^{3}$, with a uniform resolution of 10 km. The bottom boundary of the simulation box is located in the deep photosphere. The parameters of the magnetic field are chosen such that the magnetic field strength at the axis of the magnetic flux tube is $1.4\;{\rm{kG}}$ at the bottom (z = 0) and $10\;{\rm{G}}$ at the top ($z=1.84\;{\rm{Mm}}$) of the domain. This leads to the plasma $\beta =1$ height of about $0.7\;{\rm{Mm}}$ at the axis of the tube.

Zoom In Zoom Out Reset image size

Self-consistent periodic perturbations of the velocity vector, pressure, and density are introduced according to Mihalas & Mihalas (1984). For the simulations described here we used the driving frequency of $\nu =\omega /2\pi =25$ mHz. The choice of the driving frequency was motivated by the results of Shelyag et al. (2013) who showed that Alfvén waves with such periods can be generated in intergranular magnetic flux tubes by turbulent convective motions. To break the symmetry, the acoustic pulse was placed at $100\;{\rm{km}}$ off the flux tube axis in the field-free atmosphere and was limited in horizontal extent by a Gaussian of $100\;{\rm{km}}$ FWHM. The perturbation amplitude was set to 500 m s⁻¹. The duration of the simulation was $350\;{\rm{s}}$.

In order to understand the effect of plasma oscillations on chromospheric heating, we perform three simulation runs with exactly the same numerical parameters, but different driving. In the first run (denoted as AD below), we include the ambipolar term as a perturbation. This case is analogous to Khomenko & Collados (2012a). In the second run (denoted as W), we include only the wave driving, with no ambipolar term. In the third run (denoted as ADW), we include both the wave driving and the ambipolar diffusion term. In the AD run, the ambipolar term alone acting on a static flux tube produces dissipation of its static currents, which creates perturbations in all thermodynamical parameters and, in particular, in temperature. In the ADW case, both wave and ambipolar drivings produce variations of temperature and therefore it is difficult to separate the effects of waves from the effect of static current dissipation. For that we perform a flatfield-like procedure. In order to determine the heating caused by the current dissipation of the currents produced only by the waves, we subtract from ADW case the heating produced in the AD case, as well as the variations due to waves produced in the W case.

The acoustic pulse set in the field-free atmosphere generates essentially acoustic (fast) waves. These waves gradually enter the flux tube, cross the c_s = v_a level and generate fast and slow magneto-acoustic waves and Alfvén waves via the mode transformation. The latter waves propagate in the magnetic flux tube and create perpendicular currents, which then dissipate due to ion–neutral friction (see Figure 1). This process is demonstrated in Figure 2, where the logarithm of the ambipolar diffusion heating term $\mathrm{log}{\eta }_{A}{J}^{2}{\mu }_{0}/\rho$ is shown. In the first column, the heating term is shown for $100\;{\rm{s}}$ of simulation time. At this time, the perturbation has not reached the upper layers of the domain above $1\;{\rm{Mm}}$. The heating term values are below ${10}^{4}\;{\rm{J}}\;{\mathrm{kg}}^{-1}\;{{\rm{s}}}^{-1}$, with the location of the largest values at the top of the domain. The values of ${\eta }_{A}$ reach their maximum in the interior of the flux tube, while the static currents are largest at the flux tube boundaries. As a consequence, the heating by dissipation of static currents is strongest at the tube walls. Note that the wall heating is present at all heights, but is largest in the upper part of the domain, as a consequence of the strong increase of ${\eta }_{A}$ with height (see similar case in Khomenko & Collados 2012a). At $150\;{\rm{s}}$ the situation changes. The wave perturbations reach the upper part of the domain, and the maximum of the heating term is localized in the regions around the chromospheric shocks. Further on, at $200\;{\rm{s}}$ (third column of Figure 2), the ambipolar diffusion heating rate reaches ${10}^{6}\;{\rm{J}}\;{\mathrm{kg}}^{-1}\;{{\rm{s}}}^{-1}$, which is two orders of magnitude higher than for the case of dissipation of static currents (first column of the figure). The high heating rate region is now not confined to the chromospheric shock anymore, and occupies the whole upper part of the domain above $\sim 1.2\;{\rm{Mm}}$ height.

Zoom In Zoom Out Reset image size

In order to separate the different modes in the atmosphere where ${v}_{a}\gt {c}_{s}$, and study their behavior and relation to heating, we projected the velocities in the computational domain onto the three characteristic directions, defined by the following vectors (Khomenko & Cally 2012; Cally & Khomenko 2015):

$$\begin{eqnarray}&&{\hat{e}}_{{\rm{long}}}=[\mathrm{cos}\phi \mathrm{sin}\theta ,\mathrm{sin}\phi \mathrm{sin}\theta ,\mathrm{cos}\theta ];\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\hat{e}}_{| | }=[\mathrm{sin}\phi ,-\mathrm{cos}\phi ,0];\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\hat{e}}_{\perp }=[\mathrm{cos}\phi \mathrm{cos}\theta ,\mathrm{sin}\phi \mathrm{cos}\theta ,-\mathrm{sin}\theta ].\end{eqnarray} \tag{ 15 }$$

Here, θ is the inclination and ϕ is the azimuth of the magnetic field. The first vector ${\hat{e}}_{{\rm{long}}}$ is directed along the magnetic field lines. The second vector ${\hat{e}}_{| | }$ gives the direction perpendicular to ${\boldsymbol{B}}$ in the plane containing the field lines and gravity. The third vector is perpendicular to the first two. These projections allow to separate the fast wave in the ${\hat{e}}_{\perp }$ projection and the Alfvén wave in the ${\hat{e}}_{| | }$ projection. This was used to analyze the relative influence of slow and Alfvén waves on the perpendicular current, and, subsequently, on the ambipolar heating in the simulation. Figure 3 shows the occurrence plots of the perpendicular current ${J}_{\perp }$ versus the magnetic field fluctuations along ${\hat{e}}_{| | }$ (Alfvén wave, left panel), and along ${\hat{e}}_{\perp }$ (fast wave, right panel). Clearly, Alfvén waves have a stronger effect on the chromospheric non-ideal plasma heating. Furthermore, the perpendicular current increases with the increase of Alfvén wave perturbation in the magnetic field, while the occurrence plot shows no significant dependence of the perpendicular current on the fast wave magnetic field perturbation amplitude.

Zoom In Zoom Out Reset image size

Slow waves, aligned with ${\hat{e}}_{{\rm{long}}}$, are also produced by the oscillatory source in our simulations. However, they do not perturb the magnetic field in the perpendicular to the field line direction, therefore have only a minor effect on the energy balance in the modeled chromosphere.

Through the heating term in the energy equation Equation (3) associated with ${\eta }_{A}$, the wave energy is converted into thermal energy. This process is demonstrated in Figure 4. In this plot, the time dependences of the volume-averaged thermal (black curve), kinetic (blue curve) and magnetic (red curve) energy densities per unit volume are shown. The averaging of the energy density values is carried over the volume around the tube axis, extending in the simulated chromosphere between 1.3 and $1.65\;{\rm{Mm}}$. The volume boundaries lie within the plasma $\beta \lt 1$ surface, which also coincides with the region where ${\eta }_{A}$ is maximal. In the left panel of this figure, it is straightforward to see that the thermal energy density in the AD simulation increases over time due to the dissipation of static currents in the initial non-current-free magnetic field configuration. Correspondingly, the magnetic field energy density decreases over time.

Zoom In Zoom Out Reset image size

The right panel of Figure 4 shows the difference between the corresponding energy densities in the ADW simulation and W simulation. The energy densities from the AD run (left column of Figure 4) were also subtracted. It is evident from the plot that the dissipation of the currents induced by the oscillatory source leads to the decrease of the kinetic and magnetic energy density, as the magnetic energy gain due to the source action is the same in the ambipolar and ideal simulations, while the additional dissipation present in the ambipolar simulation decreases the magnetic field energy and converts it into the thermal energy. The average thermal energy per unit volume being deposited to the low-β chromospheric plasma due to wave dissipation is about 10–20 times greater than the one deposited due to the dissipation of static currents.

The computed spatially averaged thermal energy density difference due to dissipation of the perpendicular currents caused by waves corresponds to the additional temporally averaged thermal energy flux of $\sim 100\;{\rm{J}}\;{{\rm{m}}}^{-2}\;{{\rm{s}}}^{-1}$ through the magnetized chromosphere, assuming uniform heating distribution for simplicity (similar to van Ballegooijen et al. 2011). While this estimate is about one order of magnitude less than generally accepted value required to compensate the radiative energy losses ($2\div4\cdot {10}^{3}\;{\rm{J}}\;{{\rm{m}}}^{-2}\;{{\rm{s}}}^{-1}$, Osterbrock 1961; Shelyag et al. 2012; Arber et al. 2016), the Poynting flux in our simulations at the chromospheric level (z = 0–0.5 Mm) is equal to $\sim 700\;{\rm{J}}\;{{\rm{m}}}^{-2}\;{{\rm{s}}}^{-1}$, which is also nearly an order of magnitude less than that value. Therefore, our result suggests that waves, produced by photospheric motions, are an important energy source to provide heating of the solar chromosphere.

Further evidence is provided by the measurements of Poynting flux absorption. The left panel of Figure 5 shows horizontally averaged Poynting flux in the simulations with the oscillatory source and with (black curves) and without (red curves) ambipolar diffusion. The averaging is carried out over a narrow region within $200\;\mathrm{km}$ range around the tube axis. The Poynting flux in the W simulation decreases with height due to the tube expansion out of the averaging region boundaries, and does not behave differently from the ADW simulation in the lower solar atmosphere. Starting from $1.1\;{\rm{Mm}}$ height, the ideal and non-ideal simulation curves diverge. The decrease in the Poynting flux by the factor of up to 0.8 is produced by ambipolar diffusion only. The Poynting flux absorption dependence on height is shown in the right panel of Figure 5. The spatially averaged absorption coefficient is $d(1-{F}_{{ADW}}/{F}_{W})/{dz}=0.47\quad {\mathrm{Mm}}^{-1}$ for the tube axis region.

Zoom In Zoom Out Reset image size

In this paper we investigated chromospheric heating by oscillations in the chromospheric magnetic fields. An initial study was performed using a single magnetic flux tube, rooted in the photosphere in magneto-hydrostatic equilibrium. A monochromatic oscillatory perturbation was placed off the tube axis in the deep photosphere to generate magneto-acoustic oscillations in the magnetic flux tube. The currents, corresponding to fast magneto-acoustic and Alfvén waves, have been shown to be dissipated by the ambipolar diffusion mechanism in the simulated chromosphere. It has been shown that for a source amplitude of $500\;{\rm{m}}\;{{\rm{s}}}^{-1}$, which mimicks photospheric motions, the heating due to the perturbation is two orders of magnitude larger than that due to the static currents dissipation. The perturbation we used generates all types of MHD waves in the magnetic flux tube. It has been shown, however, that the Alfvénic perturbation component has a stronger effect on the heating.

Generally, it is expected that oscillations with higher (1–10 Hz) frequencies are more efficiently dissipated into heat compared to lower (similar to the ones used in this paper) frequency oscillations (e.g., Arber et al. 2016; Soler et al. 2015). On the other hand, a power law is expected for the velocity power spectrum in the solar photosphere, including its high-frequency part, currently not observed (Goldreich et al. 1994; Musielak et al. 1994; Stein & Nordlund 2001; Fossum & Carlsson 2005), leading to less power at higher frequencies. Consequently, it is currently difficult to draw a unique conclusion on a sole mechanism for the chromospheric heating based on a single simulation. Nevertheless, this study should be understood as an initial proof of concept of the chromospheric energy balance based on MHD wave absorption in non-uniform three-dimensional magnetized chromospheric plasmas, with a more detailed study of the absorption on the source frequency, amplitude and simulation resolution to follow.

Furthermore, as the simulation has to resolve the smallest spatial and temporal scales in the system, which in the case of frequencies of the order of $10\;{\rm{Hz}}$ is defined by the oscillation wavelength and period, an extremely high spatial resolution of the order of few hundred meters would be necessary to directly test high-frequency currents dissipation. Therefore, even with currently available computational resources some extrapolation would be required.

E.K. was partially supported by the Spanish Ministry of Science through projects AYA2010-18029, AYA2011-24808, and AYA2014-55078-P. This work contributes to the deliverables identified in FP7 European Research Council grant agreement 277829, "Magnetic connectivity through the Solar Partially Ionized Atmosphere." This research was supported by the SOLARNET project, funded by the European Commission's FP7 Capacities Programme under the Grant Agreement 312495. This research was undertaken with the assistance of BSC-CNS MareNostrum facilities, the resources provided at the Centre for Astrophysics & Supercomputing of Swinburne University of Technology (Australia), NCI National Facility systems at the Australian National University through the National Computational Merit Allocation Scheme supported by the Australian Government, and at the Multi-modal Australian ScienceS Imaging and Visualization Environment (MASSIVE).