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

Qusai Al Shidi  Ofer Cohen  Paul Song  Jiannan Tu  Lowell Center for Space
Science and Technology, University of Massachusetts Lowell, 600 Suffolk St.,
Lowell, MA 01854, USA

###### keywords:

Chromosphere, Models; Jets

##  1 Introduction

The base of the solar atmosphere (right above the photosphere), referred to as
the chromosphere, is highly dynamic and is strongly influenced by magnetic
field activities (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
networks. It is believed that the network regions correspond to strong
magnetic field at the photosphere altitude. 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_{\\_}B$, plus thermal pressure, $p_{\\_}T$, i.e., $p_{\\_}T+p_{\\_}B$, is
balanced at each height (Gabriel, 1976). This geometry is reflected in a
significant altitudinal change in the plasma beta, which is the ratio of the
thermal to magnetic pressures, $\beta=p_{\\_}T/p_{\\_}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, known as the coronal heating problem and
it is a long-standing problem in solar physics.

In the chromosphere shocks 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 is not well
known. Spicules are categorized into two types (de Pontieu et al., 2007). Type
I spicules have longer lifetimes and are typically larger (5 min and $\sim$15
Mm). Type II spicules have shorter lifetimes and are smaller ($\leq$2 min and
$\sim$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 (Mathioudakis, Jess, and Erdélyi, 2013; McIntosh et al.,
2011; Jess et al., 2009; Ofman and Wang, 2008; Aschwanden et al., 1999;
Nakariakov et al., 1999). These observations of the solar limbs show Alfvén
waves with amplitudes in the range of 5-25 km/s and phase speeds 250-500 km/s.
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 causes 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. 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 (Kuźma et al., 2017; Ni et al., 2016; Brady and Arber, 2016; Ni et
al., 2016; Gudiksen et al., 2011). 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 (Martinez-Sykora, 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 propogation using a collisional two-fluid approach as
well. Three-fluid simulations, with electrons being the third fluid, were made
to study reconnection events in the chromosphere (Leake, Lukin, and Linton,
2013). In these highly localized small regions electron motion is seperate
from ion and neutral motion and its effects cannot be neglected. The
previously mentioned work uses the Bifrost (Gudiksen et al., 2011) code which
employs a very large artificial viscosity for stability. 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).

The aim for this paper 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 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 sec:model the details of the model and code are explained. It is
then followed by section sec: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

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 the 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 three-fluid HiFi code used for studying reconnetion in
the chromosphere (Leake et al., 2012). Except 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-states 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:

$\frac{\partial\rho_{\\_}i}{\partial t}+\nabla\cdot\rho_{\\_}i\textbf{V}=0$
(1)

$\frac{\partial\rho_{\\_}n}{\partial t}+\nabla\cdot\rho_{\\_}n\textbf{U}=0$
(2)

$\ilabel{eq:plasmamomentum}\frac{\partial\rho_{\\_}i\textbf{V}}{\partial
t}=-\nabla\cdot[\rho_{\\_}i\textbf{VV}+n_{\\_}ik_{\\_}bT_{\\_}i+\frac{B^{2}}{2\mu_{\\_}0}\textbf{I}-\frac{\textbf{BB}}{\mu_{\\_}0}]+\rho_{\\_}i\textbf{g}-\rho_{\\_}i\nu_{\\_}{in}(\textbf{V}-\textbf{U})$
(3)

$\ilabel{eq:neutralmomentum}\frac{\partial\rho_{\\_}n\textbf{U}}{\partial
t}=-\nabla\cdot[\rho_{\\_}n\textbf{UU}+n_{\\_}nk_{\\_}bT_{\\_}n]+\rho_{\\_}n\textbf{g}+\rho_{\\_}i\nu_{\\_}{in}(\textbf{V}-\textbf{U})$
(4)

$\frac{\partial\textbf{B}}{\partial
t}=-\nabla\cdot(\textbf{VB}-\textbf{BV})-\nabla\times(\eta\textbf{J}+\frac{1}{n_{\\_}ee}\textbf{J}\times\textbf{B})$
(5)

$\ilabel{eq:energyion}\begin{aligned} \frac{\partial
k_{\\_}bT_{\\_}i}{\partial
t}=-\nabla\cdot(k_{\\_}bT_{\\_}i\textbf{V})+\frac{1}{3}k_{\\_}bT_{\\_}i\nabla\cdot\textbf{V}-\frac{2}{3}\frac{\nabla\cdot(\textbf{q}_{\\_}i+\eta\textbf{B}\times\textbf{J})}{n_{\\_}i}\\\
+\frac{2}{3}\frac{m_{\\_}i}{m_{\\_}i+m_{\\_}n}\nu_{\\_}{in}[3k_{\\_}b(T_{\\_}n-T_{\\_}i)+m_{\\_}n|\textbf{V}-\textbf{U}|^{2}]\end{aligned}$
(6) $\ilabel{eq:energyneutral}\begin{aligned} \frac{\partial
k_{\\_}bT_{\\_}n}{\partial
t}=-\nabla\cdot(k_{\\_}bT_{\\_}n\textbf{U})+\frac{1}{3}k_{\\_}bT_{\\_}n\nabla\cdot\textbf{U}-\frac{2}{3}\frac{\nabla\cdot\textbf{q}_{\\_}n}{n_{\\_}n}\\\
+\frac{2}{3}\frac{m_{\\_}n}{m_{\\_}i+m_{\\_}n}\nu_{\\_}{in}[3k_{\\_}b(T_{\\_}i-T_{\\_}n)+m_{\\_}i|\textbf{V}-\textbf{U}|^{2}]\end{aligned}$
(7)

$\nabla\times\textbf{B}=\mu_{\\_}0\textbf{J}$ (8)

$\textbf{q}=-\kappa\nabla T$ (9)

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. V, U, and B are ion velocity, neutral
velocity, and the magnetic field, respectively, with 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=m_{\\_}e(\nu_{\\_}{en}+\nu_{\\_}{ei})/e^{2}n_{\\_}e$ is the
ohmic resistivity. We use the following collision rates (see Pontieu, Martens,
and Hudson, 2001, Appendix A)

$\displaystyle\nu_{\\_}{in}=7.4\times 10^{-11}n_{\\_}nT^{1/2}$
$\displaystyle\nu_{\\_}{en}=1.95\times 10^{-10}n_{\\_}nT^{1/2}$ (10)
$\displaystyle\nu_{\\_}{ei}=3.759\times
10^{-6}n_{\\_}eT_{\\_}e^{-3/2}\ln\Lambda$

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

An in-depth derivation of weakly ionized plasma equations of different forms
can be found in 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 can
be found in 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 divergenceless 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:

$\ilabel{eq:bglm}\frac{\partial\textbf{B}}{\partial
t}=-\nabla\cdot(\textbf{VB}-\textbf{BV})-\nabla\times(\eta\textbf{J}+\frac{1}{n_{\\_}ee}\textbf{J}\times\textbf{B})+\nabla\Phi$
(11) $\ilabel{eq:glm}\frac{\partial\Phi}{\partial
t}+c_{\\_}h^{2}\nabla\cdot\textbf{B}=-\frac{c_{\\_}h^{2}}{c_{\\_}p^{2}}\Phi$
(12)

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

It can be easily shown by taking the divergence of equation eq:bglm that the
introduction of divergence in the magnetic field is correlated with the
multiplier.

$\frac{\partial\nabla\cdot\textbf{B}}{\partial t}=\nabla^{2}\Phi$ (13)

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_{\\_}h=\max(V_{\\_}A)$, and the ratio
$\alpha_{\\_}p=\Delta x\frac{c_{\\_}h}{c_{\\_}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:

$f^{*}_{\\_}{j-1/2}=\frac{1}{2}([f(u^{r}_{\\_}{j-1/2})+f(u^{l}_{\\_}{j-1/2})]-\lambda_{\\_}{j-1/2}[u^{r}_{\\_}{j-1/2}-u^{l}_{\\_}{j-1/2}])$
(14)
$f^{*}_{\\_}{j+1/2}=\frac{1}{2}([f(u^{r}_{\\_}{j+1/2})+f(u^{l}_{\\_}{j+1/2})]-\lambda_{\\_}{j+1/2}[u^{r}_{\\_}{j+1/2}-u^{l}_{\\_}{j+1/2}])$
(15)

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:

$u^{l}_{\\_}{j-1/2}=u_{\\_}{j-1}+\frac{1}{2}\phi(r_{\\_}{j-1})(u_{\\_}j-u_{\\_}{j-1})$
(16)
$u^{l}_{\\_}{j+1/2}=u_{\\_}{j}+\frac{1}{2}\phi(r_{\\_}{j})(u_{\\_}j-u_{\\_}{j-1})$
(17)
$u^{r}_{\\_}{j-1/2}=u_{\\_}{j}+\frac{1}{2}\phi(r_{\\_}{j})(u_{\\_}{j+1}-u_{\\_}{j})$
(18)
$u^{r}_{\\_}{j+1/2}=u_{\\_}{j+1}+\frac{1}{2}\phi(r_{\\_}{j+1})(u_{\\_}{j+1}-u_{\\_}{j})$
(19) $r_{\\_}j=\frac{u_{\\_}j-u_{\\_}{j-1}}{u_{\\_}{j+1}-u_{\\_}j}$ (20)
$\ilabel{minmod}\phi(r)=\max(0,\min(1,r)),\lim_{\\_}{r\to\infty}\phi(r)=1$
(21)

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.

$u^{n+1}+R_{\\_}i(u^{n+1})\Delta t=u^{n}+R_{\\_}e(u^{n})\Delta t$ (22)

The above equation is split into an implicit left side and explicit right
side. Here $R_{\\_}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
eq:plasmamomentum and eq: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 makes a fully-implicit code ill-
conditioned. This means a linear system is solved in each time step.
$R_{\\_}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 eq: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:

$c_{\\_}{fast}=\frac{1}{\sqrt{2}}[c_{\\_}{sound}^{2}+c_{\\_}{Alfven}^{2}+\sqrt{(c_{\\_}{sound}^{2}+c_{\\_}{Alfven}^{2})^{2}-(2c_{\\_}{sound}c_{\\_}{Alfven})^{2}}]^{1/2}$
(23)

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:

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

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 The Simulation Domain

In the horizontal dimension, we simulate a supergranule similar to that
described in 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 in the
corona height 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 mode is 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 condition.
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:

$\displaystyle
V_{\\_}{\perp}(t,z=0)=V_{\\_}b[\sum_{\\_}{k=0}^{s}{(\omega_{\\_}k/\omega_{\\_}p)}^{5/6}\cos{(\omega_{\\_}kt+\Phi_{\\_}k)}$
(25)
$\displaystyle+\sum_{\\_}{j=s+1}^{m}{(\omega_{\\_}p/\omega_{\\_}j)}^{5/6}\cos{(\omega_{\\_}jt+\Phi_{\\_}j)}]$

where $V_{\\_}{\perp}$ will be the imposed speed on the lower boundary in the
direction normal to the $x-z$ plane. $V_{\\_}b$ is the amplitude of the
broadband spectrum, 2 km/s, 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 filled
logarithmically, $\omega_{\\_}p$ is the peak angular frequency at 1.65 mHz.
The driver also has random phases $\Phi_{\\_}j,\Phi_{\\_}k$ which is a
multiple of 2$\pi\times$ by random number between 0 and s+m. 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 densities $n(z)=n(z=0)e^{-m_{\\_}{i}gz/kT_{\\_}0}$ and a
more realistic neutral number density profile is filled based on the above
mentioned model. 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 6000K to
8000 K, between the photosphere and the transition region. At later times, the
temperatures are determined from the energy equations eq:energyion and
eq: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, 1kG, for the network region and weak fieild, 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. The weak region part of 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 Neumann conditions $\frac{\partial
B_{\\_}{z,\perp}}{\partial z}\vline_{\\_}{z=0km}=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.

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=2100km}=0$).

In the simulations and results shown in this paper 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 256x128 where $\Delta x=58.6$ km and $\Delta z=16.4$ km. This is the
largest resolution chosen for speed and memory considering the usage of a
dedicated scientific computing GPU (see sec:code). This is larger than the
length scale of the collisionally reduced Alfvén wavelength
($\frac{V_{\\_}A}{\nu_{\\_}{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.

###  2.6 Code Architecture

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 (Karl Rupp, 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 minmod.

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

##  3 Code Validation

To validate the code two tests have been run to ensure it is producing
accurate results. The first is for the neutral fluid which doesn’t 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 propogation 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 (eq:energyneutral) has been modified to a pressure
equation similar to Sod’s to ensure reproducible results, every other equation
has not been changed.

Figure fig:shocktube shows the results and they mimic the results shown in the
paper. The runs were done in arbitrary units. The initial conditions are
$\rho_{\\_}l=1$, $V_{\\_}l=0$, $p_{\\_}l=1$ and $\rho_{\\_}r=0.125$,
$V_{\\_}r=0$, $p_{\\_}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 paper 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{(\frac{20\pi}{L}z)}$ 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 neutral velocity is
similar or zero. The three tests are three different order of magnitudes of
the collision rates. It is important to note the MolMHD code used in this
paper 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 fig:soler are overall similar.

##  4 Results

sec:results

Figure fig:granulegeometry shows a snapshot on the chromospheric solution
after $348~{}ms$. It shows that our simulation produces magnetic field
topology that is characterized by 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$ at the near the upper boundary of the simulation
domain, near the transition region. This topology is obtained fairly quickly,
in around $\sim$350ms, 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 fig: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 20s and the jet speeds up and shoots up
outside of the chromosphere into the corona in about 5s.

In Figure fig:prespicule, we show the distribution of the two heating terms -
the ohmic heating, $\eta J^{2}$, and the frictional heating,
$\nu_{\\_}{in}\rho_{\\_}i|\textbf{V}-\textbf{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 fig:midspicule shows the neutral temperature, the frictional heating,
and the velocity difference, $|\textbf{V}-\textbf{U}|^{2}$ during the
emergence of the jet to the corona (at $t=26.8~{}s$). The plots show the
formation of an appreciable shock which separates a low number density, high
temperature and strong heating region with a high number density, lower
temperature and weaker heating with the rest of the medium. The figure also
shows a very clear difference in heating rates on both 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 fig: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 fig: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 joule
heating.

Figures fig:shockwavefront1d, fig:spiculefield1d, fig:spiculefriction1d and
fig:spiculespeed1d show one-dimensional extractions of the solution during the
shock emergence along a vertical line 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. The last plot shows a
velocity profile of the jet.

###  4.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 in 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.

Largest Alfvén speeds seen in the simulation is 11 Mm/s which means that in
the 30s of simulation, magnetic effects can travel through the domain at least
3 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
fig:prespicule and fig:midspicule and is due to a high temperatures and
density which causes the two fluids to move in unison and therefore
$\rho_{\\_}{in}\nu_{\\_}{in}\|V-U\|=0$. The cause of the frictional heating is
not Alfvén wave damping but the imbalanced initial conditions and the
horizontal J $\times$ 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 which 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 unison
speed ring larger than that of upstream causing a shockwave 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. 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 fig: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
shockwave front’s direction of motion. Figure fig: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,
$\frac{\rho_{\\_}i\nu_{\\_}{in}\|\textbf{V}-\textbf{U}\|^{2}}{n_{\\_}n}$, seen
in Figure fig:spiculefriction1d coincides exactly at the shock showing
enhanced heating behind the shockfront. The cause of the shock climb is the
thermal pressure behind and in front of the shock.

The jet’s temperatures and speed are much larger than observed. This is mostly
due to the fact that in our model we have not included the radiative loss
effect in the energy equations. 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 estimate of
the 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 so it is a super-Alfvénic jet.

Kuźma et al. (2017) made a numerical simulation of a solar spicule by imposing
a vertical pulse in the solar limb. In contrast, the shock shown in our
simulation self-consistently evolves without forcing 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 our simulation, but that is 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 than 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 evolution of the speed differences,
$d/dt|\textbf{V}-\textbf{U}|$, which leads to frictional heating is caused and
maintained by magnetic pressure and thermal pressure differences between the
fluids. Tu and Song (2013), who studied a time-averaged heating in the
chromosphere, found that ohmic heating dominated the lower chromosphere and
frictional heating dominated 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 in less than 100 G, but our
model includes region of 1 kG field. The one order of magnitude higher in the
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
dominates the charged particles in this small layer of the sun. By simply
adding collision 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, A., Takasao, S.,
and Nakamura, N. (2016) investigated 1D slow-mode shocks in chromospheric
conditions and how magnetic energy can be converted to frictional heating in
these cases.

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,
including the chromosphere, and these processes will be implemented in the
model.

##  5 Conclusion

In this paper, we present a self-consistent, two-fluid MHD model to study the
chromosphere. The model includes ion-neutral interaction, 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 presents 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, radiation and new observations from the
upcoming Parker Solar Probe mission will lead to an improved model of the
chromosphere.

## Acknowledgments

We thank an unknown referee for his/her useful comments and suggestions. OC
were 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.

## References

  * Alissandrakis et al. (2018) Alissandrakis, C.E., Vial, J.-C., Koukras, A., Buchlin, E., Chane-Yook, M.: 2018, IRIS Observations of Spicules and Structures Near the Solar Limb. Sol. Phys. 293, 20\. DOI. ADS.
  * Aschwanden et al. (1999) Aschwanden, M.J., Fletcher, L., Schrijver, C.J., Alexander, D.: 1999, Coronal Loop Oscillations Observed with the Transition Region and Coronal Explorer. ApJ 520, 880\. DOI. ADS.
  * Avrett and Loeser (2008) Avrett, E.H., Loeser, R.: 2008, Models of the solar chromosphere and transition region from sumer and hrts observations: Formation of the extreme-ultraviolet spectrum of hydrogen, carbon, and oxygen. The Astrophysical Journal Supplement Series 175(1), 229\. http://stacks.iop.org/0067-0049/175/i=1/a=229.
  * Brady and Arber (2016) Brady, C.S., Arber, T.D.: 2016, Simulations of alfvén and kink wave driving of the solar chromosphere: Efficient heating and spicule launching. The Astrophysical Journal 829(2), 80\. http://stacks.iop.org/0004-637X/829/i=2/a=80.
  * De Pontieu et al. (2004) De Pontieu, B., Erdelyi, R., De Moortel, I., Metcalf, T.: 2004, Photospheric Oscillations in the Solar Atmosphere: Driving Chromospheric Spicules and Coronal Waves. AGU Fall Meeting Abstracts, SH13A. ADS.
  * de Pontieu et al. (2007) de Pontieu, B., McIntosh, S., Hansteen, V.H., Carlsson, M., Schrijver, C.J., Tarbell, T.D., Title, A.M., Shine, R.A., Suematsu, Y., Tsuneta, S., Katsukawa, Y., Ichimoto, K., Shimizu, T., Nagata, S.: 2007, A Tale of Two Spicules: The Impact of Spicules on the Magnetic Chromosphere. PASJ 59, S655. DOI. ADS.
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., Munz, C.-D., Schnitzer, T., Wesenberg, M.: 2002, Hyperbolic Divergence Cleaning for the MHD Equations. Journal of Computational Physics 175, 645\. DOI. ADS.
  * Fox et al. (2016) Fox, N.J., Velli, M.C., Bale, S.D., Decker, R., Driesman, A., Howard, R.A., Kasper, J.C., Kinnison, J., Kusterer, M., Lario, D., Lockwood, M.K., McComas, D.J., Raouafi, N.E., Szabo, A.: 2016, The Solar Probe Plus Mission: Humanity’s First Visit to Our Star. Space Sci. Rev. 204, 7\. DOI. ADS.
  * Gabriel (1976) Gabriel, A.H.: 1976, A magnetic model of the solar transition region. Philosophical Transactions of the Royal Society of London Series A 281, 339\. DOI. ADS.
  * Gudiksen et al. (2011) Gudiksen, B.V., Carlsson, M., Hansteen, V.H., Hayek, W., Leenaarts, J., Martínez-Sykora, J.: 2011, The stellar atmosphere simulation code Bifrost. Code description and validation. A&A 531, A154. DOI. ADS.
  * Hasan (2008) Hasan, S.S.: 2008, Chromospheric dynamics. Advances in Space Research 42, 86\. DOI. ADS.
  * Hillier, A., Takasao, S., and Nakamura, N. (2016) Hillier, A., Takasao, S., Nakamura, N.: 2016, The formation and evolution of reconnection-driven, slow-mode shocks in a partially ionised plasma. A&A 591, A112. DOI. https://doi.org/10.1051/0004-6361/201628215.
  * Jess et al. (2009) Jess, D.B., Mathioudakis, M., Erdélyi, R., Crockett, P.J., Keenan, F.P., Christian, D.J.: 2009, Alfvén Waves in the Lower Solar Atmosphere. Science 323, 1582\. DOI. ADS.
  * Karl Rupp (2016) Karl Rupp, F.R. Philippe Tillet: 2016, Viennacl - linear algebra library for multi- and many-core architectures. http://viennacl.sourceforge.net/.
  * Khomenko (2017) Khomenko, E.: 2017, On the effects of ion-neutral interactions in solar plasmas. Plasma Physics and Controlled Fusion 59(1), 014038\. http://stacks.iop.org/0741-3335/59/i=1/a=014038.
  * Kurganov and Tadmor (2000) Kurganov, A., Tadmor, E.: 2000, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations. Journal of Computational Physics 160(1), 241 . DOI. http://www.sciencedirect.com/science/article/pii/S0021999100964593.
  * Kuźma et al. (2017) Kuźma, B., Murawski, K., Zaqarashvili, T.V., Konkol, P., Mignone, A.: 2017, Numerical simulations of solar spicules: Adiabatic and non-adiabatic studies. A&A 597, A133. DOI. ADS.
  * Kuźma et al. (2017) Kuźma, B., Murawski, K., Kayshap, P., Wójcik, D., Srivastava, A.K., Dwivedi, B.N.: 2017, Two-fluid numerical simulations of solar spicules. The Astrophysical Journal 849(2), 78\. http://stacks.iop.org/0004-637X/849/i=2/a=78.
  * Leake, Lukin, and Linton (2013) Leake, J.E., Lukin, V.S., Linton, M.G.: 2013, Magnetic reconnection in a weakly ionized plasma. Physics of Plasmas 20(6), 061202\. DOI. https://doi.org/10.1063/1.4811140.
  * Leake et al. (2012) Leake, J.E., Lukin, V.S., Linton, M.G., Meier, E.T.: 2012, Multi-fluid simulations of chromospheric magnetic reconnection in a weakly ionized reacting plasma. The Astrophysical Journal 760(2), 109\. http://stacks.iop.org/0004-637X/760/i=2/a=109.
  * Leake et al. (2014) Leake, J.E., DeVore, C.R., Thayer, J.P., Burns, A.G., Crowley, G., Gilbert, H.R., Huba, J.D., Krall, J., Linton, M.G., Lukin, V.S., Wang, W.: 2014, Ionized plasma and neutral gas coupling in the sun’s chromosphere and earth’s ionosphere/thermosphere. Space Science Reviews 184(1), 107\. DOI. https://doi.org/10.1007/s11214-014-0103-1.
  * Maneva et al. (2017) Maneva, Y.G., Laguna, A.A., Lani, A., Poedts, S.: 2017, Multi-fluid modeling of magnetosonic wave propagation in the solar chromosphere: Effects of impact ionization and radiative recombination. The Astrophysical Journal 836(2), 197\. http://stacks.iop.org/0004-637X/836/i=2/a=197.
  * Martinez-Sykora, Pontieu, and Hansteen (2012) Martinez-Sykora, J., Pontieu, B.D., Hansteen, V.: 2012, Two-dimensional radiative magnetohydrodynamic simulations of the importance of partial ionization in the chromosphere. The Astrophysical Journal 753(2), 161\. http://stacks.iop.org/0004-637X/753/i=2/a=161.
  * Martinez-Sykora et al. (2017) Martinez-Sykora, J., Pontieu, B.D., Carlsson, M., Hansteen, V.H., Nóbrega-Siverio, D., Gudiksen, B.V.: 2017, Two-dimensional radiative magnetohydrodynamic simulations of partial ionization in the chromosphere. ii. dynamics and energetics of the low solar atmosphere. The Astrophysical Journal 847(1), 36\. http://stacks.iop.org/0004-637X/847/i=1/a=36.
  * Mathioudakis, Jess, and Erdélyi (2013) Mathioudakis, M., Jess, D.B., Erdélyi, R.: 2013, Alfvén Waves in the Solar Atmosphere. From Theory to Observations. Space Sci. Rev. 175, 1\. DOI. ADS.
  * McIntosh et al. (2011) McIntosh, S.W., de Pontieu, B., Carlsson, M., Hansteen, V., Boerner, P., Goossens, M.: 2011, Alfvénic waves with sufficient energy to power the quiet solar corona and fast solar wind. Nature 475, 477\. DOI. ADS.
  * Mignone, Tzeferacos, and Bodo (2010) Mignone, A., Tzeferacos, P., Bodo, G.: 2010, High-order conservative finite difference GLM-MHD schemes for cell-centered MHD. Journal of Computational Physics 229, 5896\. DOI. ADS.
  * Nakariakov et al. (1999) Nakariakov, V.M., Ofman, L., Deluca, E.E., Roberts, B., Davila, J.M.: 1999, TRACE observation of damped coronal loop oscillations: Implications for coronal heating. Science 285, 862\. DOI. ADS.
  * Ni et al. (2016) Ni, L., Lin, J., Roussev, I.I., Schmieder, B.: 2016, Heating Mechanisms in the Low Solar Atmosphere through Magnetic Reconnection in Current Sheets. ApJ 832, 195\. DOI. ADS.
  * Nishizuka et al. (2011) Nishizuka, N., Nakamura, T., Kawate, T., Singh, K.A.P., Shibata, K.: 2011, Statistical study of chromospheric anemone jets observed with hinode/sot. The Astrophysical Journal 731(1), 43\. http://stacks.iop.org/0004-637X/731/i=1/a=43.
  * Ofman and Wang (2008) Ofman, L., Wang, T.J.: 2008, Hinode observations of transverse waves with flows in coronal loops. A&A 482, L9. DOI. ADS.
  * Pontieu, Martens, and Hudson (2001) Pontieu, B.D., Martens, P.C.H., Hudson, H.S.: 2001, Chromospheric damping of alfvén waves. The Astrophysical Journal 558(2), 859\. http://stacks.iop.org/0004-637X/558/i=2/a=859.
  * Roe (1986) Roe, P.L.: 1986, Characteristic-based schemes for the euler equations. Annual Review of Fluid Mechanics 18(1), 337\. DOI. https://doi.org/10.1146/annurev.fl.18.010186.002005.
  * Rutten (2007) Rutten, R.J.: 2007, Observing the Solar Chromosphere. In: Heinzel, P., Dorotovič, I., Rutten, R.J. (eds.) The Physics of Chromospheric Plasmas, Astronomical Society of the Pacific Conference Series 368, Heinzel. ADS.
  * Sanderson and Curtin (2016) Sanderson, C., Curtin, R.: 2016, Armadillo: a template-based c++ library for linear algebra. http://arma.sourceforge.net/.
  * Shimojo and Shibata (2000) Shimojo, M., Shibata, K.: 2000, Physical parameters of solar x-ray jets. The Astrophysical Journal 542(2), 1100\. http://stacks.iop.org/0004-637X/542/i=2/a=1100.
  * Sod (1978) Sod, G.A.: 1978, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. Journal of Computational Physics 27(1), 1 . DOI. http://www.sciencedirect.com/science/article/pii/0021999178900232.
  * Soler et al. (2013) Soler, R., Carbonell, M., Ballester, J.L., Terradas, J.: 2013, Alfvén Waves in a Partially Ionized Two-fluid Plasma. ApJ 767, 171\. DOI. ADS.
  * Song (2017) Song, P.: 2017, A model of the solar chromosphere: Structure and internal circulation. Astrophysical Journal 846(2). DOI.
  * Song and Vasyliūnas (2011) Song, P., Vasyliūnas, V.M.: 2011, Heating of the solar atmosphere by strong damping of Alfvén waves. Journal of Geophysical Research (Space Physics) 116, A09104. DOI. ADS.
  * Song and Vasyliūnas (2014) Song, P., Vasyliūnas, V.M.: 2014, Effect of horizontally inhomogeneous heating on flow and magnetic field in the chromosphere of the sun. The Astrophysical Journal Letters 796(2), L23. http://stacks.iop.org/2041-8205/796/i=2/a=L23.
  * Tu and Song (2013) Tu, J., Song, P.: 2013, A study of alfvén wave propagation and heating the chromosphere. The Astrophysical Journal 777(1), 53\. http://stacks.iop.org/0004-637X/777/i=1/a=53.
  * Tu and Song (2016) Tu, J., Song, P.: 2016, A two-dimensional global simulation study of inductive-dynamic magnetosphere-ionosphere coupling. Journal of Geophysical Research (Space Physics) 121, 11\. DOI. ADS.
  * Vasyliūnas and Song (2005) Vasyliūnas, V.M., Song, P.: 2005, Meaning of ionospheric Joule heating. Journal of Geophysical Research (Space Physics) 110, A02301. DOI. ADS.
  * Wójcik, Murawski, and Musielak (2018) Wójcik, D., Murawski, K., Musielak, Z.E.: 2018, Acoustic waves in two-fluid solar atmosphere model: cut-off periods, chromospheric cavity, and wave tunnelling. Monthly Notices of the Royal Astronomical Society 481(1), 262\. DOI. http://dx.doi.org/10.1093/mnras/sty2306.

Figure 1.: Neutral fluid shock tube tests. The shocks are being captured
appropriately due to the TVD-MUSCL scheme used.fig:shocktube Figure 2.: The
Soler et al. (2013) initial value problem tests. Velocity is taken at $z=0$ of
the standing Alfv́en wave. The solid red is the ion velocity, the dashed blue
line is neutral velocity. The top row is where the initial neutral velocity is
the same as the ion velocity and the bottom row is where the initial neutral
velocity is at zero. The results resemble reasonably well the tests in the
Soler et al. (2013) paper.fig:soler Figure 3.: Formation of a canopy seen
after around 348ms (13000 time steps). The initial condition of a concentrated
strong field region (1kG) and weak field region (5G) creates vertical
asymmetry as plasma beta drastically decreases. The left and right boundaries
are mirrored so effects of other granules behaving the same way are
considered. Because of the short time period this is not a relaxed system but
simply a consequence of the pressure imbalance frequently seen in the
chromosphere.fig:granulegeometry Figure 4.: A zoom into the shockwave front of
the jet. The arrows represents the direction and magnitude of ion flow. The
color bar represents the neutral number density in cm-3fig:spiculeshock Figure
5.: Heating of the shock origin before the formation of the jet at t=18.9s.
Top image is the frictional heating per neutral particle. Bottom image is the
ohmic heating per neutral particle. The frictional heating is orders of
magnitude larger than the heating due to currents.fig:prespicule Figure 6.:
These figures show the jet in the middle of its journey to the corona at
t=26.8s. From top to bottom, the temperature of the neutral plasma, the
frictional heating per neutral particle, and the difference in speeds
$\|\textbf{V}-\textbf{U}\|^{2}$ of the charged and neutral plasma. The shock
can be clearly seen in the bottom-left corner of the simulation where the
magnetic flux concentration is strongest. Due to large plasma beta a ‘ring’ of
slow moving plasma is formed and persists in the shock allowing the shock to
start forming.fig:midspicule Figure 7.: Time series of the shock. From left to
right and top to bottom the times are t=26.751s, 26.785s, 26.831s, 26.842s,
26.844s, 26.849s. The white streamlines denote the magnetic field lines and
the number density of the neutrals are in cm-3. This shows them magnetic field
geometry changing due to the jet and not the other way round. The jet is not
generated because of the magnetic field.fig:spiculetimeseries Figure 8.: Time
series of frictional heating in the jet. From left to right and top to bottom
the times are t=26.751s, 26.785s, 26.831s, 26.842s, 26.844s, 26.849s. The
white streamlines denote flow.fig:heatingtimeseries Figure 9.: Ion number
density snapshots of shockwave front at line of extraction x=0.6Mm. Ion number
density is in units of $10^{11}$ cm-3. The wavefront can be clearly seen
climbing the Chromosphere.fig:shockwavefront1d Figure 10.: Blue line is ion
number density and the black and gray dashed lines are $B_{\\_}x$ and
$B_{\\_}z$ respectively at t=26.83s and vertical line at x=0.6Mm. Ion number
density is in units of $10^{11}$ cm-3. The magnetic field is changed behind
the shockwave front.fig:spiculefield1d Figure 11.: Blue line is ion number
density and red dashed line is frictional heating per neutral particle at
t=26.83s and vertical line x=0.6Mm. Ion number density is in units of
$10^{11}$ cm-3. The frictional heating per neutral particle behind the
wavefront is orders of magnitude higher than in front of
it.fig:spiculefriction1d Figure 12.: Blue line is ion number density and black
dashed line is vertical ion velocity at t=26.83s and vertical line x=0.6Mm.
Ion number density is in units of $10^{11}$ cm-3. The jet speed is
supersonic.fig:spiculespeed1d