Time-Dependent Two-Fluid Magnetohydrodynamic Model and Simulation of the Chromosphere

TL;DR

This paper introduces a two-fluid MHD model for the solar chromosphere that captures ion-neutral interactions, demonstrating that jets can be thermally driven by frictional heating rather than magnetic reconnection.

Contribution

The study develops a novel two-fluid MHD simulation including ion-neutral interactions, revealing thermal jet driving mechanisms in the chromosphere.

Findings

Magnetic canopy forms quickly and is disrupted by jets.

Jets are driven mainly by frictional heating, not magnetic reconnection.

The model self-consistently produces shocks in the chromosphere.

Abstract

The sun's chromosphere is a highly dynamic, partially-ionized region where spicules (hot jets of plasma) form. Here we present a two-fluid MHD model to study the chromosphere, which includes ion-neutral interaction and frictional heating. Our simulation recovers a magnetic canopy shape that forms quickly, but is also quickly disrupted by the formation of a jet. Our simulation produces a shock self-consistently, where the jet is driven by the frictional heating, which is much greater than the ohmic heating. Thus, our simulation demonstrates that the jet could be driven purely by thermal effects due to ion-neutral collisions and not by magnetic reconnection. We plan to improve the model to include photo-chemical effects and radiation.

Equations

$$\frac{\partial {\rho }_{\text{i}}}{\partial t}+\nabla \cdot {\rho }_{\text{i}}\text{mathbfit}V=0$$

$$\frac{\partial {\rho }_{\text{n}}}{\partial t}+\nabla \cdot {\rho }_{\text{n}}\text{mathbfit}U=0$$

$$\text{ilabel}eq:plasmamomentum\frac{\partial {\rho }_{\text{i}}\text{mathbfit}V}{\partial t}=-\nabla \cdot [{\rho }_{\text{i}}\text{mathbfit}V\text{mathbfit}V+n_{\text{i}}{\text{k}}_{\text{b}}{T}_{\text{i}}+\frac{B^2}{2{\mu }_0}\text{mathbfit}I-\frac{\text{mathbfit}BB}{{\mu }_0}]+{\rho }_{\text{i}}\text{mathbfit}g-{\rho }_{\text{i}}{\nu }_{\text{in}}(\text{mathbfit}V-\text{mathbfit}U)$$

$$\text{ilabel}eq:neutralmomentum\frac{\partial {\rho }_{\text{n}}\text{mathbfit}U}{\partial t}=-\nabla \cdot [{\rho }_{\text{n}}\text{mathbfit}U\text{mathbfit}U+n_{\text{n}}{\text{k}}_{\text{b}}{T}_{\text{n}}]+{\rho }_{\text{n}}\text{mathbfit}g+{\rho }_{\text{i}}{\nu }_{\text{in}}(\text{mathbfit}V-\text{mathbfit}U)$$

$$\frac{\partial \text{mathbfit}B}{\partial t}=-\nabla \cdot (\text{mathbfit}V\text{mathbfit}B-\text{mathbfit}B\text{mathbfit}V)-\nabla \times (\eta \text{mathbfit}J+\frac{1}{n_{\text{e}}e}\text{mathbfit}J\times \text{mathbfit}B)$$

$$\text{ilabel}eq:energyion\begin{array}{r}\frac{\partial {\text{k}}_{\text{b}}{T}_{\text{i}}}{\partial t}=-\nabla \cdot ({\text{k}}_{\text{b}}{T}_{\text{i}}\text{mathbfit}V)+\frac{1}{3}{\text{k}}_{\text{b}}{T}_{\text{i}}\nabla \cdot \text{mathbfit}V-\frac{2}{3}\frac{\nabla \cdot ({\text{mathbfit}q}_{\text{i}}+\eta \text{mathbfit}B\times \text{mathbfit}J)}{n_{\text{i}}}\\ +\frac{2}{3}\frac{{\text{m}}_{\text{i}}}{{\text{m}}_{\text{i}}+{\text{m}}_{\text{n}}}{\nu }_{\text{in}}[3{\text{k}}_{\text{b}}(T_{\text{n}}-T_{\text{i}})+{\text{m}}_{\text{n}}|\text{mathbfit}V-\text{mathbfit}U{|}^2]\end{array}$$

$$\text{ilabel}eq:energyneutral\begin{array}{r}\frac{\partial {\text{k}}_{\text{b}}{T}_{\text{n}}}{\partial t}=-\nabla \cdot ({\text{k}}_{\text{b}}{T}_{\text{n}}\text{mathbfit}U)+\frac{1}{3}{\text{k}}_{\text{b}}{T}_{\text{n}}\nabla \cdot \text{mathbfit}U-\frac{2}{3}\frac{\nabla \cdot {\text{mathbfit}q}_{\text{n}}}{n_{\text{n}}}\\ +\frac{2}{3}\frac{{\text{m}}_{\text{n}}}{{\text{m}}_{\text{i}}+{\text{m}}_{\text{n}}}{\nu }_{\text{in}}[3{\text{k}}_{\text{b}}(T_{\text{i}}-T_{\text{n}})+{\text{m}}_{\text{i}}|\text{mathbfit}V-\text{mathbfit}U{|}^2]\end{array}$$

$$\nabla \times \text{mathbfit}B={\mu }_0\text{mathbfit}J$$

$$\text{mathbfit}q=-\kappa \nabla T$$

$${\nu }_{\text{in}}$$

$$\text{ilabel}eq:bglm\frac{\partial \text{mathbfit}B}{\partial t}=-\nabla \cdot (\text{mathbfit}V\text{mathbfit}B-\text{mathbfit}B\text{mathbfit}V)-\nabla \times (\eta \text{mathbfit}J+\frac{1}{n_{\text{e}}\text{e}}\text{mathbfit}J\times \text{mathbfit}B)+\nabla \Phi$$

$$\text{ilabel}eq:glm\frac{\partial \Phi }{\partial t}+c_{\text{h}}^2\nabla \cdot \text{mathbfit}B=-\frac{c_{\text{h}}^2}{c_{\text{p}}^2}\Phi$$

$$\frac{\partial \nabla \cdot \text{mathbfit}B}{\partial t}={\nabla }^2\Phi$$

$$f_{j-1/2}^{\ast }=\frac{1}{2}([f(u_{j-1/2}^{\text{r}})+f(u_{j-1/2}^{\text{l}})]-{\lambda }_{j-1/2}[u_{j-1/2}^{\text{r}}-u_{j-1/2}^{\text{l}}])$$

$$f_{j+1/2}^{\ast }=\frac{1}{2}([f(u_{j+1/2}^{\text{r}})+f(u_{j+1/2}^{\text{l}})]-{\lambda }_{j+1/2}[u_{j+1/2}^{\text{r}}-u_{j+1/2}^{\text{l}}])$$

$$u_{j-1/2}^{\text{l}}=u_{j-1}+\frac{1}{2}\varphi (r_{j-1})(u_j-u_{j-1})$$

$$u_{j+1/2}^{\text{l}}=u_j+\frac{1}{2}\varphi (r_j)(u_j-u_{j-1})$$

$$u_{j-1/2}^{\text{r}}=u_j+\frac{1}{2}\varphi (r_j)(u_{j+1}-u_j)$$

$$u_{j+1/2}^{\text{r}}=u_{j+1}+\frac{1}{2}\varphi (r_{j+1})(u_{j+1}-u_j)$$

$$r_j=\frac{u_j-u_{j-1}}{u_{j+1}-u_j}$$

$$\text{ilabel}minmod\varphi (r)=\max (0,\min (1,r)),\underset{r\to \infty }{\lim }\varphi (r)=1$$

$$u^{n+1}+R_{\text{I}}(u^{n+1})\Delta t=u^n+R_{\text{E}})\Delta t$$

$$c_{\text{fast}}=\frac{1}{\sqrt{2}}[c_{\text{sound}}^2+c_{\text{Alfveˊn}}^2+\sqrt{(c_{\text{sound}}^2+c_{\text{Alfveˊn}}^2{)}^2-(2c_{\text{sound}}{c}_{\text{Alfveˊn}}{)}^2}{]}^{1/2}$$

$$\Delta t=0.1\frac{\Delta z}{\max (V+c_{\text{fast}})}$$

$$V_{\perp }(t,z=0)=V_{\text{b}}[\sum _{k=0}^s{({\omega }_k/{\omega }_{\text{p}})}^{5/6}\cos ({\omega }_{k}t+{\Phi }_k)$$

$$+\sum _{j=s+1}^m{({\omega }_{\text{p}}/{\omega }_j)}^{5/6}\cos ({\omega }_{j}t+{\Phi }_j)],$$

Full text

\setlastpage\inarticletrue{opening}

1 Introduction

The base of the outer solar atmosphere (right above the photosphere), referred to as the chromosphere, is highly dynamic and is strongly influenced by magnetic-field activity (Hasan, 2008). It is partially ionized, with as much as a thousand times more neutral hydrogen than ions (Alissandrakis et al., 2018). Song and Vasyliūnas (2011) and Leake et al. (2014) discussed in detail the effects of ion–neutral coupling in the chromosphere by considering neutrals and collisions into the momentum and Maxwell equations.

A local-scale configuration commonly used in chromosphere models are so-called supergranules, convection cells of 30 Mm in the horizontal directions. Supergranules are separated by strong magnetic-field regions, referred to as a “network”. It is believed that the network regions correspond to strong magnetic-field at the photosphere. Because the thermal-pressure decreases with height while the total magnetic flux is conserved, the magnetic-field tends to expand in higher altitudes to form a canopy or “wine glass” shaped magnetic geometry so that total pressure, the magnetic pressure, $p_{\text{B}}$, plus thermal-pressure, $p_{\text{T}}$, i.e. $p_{\text{T}}+p_{\text{B}}$, is balanced at each height (Gabriel, 1976). This geometry is reflected in a significant altitude change in the plasma beta, which is the ratio of the thermal to magnetic pressures, $\beta=p_{\text{T}}/p_{\text{B}}$. The chromosphere has much lower temperature than that of the hot corona. The average temperature of the chromosphere is around 6000 K while the coronal temperature is over million degrees (Avrett and Loeser, 2008). This dramatic change in the temperature, taking place in a very narrow region of a few hundred kilometers referred to as the transition region, is known as the coronal heating problem and it is a long-standing problem in solar physics.

In the chromosphere, shocks are often observed and jets are commonplace in the form of spicules. Spicules are local, hot, upward jets of plasma observed in the corona. The origination and mechanism behind spicule formation are not well known. Spicules are categorized into two types (De Pontieu et al., 2007). Type I spicules have longer lifetimes and are typically larger (five minutes and $\approx$15 Mm). Type II spicules have shorter lifetimes and are smaller (less than two minutes and $\approx$6 Mm). De Pontieu et al. (2004) suggested that photospheric oscillations give birth to spicules. It was also suggested that magnetic reconnection may be the cause of Type II spicules (De Pontieu et al., 2007).

There is strong observational evidence for the ubiquity of Alfvén waves in the Sun’s atmosphere (Aschwanden et al., 1999; Nakariakov et al., 1999; Ofman and Wang, 2008; Jess et al., 2009; McIntosh et al., 2011; Mathioudakis, Jess, and Erdélyi, 2013). These observations of the solar limbs show Alfvén waves with amplitudes in the range of 5 – 25 km s*-1* and phase speeds 250 – 500 km s*-1*. De Pontieu et al. (2007) have shown that Alfvén waves are carried by spicules with significant amplitudes of order 20 km s*-1*. It has been suggested that Alfvén waves could be damped in the chromosphere due to collisions between plasma ions and neutrals and cause sufficient heating (Song and Vasyliūnas, 2011) and (Song and Vasyliūnas, 2014). Song (2017) further proposed that the nonuniform heating associated with the Alfvén waves may explain the formation of the spicules. Brady and Arber (2016) studied heating due to Alfvén and kink waves turning into shocks and dissipating while still providing sufficient heating albeit a different mechanism. Khomenko (2017) looks at other waves as well as Alfvén waves and their effects in weakly ionized solar plasma.

Recent simulations and models of the chromosphere have adopted a single-fluid approach (Gudiksen et al., 2011; Brady and Arber, 2016; Ni et al., 2016; Kuźma et al., 2017). In order to account for neutrals’ effect on the plasma and electric field an “ambipolar diffusion” term is added to the Generalized Ohm’s Law to account for collisions (Leake and Arber, 2006; Arber, Haynes, and Leake, 2007; Martinez-Sykora, De Pontieu, and Hansteen, 2012; Martinez-Sykora et al., 2017). This is useful to study the highly magnetized, low-$\beta$, environment of the upper chromosphere. A two-fluid approach takes into account separate fluid motions of plasma and neutrals especially in highly collisional environments like the low solar atmosphere in regions of high magnetic gradients. Soler et al. (2013) studied Alfvén-wave propagation in partially ionized plasma while taking a generalized approach with free parameters with the lower solar atmosphere case as an example. Maneva et al. (2017) and Wójcik, Murawski, and Musielak (2018) studied acoustic-wave propagation using a collisional two-fluid approach as well. Three-fluid simulations, with electrons being the third fluid, were carried out to study reconnection events in the chromosphere (Leake, Lukin, and Linton, 2013). In these highly localized small regions, electron motion is separate from ion and neutral motion and its effects cannot be neglected. The Bifrost (Gudiksen et al., 2011) code employs a very large artificial viscosity for stability and uses the one-fluid approach. By self-consistently solving plasma, neutral fluids, and Maxwell’s equations, Tu and Song (2013) have shown that due to Joule dissipation Alfvén waves can provide strong enough chromospheric-heating and that the heating rate is consistent with the theoretical model of Song and Vasyliūnas (2011). Alternatively, it was also shown that the coupling of Alfvén waves with magneto-acoustic waves leads to shock heating that causes sufficient chromospheric-heating (Brady and Arber, 2016).

The aim for this article is to introduce a two-dimensional, self-consistent, and conservative magnetohydrodynamic (MHD) model for the chromosphere that includes heating. In particular, the model accounts for the neutrals as an additional fluid with its interaction to the ions as advocated by Song and Vasyliūnas (2011). The heating presented in this model is collisional in nature (Vasyliūnas and Song, 2005), and it is separated into two terms: an ohmic term due to electron collisions and frictional term due to ion–neutral collisions. This model also shows a spontaneous eruptive event due to heating in the magnetic flux concentrations of the solar magnetic network. We show that the spontaneous eruptive event can occur due to the frictional-heating, while magnetic reconnection is not a dominant factor or a requirement.

In Section \irefsec:model and Section \irefsec:domain the details of the model and code are explained. The validation of the code is carried out in Section \irefsec:validation. It is then followed by Section \irefsec:results where expected behavior from the model, such as formation of the magnetic canopy is shown and the eruption is thoroughly investigated.

2 An MHD Model for the chromosphere

\ilabel

sec:model

2.1 Governing Equations

In the model presented here, a two-fluid, collisional approach is used to account for the friction between ions and neutrals. The first fluid, the ion fluid, is similar to the traditional, single-fluid plasma, and assumes charge quasi-neutrality. With Maxwell’s equations, this fluid is used to calculate self-consistently the evolution of the plasma motion, magnetic fields, and currents in the system. The second, neutral-fluid, accounts for the neutral particles in the system. The two fluids interact by collisions. The friction between ions and neutrals affects the neutral fluid motion so that the neutrals experience the electromagnetic effects.

This is similar to the two-fluid HiFi code used for studying reconnection in the chromosphere (Leake et al., 2012) but without the electron pressure. Electrons are not treated as a separated fluid in this model because the electron’s mass is more than three orders of magnitude smaller than ion or neutral particle mass, so its contribution to heating and the equations, considering the scales, is negligible. The effects of its motion are included in the current terms. We use empirical number-density model that comes from observations assuming steady-state in the chromosphere (Avrett and Loeser, 2008). Here, we do not account for ionization and recombination which we plan to implement photo-chemistry in future work.

Our model uses Cartesian, 2.5 dimensionality, which means that all the quantities are functions of horizontal [$x$] and vertical [$z$] directions but uniform in the $y$-direction, while vector quantities have all three components including the third component that is perpendicular [$\perp$] to the $x$–$z$-plane. The perpendicular components are used to generate Alfvén waves at the photospheric boundary.

The governing equations are adapted from Tu and Song (2016) for the two-fluid (plasma and neutral) system as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here i, e, and n subscripts denote proton, electron and neutral, respectively. The primitive scalars m, $\rho$, $T$ are particle mass, mass density, and temperature, respectively. ${\mathbfit V}$, ${\mathbfit U}$, and ${\mathbfit B}$ are ion velocity, neutral velocity, and the magnetic-field, respectively, with ${\mathbfit g}$ denoting the solar gravitational acceleration near the solar surface. $\kappa$ is the Spitzer conductivity. $\nu_{ij}$ denotes collision rate of species $i$ to species $j$, and $\eta=\text{m}_{\text{e}}(\nu_{\text{en}}+\nu_{\text{ei}})/\text{e}^{2}n_{\text{e}}$ is the ohmic resistivity. We use the following collision rates (see De Pontieu, Marten, and Hudson, 2001, Appendix A)

[TABLE]

where $\ln\Lambda$ is the Coulomb logarithm taken as 10 for simplicity and $n$ is in cm*-3* and $T=\frac{T_{\text{i}}+T_{\text{n}}}{2}$ in K.

An in-depth derivation of weakly ionized plasma equations of different forms has been given by Leake et al. (2014). In the case we have chosen for the solar chromosphere and its two-fluid nature, the exchange of heat due to frictional-heating is proportional to the mass of the species, for this reason the electron species is ignored. A detailed derivation of frictional-heating has been given by Vasyliūnas and Song (2005). The thermal-energy equations are not total energy equations so this is the only non-conservative variable used; however, this was chosen to allow for the semi-implicit scheme to implicitly solve collisional terms in regions of the chromosphere where the collision frequency is much larger than the Alfv́en frequency.

In this work, we do not include an electron pressure in our momentum equation, as we expect that it will only matter in regions of small length scales and large current. We are currently working on adding the electron pressure to the code, and we plan future investigations to study its role and importance.

2.2 Divergence Cleaning

In order to maintain a divergence-less solution, the simulation employs a general Lagrange multiplier method (Dedner et al., 2002). This introduces a new unknown, the Lagrange multiplier $\Phi$, into the model:

[TABLE]

[TABLE]

The scalars $c_{\text{h}}$ and $c_{\text{p}}$ stand for the hyperbolic and parabolic speeds of the divergence errors.

It can be easily shown by taking the divergence of Equation \irefeq:bglm that the introduction of divergence in the magnetic-field is correlated with the multiplier.

[TABLE]

The divergence errors are assigned into hyperbolic and parabolic parts, where the errors in divergence advects and dissipates away from where they were created. The hyperbolic speed is chosen to be the fastest speed in the domain, the Alfvén speed $c_{\text{h}}=\max(V_{\text{A}})$, and the ratio $\alpha_{\text{p}}=\Delta x\frac{c_{\text{h}}}{c_{\text{p}}^{2}}$ (see Mignone, Tzeferacos, and Bodo, 2010, p. 5899, 5911), also known as the damping parameter, is kept at a constant of 0.5.

This method keeps the divergence errors bounded to the order of the grid size in the simulation.

2.3 Numerical Scheme

The model’s finite-volume method and second-order scheme is the Total Variation Diminishing Monotonic Upwind Scheme for Conservation Laws (TVD-MUSCL) to avoid spurious oscillations and capture shocks that may arise in the solution. It is a cell-centered numerical scheme with extrapolated cell-edge variables. The extrapolated variables are calculated using a linear reconstruction with flux limiters. The nonlinear flux is handled using a high-resolution scheme by Kurganov and Tadmor (2000) of the form:

[TABLE]

[TABLE]

where $u$ is the primitive variable, $f$ is the nonlinear conservative flux, $l$ and $r$ are the left and right extrapolated states, $j$ is the index of one of the two dimensions, and $\lambda$ is the largest eigenvalue (characteristic speed) of the corresponding $u$.

The left and right states are calculated using the following linear reconstruction:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\phi$ is the flux limiter and the simple minmod flux limiter (Roe, 1986) is chosen for simplicity and speed. The flux limiter ensures the scheme is total-variation-diminishing. With the method above the numerical scheme is second-order when the solution is smooth and first-order upwind when a shock is present.

2.4 Semi-Implicit Time Integration

The time integration used for the simulation is the simple semi-implicit Eulerian time step.

[TABLE]

The above equation is split into an implicit left side and explicit right side. Here $R_{\text{I}}$ are the implicit terms and consists of the collisional and diffusive terms, like collision rate, resistivity and heat flux. The collision time scales are typically many orders of magnitude shorter than Alfvén time scales so these terms are dominant. The momentum equations (Equations \irefeq:plasmamomentum and \irefeq:neutralmomentum) have collisional terms for each fluid which introduce off-diagonal terms that must be solved implicitly. While the case studied is the chromospheric network, which has sharp magnetic gradients those off-hand terms are also large and this makes a fully implicit code ill-conditioned. This means a linear system is solved in each time step. $R_{\text{E}}$ contains the explicit terms, which are the source terms and the numerical fluxes calculated as shown above.

The simulation spans hundreds of thousands of time steps, which is why a simple time integration method is used. Conversely, it can be shown using Von Neumann analysis that the general Lagrange multiplier method used in Equation \irefeq:glm is unstable in other semi-implicit methods.

In order to satisfy Courant-Friedrichs-Lewy (CFL) condition is satisfied by finding the fast wave speed and local bulk flow velocities:

[TABLE]

Where $c$ are the characteristic speeds. This is calculated for both plasma and neutrals (where neutral’s is simply the sound speed). Then with a CFL of 0.1 the time steps become:

[TABLE]

Here the maximum of the local speed and fast wave speed sum for each cell is considered. Although the scheme is stable under a higher CFL, a CFL of 0.1 was chosen for accuracy.

2.5 Code Architecture

\ilabel

sec:code

The code is written in C++ to make prototyping easy and still have low-level access. The Armadillo library (Sanderson and Curtin, 2016) is used for matrix building and CPU related code as it makes writing linear algebra code more expressive while still calling BLAS subroutines directly. To make use of GPUs while still being brand-agnostic, the GPU library ViennaCL (Rupp, Tillet, and Rudolf, 2016) is used. The library is written in OpenCL and provides iterative solvers and preconditioners, which the code takes advantage of due to its semi-implicit nature. ViennaCL provides high and low-level access to the GPU. Low-level access is needed to make a custom OpenCL kernel of the flux limiter, and the minmod function shown in \irefminmod.

Convergence tests show that the results are consistent with larger sizes (which the GPU’s memory can hold). This grid size is chosen because it is the smallest size that does not require a trade-off in terms of the physics captured by the model.

3 Code Validation

\ilabel

sec:validation

To validate the code, two tests have been run to ensure that it is producing accurate results. The first is for the neutral fluid, which does not experience Lorentz forces on its own, but does experience pressure forces. The Sod’s shock-tube test is used here to demonstrate the code’s shock capturing capabilities. The second is the initial value problem tests done by Soler et al. (2013) to show Alfv́en-wave propagation in a partially ionized collisional medium. This is to demonstrate the collisions and the electro-magnetic forces working as they are supposed to.

3.1 Sod’s shock-tube

The shock-tube tests studied by Sod (1978) are often used to tests against shock capturing schemes. The shock capturing scheme used in this code is the TVD-MUSCL scheme which is known to perform well against these tests. The fluid used in this case is the neutral fluid in the absence of collisions. The thermal-energy Equation \irefeq:energyneutral has been modified to a pressure equation similar to Sod’s to ensure reproducible results; every other equation has not been changed.

Figure \ireffig:shocktube shows the results, and they mimic the results shown in the article by Sod (1978). The runs were done in arbitrary units. The initial conditions are $\rho_{\text{l}}=1$, $V_{\text{l}}=0$, $p_{\text{l}}=1$ and $\rho_{\text{r}}=0.125$, $V_{\text{r}}=0$, $p_{\text{l}}=0.1$. The number of grids chosen were 1001 and the results are taken at time $t=0.1$. The boundary-conditions were all static Dirichlet boundary-conditions matching the initial conditions for either side. The five regions mentioned in the original article by Sod (1978) can be clearly seen in the density plot.

3.2 Initial-Value Problem

Soler et al. (2013) ran many tests of Alfvén-wave propagation in a partially ionized medium. The time-dependent test that was presented was an initial-value-problem where the initial value of the velocity is $v_{\perp}(z)=\cos{(20\pi z/L)}$ in the region of $z=[-10L,10L]$ corresponding to a fundamental mode of a standing Alfvén-wave. Three tests are shown of two different cases, where the initial value of the neutral velocity is similar or zero. The three tests are three different orders of magnitude of the collision rates. It is important to note the MolMHD code used in this article is a fundamentally different set of equations, using the primitive variables instead of the conservative ones. However, despite these differences, the results shown in Figure \ireffig:soler are similar overall.

4 The Simulation Domain

\ilabel

sec:domain

In the horizontal dimension, we simulate a supergranule similar to that described by Song (2017). Because of symmetry, our simulation domain contains half of a supergranule, from the center of a network with strong-field to the center of internetwork at the photospheric boundary. If we further assume that the neighboring supergranules have the same polarity as the one in the simulation domain, there is no reconnection between the supergranules of interest. In this case the normal component of the field [$B_{x}$] is zero at the left and right boundaries. The right and left boundaries act as mirrors to reflect any normal component of the field. Similarly, the normal component of the velocity is also set to zero at the right and left boundaries, i.e. there is no horizontal flow from one supergranule to another at the chromospheric altitude although inter-supergranule flow is possible at coronal heights and in or below the photosphere.

Since we only simulate the chromosphere and treat the transition region above and photosphere below as upper and lower boundaries, respectively, the photosphere below is made up as a Dirichlet boundary-conditions. Below the surface of the photosphere, the Sun is highly dynamic. The assumption made in this model is that the Alfvén waves are generated at the photosphere and propagate while undergoing damping towards the transition region. For this reason, to imitate the energy flux and activity of the photosphere, a driver at the bottom of the simulation domain is made as a Dirichlet boundary-conditions. Using the broadband-spectrum-driving method of Tu and Song (2013), the perpendicular velocity is a continuous oscillation with a broadband time series of the form:

[TABLE]

where $V_{\perp}$ will be the imposed speed on the lower boundary in the direction normal to the $x$–$z$ plane. $V_{\text{b}}$ is the amplitude of the broadband spectrum, 2 km s*-1*, to ensure an appropriate photospheric flux of $2\times 10^{7}$ erg cm*-2* s*-1* (see Tu and Song, 2013, section 3.2), $\omega_{j}$ and $\omega_{k}$ are angular frequencies, $\omega_{\text{p}}$ is the peak angular frequency at 1.65 mHz. The driver also has random phases $\Phi_{j},\Phi_{k}$ each of which is a multiple of 2$\pi\times$ a random number between 0 and $s+m$. Here $s$ is the number of frequencies chosen smaller than the peak frequency and $m$ is the number of frequencies larger than the peak frequency. For the case of this simulation $s=100$ and $m=500$.

The number-density of the ions and neutrals are specified at the boundary as the photospheric values from Avrett and Loeser (2008) of 1020 m*-3* and 1023 m*-3*, respectively. As an initial condition, a stratified hydrostatic profile is made for the ion and neutral densities $n(z)=n(z=0)e^{-\text{m}_{\text{i}}\text{g}z/\text{k}_{b}T_{0}}$; this can be seen in Figure \ireffig:ic. Here $T_{0}$ = 7000 K which is the temperature imposed on the photosphere and also the initial condition for both ionized plasma and neutrals in the entire domain. An isothermal approach was chosen initially because the chromosphere’s temperature does not change much, between 6000 K to 8000 K, between the photosphere and the transition region. At later times, the temperatures are determined from the energy Equations \irefeq:energyion and \irefeq:energyneutral.

A requirement in some MHD codes with zero magnetic-field divergence constraints is that the initial condition must also be divergence free. For this reason $B_{x}(t=0)$ and $B_{\perp}(t=0)$ is chosen to be zero while $B_{z}(x,t=0)$ varies in the horizontal $x$-direction. It consists of regions of strong-field, $\approx$1 kG, for the network region and weak-field, $\approx$50 G, for the internetwork region with a sharp linear drop in between. The 1 kG region only takes up 0.05 Mm in the horizontal direction and is in the left part of the domain, while the weak region part is in the rest of the domain. The goal is to allow this initial condition to quickly relax to a canopy because of the large $\beta$ change due to height without solving a force-free field. At the bottom boundary, a Dirichlet condition is applied to $B_{z}(x)$ but for the other two components the field is set to Neumann conditions, $\frac{\partial B_{z,\perp}}{\partial z}\vline_{z=0\,\text{km}}=0$. This configuration resembles the stem of the canopy and the magnitude of the supergranule strengths at the photosphere. The bottom boundary’s $B_{z}(x)$ is an unchanging Dirichlet condition to mimic the convection cell underneath. This is also shown in Figure \ireffig:ic.

At the top boundary which represents the transition region and the corona above it, the model uses Neumann conditions to allow free-flow of the plasma through it (i.e. $\frac{\partial\textbf{u}}{\partial z}\vline_{z=2100\,\text{km}}=0$).

In the simulations and results shown in this article the domain is 15 Mm horizontally (half the size of the magnetic internetwork) and 2.1 Mm vertically (typical thickness of the chromosphere). It is a uniform Cartesian grid of 256$\times$128 where $\Delta x=58.6$ km and $\Delta z=16.4$ km. This is the largest cell size chosen for speed and memory considering the usage of a dedicated scientific computing GPU (see Section \irefsec:code). This is larger than the length scale of the collisionally reduced Alfvén-wavelength ($\frac{V_{\text{A}}}{\nu_{\text{in}}}$) which means for this simulation the plasma and neutrals are coupled. Since this is a 2.5D simulation all spatial derivatives in the perpendicular direction are taken as zero ($\frac{\partial}{\partial\perp}=0$).

All momentum components are set initially as zero.

5 Results

\ilabel

sec:results

Figure \ireffig:granulegeometry shows a snapshot on the chromospheric solution after 348 ms. It shows that our simulation produces a magnetic-field topology that is characterized by a canopy or “wine glass” shape (Rutten, 2007), with a strong-field of about 1 kG at the stem on the left, and a weaker field of about 50 G near the upper boundary of the simulation domain, near the transition region. This topology is obtained fairly quickly, in around $\approx$350 ms, because of the magnetic pressure imbalance imposed in the initial conditions. This is not a relaxed system, as the canopy holds for a few seconds before an eruption starts forming.

Figure \ireffig:causality shows the force that is leading to the jet forces. It is a snapshot of $t=101$ ms, which is after the simulation start and before the jet formation.

Figure \ireffig:spiculeshock shows a formation of what resembles a jet–a parcel of plasma quickly moving upward towards the transition region. The shock spontaneously forms after around 20 seconds and the jet speeds up and shoots up outside of the chromosphere into the corona in about five seconds.

In Figure \ireffig:prespicule, we show the distribution of the two heating terms–the ohmic-heating [$\eta J^{2}$] and the frictional-heating [$\nu_{\text{in}}\rho_{\text{i}}|{\mathbfit V}-{\mathbfit U}|^{2}$] in the simulation domain, right before the eruption. The heating rates are normalized per neutral particle. It can be seen that the frictional-heating dominates the ohmic-heating by orders of magnitude, especially in the strong-field region where the shock forms.

Figure \ireffig:midspicule shows the neutral temperature, the frictional-heating, and the velocity difference [$|{\mathbfit V}-{\mathbfit U}|^{2}$] during the emergence of the jet to the corona (at $t=26.8$ seconds). The plots show the formation of an appreciable shock, which separates a low number-density, high temperature, and strong heating region from a high number-density, lower temperature, and weaker heating in the rest of the medium. The figure also shows a very clear difference in heating rates on the two sides of the shock, where the front of the shock travels at super sonic speeds. The shock propagates supersonically until the plasma is ejected when the chromosphere is depleted and starts to fill up again with flows going back down. Downstream of the shock, the plasma flow is returning due to a cavity left behind, where upstream of the shock, the chromospheric plasma is still unchanged.

Figure \ireffig:spiculetimeseries shows a time series of the evolution of the shock number-density and magnetic-field lines. It can be seen that the field lines follow the shock evolution, which means magnetic-field topology (i.e. reconnection region) is not a driver of the jet. Figure \ireffig:heatingtimeseries shows the same time series but with contours of frictional-heating and flow streamlines. It clearly shows that the jet is driven by frictional-heating.

Figures \ireffig:shockwavefront1d, \ireffig:spiculefield1d, \ireffig:spiculefriction1d, and \ireffig:spiculespeed1d show one-dimensional extractions of the solution during the shock emergence along a vertical column at $x=0.6$ Mm. These plots exemplify the cause and effect of the shock. The first plot shows the time series of the wavefront of the shock climbing in the chromosphere. The second and third show that the frictional-heating jump follows the shock, and that the magnetic-field changes correspond to the jet behind the shock. Figure \ireffig:spiculespeed1d shows a velocity profile of the jet.

5.1 Discussion

A magnetic canopy forms due to the magnetic and thermal-pressure imbalance of the strong-field regions from the photosphere and the weak-field but high-temperature regions scaling up to the transition region. This is consistent with a two-fluid analytic model that was presented by Song (2017). It shows that the model can provide numerical solutions for chromospheric conditions. However, the magnetic canopy does not relax or reach a steady state due to the nature of the collisional model and not in the timescales of the current simulation. The supergranule canopy may not last because the photospheric power flux, generated at the lower boundary, continuously adds energy into the system. Small scale processes or highly dynamic events like jets override the canopy magnetic-field because of large plasma speeds. Our results show that the jet has immediate effects on the magnetic-field due to induction and not vice versa.

The largest Alfvén speeds seen in the simulation is 11 Mm s*-1* which means that in the 30 seconds of simulation, magnetic effects can travel through the domain at least three times with constant density and even less in a hydrostatic density. Alfvén waves generated as the drivers for frictional-heating at the bottom boundary do not reflect at the top boundary enough times to warrant a time-averaged study of the entire heating of the chromosphere, nor does it have enough energy to drive the jet of the same kind as presented. The presence of the shock also inhibits any long-term ideal studies of the atmosphere since it depletes and cavities form with the jet. In a future study, we will see how Alfv́en waves affect the chromosphere with force-free initial conditions so as to not form any jets that can destroy a canopy configuration.

Although the frictional-heating is large in the entire magnetic network column, an instance of no frictional-heating creates a separation of upstream and downstream fluids preceding the shock. This can be seen in Figure \ireffig:prespicule and \ireffig:midspicule and is due to a high temperatures and density, which causes the two fluids to move in unison and therefore $\rho_{\text{in}}\nu_{\text{in}}\|V-U\|=0$. The cause of the frictional-heating is not Alfvén-wave damping but the imbalanced initial conditions and the horizontal ${\mathbfit J}\times{\mathbfit B}$ force. The lower medium heats faster than above it, and this creates the shock in the lower atmosphere that then propagates due to thermal-pressure.

It has been suggested that the driver of these kinds of spicule-like events is magnetic reconnection that originates from the vicinity of a large magnetic flux concentration in the network (De Pontieu et al., 2007). In our simulation, the frictional-heating heats a fluid downstream of the unity speed ring larger than that of upstream causing a shock wave as expected from a Type I spicule. However the simulation does not show reconnection effects with enough energy to drive a spicule. Although MHD simulations do not mimic reconnection directly, currents present in reconnection events would accelerate electrons into a collision due to the dense and hot medium surrounding it. Therefore, magnetic reconnection is directly related to ohmic-heating in the chromosphere. If comparing the dominating heating terms, frictional-heating dominates ohmic-heating, which means collisions are causing these jets to form and not reconnection-like currents, then the currents are not due to anti-parallel components but a gradient of a magnetic-field line towards the same direction. The size and extreme temperatures of the jet in the simulation closely resemble chromospheric jets (Shimojo and Shibata, 2000; Nishizuka et al., 2011).

In Figure \ireffig:shockwavefront1d the jet can be seen accumulating a shock of number-density $10^{11}$ cm*-3*. This density is consistent with IRIS observations of spicules in the corona (Alissandrakis et al., 2018). The simulation does show that changes in the magnetic-field happen behind the shock wave front’s direction of motion. Figure \ireffig:spiculefield1d implies the change in the magnetic-field is a consequence of the jet and the rising magnetic-field density in the plane of the simulation is due to compression from the shock itself. The shock in the frictional-heating term [$\rho_{\text{i}}\nu_{\text{in}}\|{\mathbfit V}-{\mathbfit U}\|^{2}/n_{\text{n}}$] seen in Figure \ireffig:spiculefriction1d coincides exactly with the shock showing enhanced heating behind the shock front. The cause of the shock traveling upward is the thermal-pressure behind and in front of the shock.

The jet’s temperatures and speed are much larger than observed on the Sun. This is mostly due to the fact that in our model we have not included the radiative loss in the energy equation. The larger-than-reality shock heating because of the greater flow speed may also further amplify the temperature. The over-estimated temperature of the downstream medium may also lead to larger collision frequencies. This growing downstream temperature has a compounding effect on the thermal-pressure gradient downstream and upstream and causes the speeds to go supersonic. Further studies will incorporate the photo-chemistry and also radiative losses. These implementations will help to see how they affect shock formation and the energy balance inside the chromosphere. It is also important to note that this shock travels at super-Alfvénic speeds.

Kuźma et al. (2017) carried out a numerical simulation of a solar spicule by imposing a vertical pulse at the solar limb. In contrast, the shock shown in our simulation self-consistently evolves without requiring a push on the plasma, and the jet is a consequence of thermal effects. Thus, it provides a complete picture of the chromospheric-heating conditions for shock formation. The spicule speeds in Kuźma et al. (2017) match more closely to the observed 25 km s*-1* than in our simulation, but that is to be expected in the case where the pulse is controlled and the simulation is single-fluid. Kuźma et al. (2017) also investigated adiabatic and non-adiabatic studies of spicules and the wave structures found in the plasma. The jets in our simulation are different from the imposed pulses of Kuźma et al. (2017) since that is a cold dense plasma that has been imposed, while in our simulation, it is a hot large plasma jet that evolved self-consistently.

Song and Vasyliūnas (2014) found that by subtracting the two fluids’ momentum equations the time evolution of the speed differences [$\text{d/dt}|{\mathbfit V}-{\mathbfit U}|$] which leads to frictional-heating, is caused and maintained by magnetic pressure and thermal-pressure differences between the fluids. This is the cause for the heating near the photosphere as shown in Figure \ireffig:causality. The Lorentz force in the horizontal direction dominates other forces and causes a difference in the species’ horizontal speeds. This leads to the frictional-heating that creates a vertical thermal-pressure difference that causes the jet’s upwards motion. Tu and Song (2013), who studied a time-averaged heating in the chromosphere, found that ohmic-heating dominated in the lower chromosphere and frictional-heating dominated in the upper chromosphere, which adequately explained the heating of the chromosphere. This picture seems to be different from the results of our simulation, as frictional-heating dominated ohmic-heating in the lower chromosphere, but this is most likely due to the field strength used in the two studies. Their study limited the field to less than 100 G, but our model includes a region of 1 kG field. The one order of magnitude higher field decreases the ohmic-to-frictional-heating ratio by two orders of magnitude.

We use our model to demonstrate the importance of collisions in the chromosphere. A model without collisions is incomplete as neutral particles dominate the charged particles in this small layer of the Sun. By simply adding collisional terms we show that the simulation replicates chromospheric conditions well. Such a model can be used as an interface between the (fully) neutral photosphere and the (fully) ionized corona. Hillier, Takasao, and Nakamura (2016) investigated 1D slow-mode shocks in chromospheric conditions and how magnetic energy can be converted to frictional-heating in these cases.

We plan to dedicate a future study to investigate how the initial and boundary conditions affect the frictional-heating and the formation of the eruptive feature. We would like to apply periodic boundary-conditions and turn off the Alfvén-wave driver, while shortening the domain to just the lower chromosphere to study the formation more carefully and reduce boundary-conditions effects.

Finally, this work serves as a step towards a better understanding of the solar chromosphere. The Parker Solar Probe mission (Fox et al., 2016) is expected to shed light on the dominant processes in the solar atmosphere, and these processes will be implemented in the model.

6 Conclusion

\ilabel

sec:conclusion

In this article, we present a self-consistent, two-fluid MHD model to study the chromosphere. The model includes ion–neutral interactions, which allows us to investigate the role of frictional-heating in the dynamics of the chromosphere. The very high collision rates in the partially ionized chromosphere stabilizes the code so that no numerical viscosity is needed.

The model obtains preliminary results which require further investigation. Our simulation suggests that, between network and internetwork regions, a jet could potentially be driven purely by thermal processes, especially the frictional-heating, and not by magnetic reconnection. Thus, it is crucial to include the neutrals in any attempt to model the lower chromosphere.

The discrepancy between observational speeds and speeds in the simulation could be explained by the need to include the effects of photo chemical reactions in the chromosphere and the radiative losses. Further refinement of the model with photo-chemical effects and radiation.

Acknowledgments

We thank an unknown referee for their useful comments and suggestions. O. Cohen was supported by NASA Living with a Star grant number NNX16AC11G. J. Tu and P. Song were supported by NSF grants AGS-1702134 and AGS-1559717.