# Numerical simulations of a two-fluid jet at a magnetic null point in a solar
arcade

J. J. González-Avilés1 K. Murawski2 and T. V. Zaqarashvili3,4,5
1Investigadores por México-CONACYT, Servicio de Clima Espacial México,
Laboratorio Nacional de Clima Espacial, Instituto de Geofísica, Unidad
Michoacán,
Universidad Nacional Autónoma de México, 58190 Morelia, Michoacán, México
2 Institute of Physics, University of Maria Curie-Skłodowska, Pl. M. Curie-
Skłodowskiej 5, PL-20-031 Lublin, Poland
3 IGAM, Institute für Physik, University of Graz, Universitätsplatz 5, A-8010
Graz, Austria
4 Ilia State University, Cholokashvili ave 5/3, Tbilisi, Georgia
5 Abastumani Astrophysical Observatory, Mount Kanobili, Abastumani, Georgia
E-mail<EMAIL_ADDRESS>
(Accepted 2022 July 15. Received 2022 June 21; in original form 2021 November
05)

###### Abstract

We study the formation and evolution of jets in the solar atmosphere using
numerical simulations of partially ionized plasma. The two-fluid
magnetohydrodynamic equations with ion+electron and neutral hydrogen
components are used in two-dimensional (2D) Cartesian geometry. Numerical
simulations show that a localized nonlinear Gaussian pulse of ion and neutral
pressures initially launched from the magnetic null point of a potential
arcade located below the transition region quickly develops into a shock due
to the decrease of density with height. The shock propagates upwards into the
solar corona and lifts the cold and dense chromospheric plasma behind in the
form of a collimated jet with an inverted-Y shape. The inverted-Y shape of
jets is connected with the topology of a magnetic null point. The pulse also
excites a nonlinear wake in the chromosphere, which leads to quasi-periodic
secondary shocks. The secondary shocks lift the chromospheric plasma upwards
and create quasi-periodic jets in the lower corona. Ion and neutral fluids
show generally similar behavior, but their relative velocity is higher near
the upper part of jets, which leads to enhanced temperature or heating due to
ion-neutral collisions. Simulations of jets with inverted-Y shape and their
heating may explain the properties of some jets observed in the solar
atmosphere.

###### keywords:

(magnetohydrodynamics) MHD – methods: numerical – Sun: atmosphere – Sun:
chromosphere

††pubyear: 2022††pagerange: Numerical simulations of a two-fluid jet at a
magnetic null point in a solar arcade–Numerical simulations of a two-fluid jet
at a magnetic null point in a solar arcade

##  1 Introduction

Nowadays, the knowledge of the solar atmosphere properties has significantly
improved, but it is still challenging to understand its dynamics and
associated phenomena fully. In particular, the fine and small-scale plasma
structures still pose many open questions to answer (Judge, 2006; Lipartito et
al., 2014; Kiss et al., 2017). Among these fine structures are spicules, which
are spiky-like gas jets mainly observed in the chromosphere (Secchi, 1877;
Sterling, 2000; De Pontieu et al., 2004), and the inverted-Y shape jet events
first identified by Shibata et al. (2007), which could manifest as bright
regions on broadband Ca II H images (Morita et al., 2010; Nishizuka et al.,
2011; Singh et al., 2012). In the context of numerical magnetohydrodynamic
(MHD) simulations of inverted-Y shape jets, Nishizuka et al. (2008) showed the
formation of giant chromospheric anemone jets with the inverted-Y shape, as
observed by Hinode. Yang et al. (2013) found that 2.5D numerical MHD
simulations reproduced many observed chromospheric anemone jets with
inverted-Y shapes. The formation of such jets are commonly related to
signatures of magnetic reconnection in sunspots (Tian et al., 2018), in the
solar photosphere (Nelson et al., 2019), in the upper chromosphere and in the
TR (Nishizuka et al., 2011). The magnetic reconnection may take place, for
instance, near magnetic null points, which are expected to be rather common in
different layers of the solar atmosphere (Close et al., 2004; Régnier et al.,
2008; Freed et al., 2015). Among the numerical simulations related to the null
points, there are studies of cool and dense blobs (Murawski et al., 2011a) as
well as MHD wave propagation and shocks (see, e.g., Kuźma et al., 2015; Sabri
et al., 2020; Pennicott & Cally, 2021). Besides, Smirnova et al. (2016)
performed 2D numerical simulations of jet-like structures by launching a
pressure pulse from the null point of a potential magnetic arcade in the
framework of the single-fluid MHD and found jets that mimic some properties of
type I and type II spicules.

The numerical modeling of small-scale jets, including spicules and inverted-Y
shape jets, has been mainly developed by adopting the single-fluid MHD
approach. However, in the lower solar atmosphere, which covers the region from
the photosphere to the upper chromosphere, a significant fraction of neutrals
must be taken into account to describe the atmosphere realistically
(Zaqarashvili et al., 2011a). The effect of neutral atoms can be studied
either in the framework of a single-fluid MHD model with Cowling conductivity
or in the two-fluid MHD approach. In the two-fluid model ions+electrons and
neutral atoms are treated as two separate fluids which are coupled by ion-
neutral collisions (Zaqarashvili et al., 2011a; Soler et al., 2013; Leake &
Linton, 2013; Oliver et al., 2016; Kuźma et al., 2017; Al Shidi et al., 2019;
Wójcik et al., 2018).

In this paper, we extend the model of Smirnova et al. (2016) and focus our
analysis on a small-scale and inverted-Y shape jet (not limited to spicules)
using the two-fluid equations that include ions+electrons and neutrals. In
particular, we simulate the formation of an impulsively generated jet by
pressure pulses initially launched from a magnetic null point of a potential
arcade. We suppose that a trigger of the pressure pulse could be magnetic
reconnection, but we do not describe the details of the excitation. Instead,
we aim to study the formation and evolution of the jet.

We describe in detail the two-fluid equations, the model of the solar
atmosphere, the magnetic field configuration, perturbations, and the numerical
methods in Section 2. In Section 3, we present the results of numerical
simulation and the analysis, and provide the final comments and conclusions in
section 4.

##  2 Model and methods

###  2.1 The system of two-fluid equations

To model the formation and evolution of a jet in the solar atmosphere, we
consider a plasma that consists of two components: ionized fluid
(ions+electrons) and neutral fluid (hydrogen atoms) that interact via ion-
neutral collisions (see, e.g. Draine et al., 1983; Smith & Sakai, 2008;
Zaqarashvili et al., 2011a; Khomenko et al., 2014; Leake et al., 2014; Kuźma
et al., 2017). In particular, we write the two-fluid equations in conservative
form as follows (see, e.g. Oliver et al., 2016):

$\displaystyle\frac{\partial\varrho_{i}}{\partial
t}+\nabla\cdot(\varrho_{i}{\bf v}_{i})$ $\displaystyle=$ $\displaystyle 0,$
(1) $\displaystyle\frac{\partial\varrho_{n}}{\partial
t}+\nabla\cdot(\varrho_{n}{\bf v}_{n})$ $\displaystyle=$ $\displaystyle 0,$
(2) $\displaystyle\frac{\partial(\varrho_{i}{\bf v}_{i})}{\partial
t}+\nabla\cdot\left(\varrho_{i}{\bf v}_{i}{\bf v}_{i}+\left(p_{i}+\frac{{\bf
B}^{2}}{2}\right){\bf I}-{\bf B}{\bf B}\right)$ $\displaystyle=$
$\displaystyle\alpha_{in}({\bf v}_{n}-{\bf v}_{i})+\varrho_{i}{\bf g},$ (3)
$\displaystyle\frac{\partial(\varrho_{n}{\bf v}_{n})}{\partial
t}+\nabla\cdot\left(\varrho_{n}{\bf v}_{n}{\bf v}_{n}+p_{n}{\bf I}\right)$
$\displaystyle=$ $\displaystyle-\alpha_{in}({\bf v}_{n}-{\bf
v}_{i})+\varrho_{n}{\bf g},$ (4) $\displaystyle\frac{\partial E_{i}}{\partial
t}+\nabla\cdot\left(\left(E_{i}+p_{i}+\frac{{\bf B}^{2}}{2}\right){\bf
v}_{i}-{\bf B}({\bf v}_{i}\cdot{\bf B})\right)$ $\displaystyle=$
$\displaystyle\alpha_{in}{\bf v}_{i}\cdot({\bf v}_{n}-{\bf
v}_{i})+\varrho_{i}{\bf g}\cdot{\bf v}_{i}+Q_{i},$ (5)
$\displaystyle\frac{\partial E_{n}}{\partial t}+\nabla\cdot[(E_{n}+p_{n}){\bf
v}_{n}]$ $\displaystyle=$ $\displaystyle\alpha_{in}{\bf v}_{n}\cdot({\bf
v}_{i}-{\bf v}_{n})+\varrho_{n}{\bf g}\cdot{\bf v}_{n}+Q_{n},$ (6)
$\displaystyle\frac{\partial{\bf B}}{\partial t}+\nabla\cdot({\bf v}_{i}{\bf
B}-{\bf B}{\bf v}_{i})$ $\displaystyle=$ $\displaystyle{\bf 0},$ (7)
$\displaystyle\nabla\cdot{\bf B}$ $\displaystyle=$ $\displaystyle 0.$ (8)

Here $\varrho_{i,n}$ are the mass densities of ions and neutrals, ${\bf
v}_{i,n}$ represent the velocities, $p_{i,n}$ are the ion+electron and neutral
gas pressures, $E_{i,n}$ are the total energy densities, $I$ is the unit
matrix tensor, and ${\bf B}$ is the magnetic field which is normalized by
$\sqrt{\mu}$ with $\mu$ being the magnetic permeability. In this paper we
consider the ratio of specific heats $\gamma=5/3$, and the gravitational
source term on the right hand side of equations (3)-(6) is given by ${\bf
g}=(0,-g,0)$ with its magnitude $g=274.78$ m s-2, which represents an average
over the solar surface. In equations (5)-(6), the total energy densities
$E_{i}$ and $E_{n}$ are given as

$\displaystyle E_{i}$ $\displaystyle=$
$\displaystyle\frac{p_{i}}{\gamma-1}+\frac{\varrho_{i}|{\bf
v}_{i}|^{2}}{2}+\frac{{|\bf B|}^{2}}{2},$ (9) $\displaystyle E_{n}$
$\displaystyle=$ $\displaystyle\frac{p_{n}}{\gamma-1}+\frac{\varrho_{n}|{\bf
v}_{n}|^{2}}{2}.$ (10)

We neglect the effect of ionization, recombination, and electron dynamics. The
latter is justified as the electron mass is small compared to the ion and
neutral masses. The effect of the interaction between ions and neutrals
depends on friction coefficient, $\alpha_{in}$, defined by assuming that the
ion-neutral momentum transfer cross-section, $\sigma_{in}$, is independent of
the relative velocity of the colliding particles. This assumption is valid
because charge exchange interactions (for which $\sigma_{in}\approx$ constant)
generally contribute more to the momentum exchange than elastic scattering.
Considering these approximations, $\alpha_{in}$ is defined (Braginskii, 1965)
as

$\alpha_{in}=\frac{4}{3}\frac{\sigma_{in}}{m_{i}+m_{n}}\sqrt{\frac{8k_{B}}{\pi}\left(\frac{T_{i}}{m_{i}}+\frac{T_{n}}{m_{n}}\right)}\varrho_{i}\varrho_{n},$
(11)

where $k_{B}$ is Boltzmann’s constant, $m_{i,n}$ are the atomic masses, and
$T_{i,n}$ are their temperatures (see, e.g. Braginskii, 1965; Chapman &
Cowling, 1970). Regarding the proton-hydrogen collision cross section, we use
$\sigma_{in}=1.4\times 10^{-19}$ m2 (Vranjes & Krstic, 2013). The frictional
interaction between ions and neutrals results in additional heat production
exchange terms, $Q_{i,n}$, in equations (5) and (6). These terms are specified
as (see, e.g., Ballester et al., 2018)

$\displaystyle Q_{i}=\alpha_{in}\left[\frac{1}{2}|{\bf v}_{i}-{\bf
v}_{n}|^{2}+\frac{3}{2}\frac{k_{B}}{m_{H}(\mu_{i}+\mu_{n})}(T_{n}-T_{i})\right],$
(12) $\displaystyle Q_{n}=\alpha_{in}\left[\frac{1}{2}|{\bf v}_{i}-{\bf
v}_{n}|^{2}+\frac{3}{2}\frac{k_{B}}{m_{H}(\mu_{i}+\mu_{n})}(T_{i}-T_{n})\right].$
(13)

Equations (1)-(8) are supplemented by the following ideal gas laws:

$p_{i}=\frac{k_{B}}{m_{i}}\varrho_{i}T_{i},\quad
p_{n}=\frac{k_{B}}{m_{n}}\varrho_{n}T_{n}.$ (14)

Here $m_{i}=m_{H}\mu_{i}$, $m_{n}=m_{H}\mu_{n}$, and $m_{H}$ represents the
hydrogen mass, which is the main ingredient of the gas, and therefore
$m_{n}\simeq m_{i}=m_{H}=m_{p}$ (with $m_{p}$ being the proton mass). The mean
masses $\mu_{i}\approx 0.58$ and $\mu_{n}\approx 1.21$, which include the
influence from heavier elements taken from the OPAL solar abundance table
model (see, e.g. Vögler et al., 2005).

The system of equations (1)-(8) is general and represent a 3D Cartesian
dimensional space, however we limit our model to a two-dimensional (2D) case,
where $y$ axis is directed vertically upwards, with transverse components
$(\partial/\partial z=0)$ of velocities and magnetic field set identically to
zero, i.e., $v_{iz}=v_{nz}=B_{z}=0$. As a result, Alfvén waves are removed
from the consideration.

###  2.2 Model of the solar atmosphere

At the initial time of simulations, we assume that the solar atmosphere is in
magnetohydrostatic equilibrium, i.e., we set the ion and neutral velocities
equal to zero (${\bf v_{i}=v_{n}=0}$). Additionally, the Lorentz force acts
directly on the ions, and it should balance the gravity force and the gas
pressure gradient,

$(\nabla\times{\bf B})\times{\bf B}-\nabla p_{i}+\varrho_{i}{\bf g}={\bf 0},$
(15)

while for neutrals the hydrostatic equation is satisfied,

$-\nabla p_{n}+\varrho_{n}{\bf g}={\bf 0}.$ (16)

Considering the ideal gas laws for ions and neutrals, given by equation (14),
and taking into account the $y$-components of magnetohydrostatic equation (15)
for a force-free magnetic field ($(\nabla\times{\bf B})\times{\bf B}={\bf 0}$)
and equation (16), we arrive to the expressions for the equilibrium gas
pressures

$\displaystyle p_{n}(y)$ $\displaystyle=$ $\displaystyle
p_{0n}\exp{\left(-\int_{y_{0}}^{y}\frac{dy^{\prime}}{\Lambda_{n}(y^{\prime})}\right)},$
(17) $\displaystyle p_{i}(y)$ $\displaystyle=$ $\displaystyle
p_{0i}\exp{\left(-\int_{y_{0}}^{y}\frac{dy^{\prime}}{\Lambda_{i}(y^{\prime})}\right)}.$
(18)

Here

$\displaystyle\Lambda_{i}(y)=\frac{k_{B}T_{i}(y)}{m_{H}\mu_{i}g}\quad\mbox{and}\quad\Lambda_{n}(y)=\frac{k_{B}T_{n}(y)}{m_{H}\mu_{n}g}$
(19)

are the pressure scale heights, and $p_{0_{i,n}}$ denote the gas pressures at
the reference base-level $y_{0}=0$ Mm, i.e., at the photosphere. The pressures
$p_{0_{i,n}}$ are the free parameters of the model. So, we can obtain
different profiles for different values of $p_{0_{i,n}}$. In this paper, we
estimate $p_{0_{i,n}}$ using the number density of ions and neutrals at
$y_{0}=0$ Mm from the C7 model of Avrett & Loeser (2008), with the ions and
neutral number densities $n_{i0}(y=y_{0})=8.397\times 10^{19}$ m-3, and
$n_{n0}(y=y_{0})=1.187\times 10^{23}$ m-3. Then, we use equation (14) and get
$p_{0_{i}}=13.1528$ Pa and $p_{0_{n}}=8912.28$ Pa. This choice gives profiles
that consistently describe the behavior of mass densities from the photosphere
up to the TR that is located at $y\approx 2.1$ Mm. The mass densities are
given by equation (14). Specifically, in this paper, we adopt the semi-
empirical model of Table 26 of Avrett & Loeser (2008) for the temperature
field. We consider that temperatures of ions and neutrals are initially equal
(at $t=0$ s), i.e., they are in thermal equilibrium; we simply set
$T_{i}=T_{n}=T$ and show on top-left of Fig. 1. We display all the equilibrium
profiles including the mass densities and gas pressures for ions and neutrals
and the ionization fraction, $\varrho_{n}/\varrho_{i}$, for a region up to
$y=10$ Mm in Fig. 1. Note that at the bottom of the photosphere ($y=0$ Mm),
the mass density of neutrals ($\sim 10^{-4}$ kg m-3) is higher than the mass
density of ions ($\approx 10^{-7}$ kg m-3), while at higher altitudes ($y>2.1$
Mm) the mass density of ions ($\sim 10^{-12}$ kg m-3) is higher than mas
density of neutrals ($\sim 10^{-14}$ kg m-3). Moreover, gas pressure is higher
for neutrals ($\sim 10^{4}$ Pa) than for ions ($\sim 10$ Pa) at $y=0$ Mm.
However, the ion pressure ($p_{i}\sim 10^{-2}$ Pa) is higher than the neutral
pressure ($p_{n}\sim 10^{-3}$ Pa) at higher altitudes. A similar relation
holds for mass densities, i.e., $\varrho_{i}>\varrho_{n}$.

As a result of the adopted profiles for the mass densities, the ionization
fraction ($\varrho_{n}/\varrho_{i}$) is about $10^{3}$ at $y=0$ Mm. In
contrast, the ionization fraction tends to $10^{-1}$ at coronal heights
($y>2.1$ Mm), so the gas is essentially fully ionized there. Although this
ionization fraction seems to be a bit too high, it could be possible to obtain
a smaller ionization fraction, , e.g. of about $10^{-6}$, in the solar corona
by starting with the profiles of mass densities of ions and neutrals in the C7
model. Therefore, we would obtain different profiles of temperatures for ions
and neutrals. However, this would also imply that $T_{n}$ is only $100$ K at
the TR (see, e.g., Kuźma et al., 2017). Alternatively, we could also modify
the amplitude of reference pressures $p_{0_{i,n}}$ to obtain a lower
ionization fraction in the solar corona, but it could not be easy to maintain
the realistic ionization fraction ($\sim 10^{3}$) in the photosphere.

In summary, by solving the hydrostatic equilibrium, we can be closer to the
strongly ionized limit in the corona, but it would depart us a bit too far
from the partially ionized limit at the photosphere, or vice versa, but not
both simultaneously. Therefore, we decided to keep initially the same
temperatures for ions and neutrals and take the pressure values at the
reference base-level $y=0$ Mm from the C7 model (Avrett & Loeser, 2008).
Despite the ionization fraction $\sim 10^{-1}$ in the low corona, the gas is
not as strongly ionized as in Fig. 1 of Alharbi et al. (2021). These
conditions are still reasonable as it is shown in the two-fluid simulations of
jets and waves in the solar atmosphere (see, e.g., Kuźma et al., 2021;
Murawski et al., 2020; Wójcik et al., 2019; Srivastava et al., 2018).

Figure 1: Equilibrium profiles for temperature (top-left), gas pressures of
ions and neutrals (top-right). Mass densities (bottom-left), and ionization
fraction, $\varrho_{n}/\varrho_{i}$, (bottom-right) vs height $y$.

###  2.3 Magnetic field configuration

In a scenario of the solar atmosphere being in magnetohydrostatic equilibrium
at initial time of the simulations, the magnetic field is force-free
($(\nabla\times{\bf B})\times{\bf B}={\bf 0}$) and it satisfies the
divergence-free condition ($\nabla\cdot{\bf B}=0$). We additionally assume
that the magnetic field is current-free ($\nabla\times{\bf B}={\bf 0}$). We
define the magnetic field components in terms of magnetic flux function,
$A(x,y)$, given as, ${\bf B}=\nabla\times(A{\bf\hat{z}})$, where
${\bf\hat{z}}$ is a unit vector along the $z-$direction. The term magnetic
flux function is commonly used for a 2D magnetic field to link with the
definition of the reconnection rate and, particularly, to the different fluxes
in the partition corresponding to the regions of the plane bounded by
separatrices, which are the global stable and unstable manifolds of hyperbolic
nulls points (see, e.g., Yeates & Hornig, 2011). Therefore, the horizontal,
$B_{x}$, and vertical, $B_{y}$, magnetic field components are:

$\displaystyle B_{x}=-\frac{\partial A(x,y)}{\partial y},\quad
B_{y}=\frac{\partial A(x,y)}{\partial x}.$ (20)

Here, $A(x,y)$ is specified as in Smirnova et al. (2016),

$A(x,y)=B_{0}\frac{x}{(kx)^{2}+(k(y+h))^{2}}-bx,$ (21)

with $B_{0}$ representing the strength of the magnetic field of the arcade,
$h$ is the depth of the point of the magnetic singularity, $k=1/h$ is the
inverse of the depth, and $b$ is the strength of a vertical ambient magnetic
field. Thus, $A(x,y)$ describes the geometry of a single magnetic arcade with
the null point located on the symmetry axis ($x=0$). See Fig. 2 (left).

We use equation (20) to express the magnetic field components in the following
forms:

$\displaystyle B_{x}(x,y)$ $\displaystyle=$
$\displaystyle\frac{2B_{0}x(y+h)}{k^{2}(x^{2}+(y+h)^{2})^{2}},$ (22)
$\displaystyle B_{y}(x,y)$ $\displaystyle=$
$\displaystyle\frac{B_{0}(h-x+y)(h+x+y)}{k^{2}(x^{2}+(y+h)^{2})^{2}}-b.$ (23)

According to Smirnova et al. (2016), for $k=1/h$, $B_{0}=B_{y}(0,0)$ has the
physical meaning of the magnetic field strength at the point $(x=0,y=0)$.
Along $y-$axes (for $x=0$ Mm) we have

$B_{y}(x=0,y)=\frac{B_{0}}{\left(\frac{y}{h}+1\right)^{2}}-b.$ (24)

Hence, the coordinate of the null point $Y$, where $B_{y}(0,Y)=0$, can be
calculated from

$B_{0}=b\left(\frac{Y}{h}+1\right)^{2}.$ (25)

In this paper, we let vary the parameters $h$, $Y$, $b$, and $B_{0}$ and
consider two different values of $y$. The chosen parameters are essentially
the same as in Smirnova et al. (2016).

The plasma $\beta$ parameter helps estimate the ratio between ion and neutral
pressures to magnetic pressure. We define the bulk plasma $\beta$ as follows:

$\beta(x,y)=\frac{p_{i}(y)+p_{n}(y)}{B^{2}(x,y)/2}.$ (26)

Here, the pressures $p_{i,n}(y)$ are given by equations (17)-(18), and
$B^{2}=(B_{x}^{2}+B_{y}^{2})$. The spatial profile of plasma $\beta$ is
displayed on the right panel of Fig. 2, where we observe that $\beta>1$ in the
lower atmosphere (the photosphere and the chromosphere) and at the null point
$\beta\rightarrow\infty$; Otherwise, $\beta<1$ in the solar corona ($y>2.1$
Mm). Such behavior of plasma $\beta$ is consistent to our expectations.

Figure 2: Magnetic field lines for $b=20$ G, $h=0.3$ Mm, and $k=3.33$ Mm-1
(left) and bulk plasma $\beta$ (right) at the initial state ($t=0$ s). In
these plots the location of the null point is at $y=$1.6 Mm.

###  2.4 Perturbations

We perturb the magnetohydrostatic equilibrium atmosphere, described in
Sections 2.2 and 2.3, initially (at $t=0$ s), by localized Gaussian pulses in
ion and neutral gas pressures (see, e.g. Kuźma et al., 2017), given as

$p_{i,n}(x,y,t=0)=p_{i,n}(y)\left[1+A_{p}\exp\left(-\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{w^{2}}\right)\right].$
(27)

Here $A_{p}$ denotes the relative amplitude of the pulses, $(x_{0},y_{0})$ is
their initial position, and $w$ their width. We set and hold fixed the value
of $x_{0}$ to 0.0 Mm and let $y_{0}$ vary according to the position of the
null point $Y$. Therefore, we launch the perturbation from the upper
chromosphere ($y_{0}=1.6$ Mm) and above the TR, i.e., from the low corona
($y_{0}=2.2$ Mm). In our simulations, we set $w=0.3$ Mm and $A_{p}=8.5$, which
are the same values as used by Smirnova et al. (2016). The signals launched
from the null point location may mimic a post-reconnection event. The
perturbations produce energy deviations to the background atmosphere at the
initial time, that we estimate using the equations (5) and (6). Therefore, the
deviations in the total energy density of ions and neutrals to the background
atmosphere are

$\displaystyle E_{i}(x,y,t=0)$ $\displaystyle=$
$\displaystyle\frac{p_{i}(y)}{\gamma-1}+\frac{|{\bf
B}(x,y,t=0)|^{2}}{2}+\frac{p_{i}(y)A_{p}}{\gamma-1}\exp\left(-\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{w^{2}}\right)=E_{ib}+\Delta
E_{i},$ (28) $\displaystyle E_{i}(x,y,t=0)$ $\displaystyle=$
$\displaystyle\frac{p_{i}(y)}{\gamma-1}+\frac{p_{i}(y)A_{p}}{\gamma-1}\exp\left(-\frac{(x-x_{0})^{2}+(y-y_{0})^{2}}{w^{2}}\right)=E_{nb}+\Delta
E_{n}.$ (29)

Here the deviations in energies of ions and neutrals from their background
values $(E_{i,nb})$ are $\Delta E_{i}$, and $\Delta E_{n}$, respectively. The
latter means that the amplitude of the Gaussian pulse determines the energy
deviation generated by the perturbation.

###  2.5 Numerical methods

To solve the two-fluid equations (1)-(8) numerically, we employ the JOANNA
code (Wójcik, 2017). In all simulations, we set the Courant-Friedrichs-Levy
(CFL) number equal to 0.9 and choose the third-order strong stability
preserving Runge-Kutta (SSP-RK3) time integrator (Durran, 2010). Additionally,
we adopt the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann
solver (Miyoshi & Kusano, 2005) in combination with a linear reconstruction
and the minmod limiter. To control numerically the growth of the solenoidal
constraint condition given by equation (8), we use the extended generalized
Lagrange multiplier method (Dedner et al., 2002). This method is robust in low
plasma beta ($\sim 10^{-3}-10^{-2}$) regions, as implied in the solar corona;
see, for example, the right panel of Fig. 2.

We carry out the simulations in the domain $x\in[-5.12,5.12]$, $y\in[0,60]$,
in units of Mm. Here, $y=0$ Mm represents the bottom of the photosphere. In
the numerical simulations, we use a uniform grid within the region $(-5.12\leq
x\leq 5.12)$ Mm $\times$ $(0\leq y\leq 10.24)$ Mm, which is covered by
1024$\times$1024 grid cells, which leads to a spatial resolution of $10$ km.
Above this region, within the zone, $(10.24\leq y\leq 60)$ we implement a
stretched grid along the vertical direction $y$; in particular, we divide this
region into $64$ grid cells whose size grows with $y$. Such stretched grid
helps absorb incoming signals and minimizes reflections from the top boundary
(Srivastava et al., 2017; Kuźma et al., 2017). We impose outflow boundary
conditions at the side edges specified by $x=-5.12$ Mm and $x=5.12$ Mm. We set
all the plasma quantities to their equilibrium values at the bottom and top
boundaries delimited by $y=0$ Mm and $y=60$ Mm.

##  3 Results of numerical simulations

###  3.1 Parameters

The equilibrium magnetic field, defined by equations (22) and (23) and
displayed in Fig. 2, represents a magnetic null point in a potential arcade
that extends up to the solar corona. In the context of the solar atmosphere,
the magnetic arcades typically refer to active region structures; however, in
the quiet Sun, an arcade may represent the conditions corresponding to a
bipolar region of a network magnetic field. The maximum magnetic field
strength takes place at the base of the magnetic arcade and it reaches about
$50$ G, see, for example, Choe & Lee (1996) and Bagwell & Ma (2020). Regarding
the null points on the Sun, some estimations indicate that their number is
roughly one per supergranular cell above the height $1.5$ Mm (see, e.g., Close
et al., 2004; Régnier et al., 2008), and likely more below that height.
However, numerical simulations and extrapolations of photospheric magnetic
fields predict ubiquitous null points in the solar chromosphere and the corona
(Galsgaard & Nordlund, 1997; Longcope, 2005). Based on the above references,
from all the possible combinations of the magnetic field strength and the
location of the null point in the solar atmosphere, we select two illustrative
ones, which are labeled as Run#1 and Run#2, and displayed in Table 1. These
runs represent two magnetic field strengths in the upper range of the typical
values of the arcades; however, the magnitude of the fields is about $100$ G,
which is close to the maximum magnetic field strength at the base of the
magnetic arcade. The chosen parameters indicate two locations of the null
points at the chromosphere ($Y=1.6$ Mm) and at the corona ($Y=2.2$ Mm). In the
following subsections, we describe in detail the results corresponding to
Run#1 and some representative results for Run#2.

|  $h$ (Mm) |  $Y$ (Mm) |  $b$ (G) |  $B_{0}$ (G) |  $y_{0}$ (Mm)
---|---|---|---|---|---
Run#1 | 0.3 | 1.6 | 10 | 401 | 1.6
Run#2 | 0.3 | 2.2 | 10 | 684 | 2.2
Table 1: Parameters used for Run#1 and Run#2.

###  3.2 Morphology of the jets

In Fig. 3, we show the mass density of ions, $\varrho_{i}$, (top panels) and
neutrals $\varrho_{n}$ (bottom panels), expressed in kg m-3 and corresponding
to Run#1. In the top-left, we show a snapshot of $\varrho_{i}$ at $t=30$ s.
Here, we identify the formation of a collimated plasma structure consisting of
two parts: (i) the bottom part that mimics an inverted shaped-V jet of
$\varrho_{i}\sim 10^{-11}$ kg m-3, which is larger than the background
$\varrho_{i}\sim 10^{-13}$ kg m-3 and it reaches a height of approximately 3
Mm; (ii) the top part, which is a low, dense plasma of ions ($\varrho_{i}\sim
10^{-12}$ kg m-3) which propagates upwards up to about $4$ Mm. The physical
mechanism of the jet formation is the following: the pressure pulse initially
launched from the solar chromosphere ($y=1.6$ Mm), propagates upwards, grows
in its amplitude due to the decrease of density with height to conserve
momentum $\varrho_{i}v_{iy}$ $c_{i}$, and steepens into a shock. This shock
propagates upward into the corona, but chromospheric plasma lags behind the
shock front, which rises to about $y=4$ Mm. The chromospheric material is
lifted so that the plasma’s rarefaction resibling behind the shock front,
leading to low pressure. As a result, the pressure gradient force works
against the gravity and pushes the chromospheric material to penetrate the
solar corona (see, e.g., Murawski & Zaqarashvili, 2010). The next snapshot
(top-middle), at $t=300$ s, shows that the chromospheric plasma propagates
upwards and becomes narrow. This behavior is consistent with the fact that
ions follow the magnetic field lines, which are vertical at the corona, as
shown in Fig. 2. In the top-right, the mass density of ions at $t=600$ s
illustrates that the plasma flows downwards because of the gravity action. At
this time, a double structure is visible at the top of the jet. This double
structure appears when the plasma falls back at the central part of the jet
and is forced to move upwards again at the jet boundaries due to the next
shock. At the bottom, we display the mass densities of neutrals for the same
three instants of times ($t=30,300,600$ s) that resemble similar distribution
to $\varrho_{i}$. However, we can indicate some differences. For example, at
$t=30$ s, the inverted-Y shape of the jet is slightly smaller ($\sim 3$ Mm)
than for ions. At $t=600$ s, the distribution of mass density of neutrals is
wider ($\sim 0.5$ mm) than the mass density of ions.

Run#1

Figure 3: Temporal evolution of log($\varrho_{i}(x,y)$) (top), and
log($\varrho_{n}(x,y)$) (bottom) corresponding to Run#1 at $t=30$ s, $300$ s,
and $600$ s from left to right.

In Fig. 4 we display snapshots of $\varrho_{i}$ at three instants of time
$t=30,300,600$ s for Run#2. In this scenario, $\varrho_{n}$ resembles
$\varrho_{i}$, so we decided not to show them. On the left, at $t=30$ s, no
collimated plasma jet formation is discernible, as in Run#1. Instead, we see a
small structure of about $10^{-11}$ kg m-3 that practically remains at its
initial position, i.e., at $y=2.2$ Mm. At $t=30$ s, we can also see that the
neutrals with a smaller value of $\varrho_{n}\approx 10^{-12}$ kg m-3 reach a
height of approximately $8$ Mm. However, this structure does not mimic any
characteristic of a jet. It is rather the response of the low, dense
background atmosphere with $\varrho_{n}\sim 10^{-13}$ kg m-13 to the initial
pressure pulses. That is, in this scenario, the location of the pressure pulse
is over the TR, at $y=2.2$ Mm, and it is present in a low-density region with
$\varrho_{n}\approx 10^{-13}$ kg m-3; these conditions are not suitable for
the formation of any jet as it takes place in Run#1. At $t=300$ s (middle),
the less dense in $\varrho_{n}$ structure becomes elongated, while the denser
part with $\varrho_{n}\sim 10^{-11}$ kg m-3 remains just over $y=2.2$ Mm.
Finally, on the right, at $t=600$ s, the mass density of neutrals is similar
to the time $t=300$ s.

Run#2

Figure 4: Temporal evolution of log($\varrho_{i}(x,y)$) corresponding to Run#2
at $t=30$ s, $300$ s, and $600$ s.

Although we impose outflow boundary conditions at the sides edges $x=\pm 5.12$
Mm, the system’s total mass was not significantly lost. We estimate the
relative percentage of mass loss between the initial and final states of the
simulations, given by $\Delta\varrho_{rel}=\frac{\Delta\varrho_{max-
tf}}{\varrho_{max-bg}}*100\%$, where $\Delta\varrho_{max-tf}$ is the maximum
difference of the mass density between the initial and final states, and
$\varrho_{max-bg}$ is the maximum value of the background density at the final
state of the simulations. In the particular case of Run#1, the relative
percentage of the loss mass of ions was about 0.1%, while it was about 0.01%
for neutrals. Therefore, the relative total mass loss was about 0.11% for
Run#1. For Run#2, the relative percentage of mass loss for ions was 0.01%,
while for neutrals, it was approximately 0.001%; then, the relative total mass
loss for Run#2 was about 0.011%. So practically, the total mass was conserved
in these two simulations.

###  3.3 Time evolution of jets

To show the evolution of the ions and neutrals, and identify their potential
differences of Fig. 3, we plot the temporal evolution of the ionization
fraction $\varrho_{n}/\varrho_{i}$ and the absolute value of vertical velocity
difference $|V_{yi}-V_{yn}|$ along vertical cuts in $y$, in Fig. 5 for Run#1
(left panels) and Run#2 (right panels). In the top-left panel, we show the
ionization fraction for Run#1, from which we infer that the temporal evolution
of ions and neutrals differs one from another since, at lower altitudes
($y\sim 2$ Mm), neutrals are more predominant than ions, while at higher
altitudes ($y>2.1$ Mm) the ions are more abundant in the evolution. In
particular, we identify with the dashed line that ions and neutrals follow a
parabolic path, i.e., there is material moving upwards that later falls back.
This parabolic path is quasi-periodic since the next shocks are coming up from
the chromosphere (see below). In the top-right, for Run#2, we see that the
ionization fraction does not vary in the time since neutrals are predominant
at lower altitudes, while ions consist of the abundant component at higher
altitudes. Meanwhile, in the bottom-left panel for Run#1, the logarithm of the
absolute difference between vertical components of ion and neutral velocities,
$\log|V_{yi}-V_{yn}|$, exhibits some differences between the kinematics of
ions and neutrals. We note that the difference between these vertical
velocities attains a minimum along the parabolic path, which means that the
vertical velocity of ions is close to the vertical velocity of neutrals. It
means that the jet carries ions and neutrals, which are lifted up as a part of
the chromospheric material that lags behind the shock. In contrast, in the
surroundings of the jet, in the solar corona, the vertical velocity of ions is
greater than that of neutrals, as it is expected for a predominant region of
ions. Finally, in the bottom-right for Run#2, the vertical velocity difference
is close to zero at lower altitudes, while it is slightly greater at higher
altitudes, implying that ions are predominant in the solar corona.

Run#1 Run#2

Figure 5: Time-distance plots for the $\log(If(x=0,y))$ (left) and
$\log(|V_{y_{i}}(x=0,y)-V_{y_{n}}(x=0,y)|$ for Run#1 (left panels) and for
Run#2 (right panels).

We can identify the global behavior of the mass density of ions for Run#1 by
the time signatures collected for $x=0$ Mm at different detection points along
the vertical direction $y$, as shown in Fig. 6. This plot essentially reveals
the shock’s presence, represented by the steep gradient in the mass density.
This shock propagates and passes $y=3$ Mm reaching the level $y=8$ Mm. In this
figure, it is clear that the shock is strong just after forming, i.e., at
$y=3$ Mm, while it becomes weak when it reaches upper heights. For example, at
$y=6$ Mm (green line), the steep gradient of ion mass density
($\varrho_{i}\sim 10^{-11}$ kg m-3) at about $100$ s is two orders of
magnitude higher than in the ambient coronal plasma ($\varrho_{i}\sim
10^{-13}$ kg m-3).

Figure 6: Time series of the mass density of ions, $\varrho_{i}$, collected at
the points given by $x=0,y=3,4,5,6,7,8$ Mm for Run#1.

It is interesting to check whether the behavior of the jets is related to the
magnetic null point and the potential arcade. Therefore, we perform an
additional simulation considering a constant and straight magnetic field with
its magnitude equal to $10$ G for the exact parameters of the pressure pulses
corresponding to Run#1, i.e., $A_{p}=8.5$ and $y_{0}=1.6$ Mm. For this
simulation, we only show the time-distance plots of the logarithm of the
ionization fraction, $\varrho_{n}/\varrho_{i}$, and the absolute value of the
vertical velocity difference $|V_{yi}-V_{yn}|$, which are displayed in Fig. 7.
For instance, on the left, the ionization fraction shows that ions and
neutrals are distributed more homogeneously over the parabolic path than the
results shown in the top-left panel of Fig. 5. This behavior is related to the
magnetic field configuration since, in the case of the null point in the
potential arcade, the ions are not guided by the magnetic field from the start
of the simulation; instead, they propagate over a free-magnetic region and
then feel the presence of the magnetic field. On the other hand, for a
vertical straight field, the ions are guided from the start of the simulation,
so the pressure pulses propagate and then steepen into a shock that lifts the
chromospheric material, which contains a large amount of neutrals. However, in
the background density of the solar corona, ions are more abundant than
neutrals. The vertical velocity difference shown in the right panel of Fig. 7,
behaves consistently as the ionization fraction does. We can see that the
vertical velocity difference is near zero over the distribution that defines
the parabolic path, while in the solar corona, the vertical velocity of ions
is greater than neutrals. Then, we can specifically distinguish the following
difference between the scenario with the magnetic null point in the potential
arcade and the vertical straight field: (i) the maximum height reached by the
mass densities of ions and neutrals of $\sim 9$ Mm is greater than the reached
in the vertical magnetic field case $\sim 7$ Mm; (ii) the parabolic path of
the mass density of ions and neutrals for the case of the null point seem to
have a greater period of oscillations $\sim 450$ s than $\sim 350$ s shown in
the vertical magnetic field case; (iii) the ions and neutrals behave more
homogeneously over the parabolic path in the case of a vertical straight field
than in the case of the null point in the potential arcade.

Run#1 for a straight magnetic field

Figure 7: Time-distance plots for the $\log(If(x=0,y))$ (left) and
$\log(|V_{y_{i}}(x=0,y)-V_{y_{n}}(x=0,y)|$ (right) corresponding to Run#1 with
a straight magnetic field configuration.

###  3.4 Time series of ions and neutrals at the null point and the cut-off
period

To identify the localized behavior of the jet at the null point, in Fig. 8, we
show the time series of $\log T_{i}$ (left), and vertical velocities
$V_{y_{i,n}}$ (right) collected at the point $(x,y)=(0,1.6)$ Mm, i.e., at the
null point, corresponding to Run#1. In these plots, we illustrate the behavior
of the neutrals since we are analyzing their behavior at the point
$(x,y)=(0,1.6)$ Mm, located at the lower chromosphere, where the dynamic of
neutrals is important. In $\log T_{i}$, we identify variations with a
periodicity of about 200 s. On the right panel, we see that time signatures
for vertical velocities of ions and neutrals behave alike and show an
oscillatory behavior of about 200 s.

Figure 8: Time signatures of $\log T_{i}$ (left), and $V_{y_{i,n}}$ (right)
collected at the point ($x=0$, $y=1.6$) Mm, corresponding to Run#1.

The plasma above the TR, for instance, at $y=5$ Mm for Run#1, exhibits
oscillations of about 450 s. This feature is initially schematized in the
time-distance plot of the ionization fraction of Fig. 5, and confirmed in the
time-distance diagram of the mass density of ions and the time series of the
vertical velocity for ions and neutrals collected at $(x=0,y=5)$ Mm, as shown
on the top-left and top-right of Fig. 9 up to $t=1800$ s. These oscillations
are due to secondary shocks through the process: excited pulse propagates
upwards and steepens into shock due to the density stratification. This shock
propagates from the chromosphere into the corona and lifts the chromospheric
plasma behind. On the other hand, local chromospheric plasma starts to
oscillate at local acoustic cut-off frequency after the propagation of pulses.
In the linear limit, these oscillations are linear and decay as time
progresses (Kuridze et al., 2009; Pylaev et al., 2017). However, in the
nonlinear limit, these oscillations lead to the excitation of quasi-periodic
or secondary shocks. These shocks lift the chromospheric plasma with a
periodicity, which depends on the initial amplitude of the pulse so that the
stronger pulse leads to a longer periodicity (Zaqarashvili et al., 2011b).

We can compare the quasi-periodic oscillations of the secondary shocks shown
on the top panels of Fig. 9 to the acoustic cut-off period for a non-
isothermal solar atmosphere in the linear regime, which following the work by
Deubner & Gough (1984) is defined as follows:

$P_{ac}=\frac{4\pi\Lambda}{c_{s}\sqrt{1+2\frac{d\Lambda}{dy}}}.$ (30)

Here $c_{s}=\sqrt{\gamma
p/\varrho}=\sqrt{\gamma\left(\varrho_{i}/\mu_{i}+\varrho_{n}/\mu_{n}\right)\frac{k_{B}T}{m_{H}}/\varrho_{i}+\varrho_{n}}$
is the total sound speed, i.e., for ions plus neutrals, and
$\Lambda=\left|p/\frac{dp}{dy}\right|$ is the single fluid scale height, with
$p=p_{i}+p_{n}$ representing the total (ions plus neutrals) gas pressure. In
the bottom panel of Fig. 9, we show the acoustic cut-off period $P_{ac}$ as
function of height at the initial time of the simulations. The cut-off period
is higher for ions than for neutrals, from $y=0$ to $y=2.1$ Mm, covering the
photosphere to the chromosphere region. In particular, the values of $P_{ac}$
range from nearly 250 s at $y=0$ Mm, to around 350 s, from $y=1$ Mm to $y=2$
Mm. The cut-off period for neutrals is smaller remaining about 200 s up to the
top chromosphere ($\sim 2$ Mm). The cut-off period for both ions and neutrals
reaches its minimum value of about 50 s at the TR $\sim y=2.1$ Mm; higher up,
the cut-off exponentially grows with height. In particular, for $y=y_{0}=1.6$
Mm, that is the initial position of pressure pulses for Run#1, $P_{ac}$
attains a value of about $350$ s, as shown in the bottom of Fig. 9. The
detected oscillation period of around 450 s (see the top panels of Fig. 9) is
longer due to the nonlinear consideration as shown by Zaqarashvili et al.
(2011b). Our results are consistent with the classical models of rebound
(secondary) shocks for spicules jet formation (see, e.g., Hollweg, 1982;
Sterling & Hollweg, 1988; Murawski & Zaqarashvili, 2010; Murawski et al.,
2011b).

Figure 9: Top: Time-distance plots of $\log(\varrho_{i}(x=0,y))$ (left) and
time series of $V_{y_{i,n}}$ collected at the point $(x=0,y=5)$ Mm (right).
Bottom: The acoustic cut-off period, $P_{ac}$ for ion-neutral fluid (solid
line), ionized fluid (dashed line), and neutrals (dot dashed line), as
functions of height at initial time.

###  3.5 Temperature and the inverted-Y shape jets

To identify the evolution and behavior of ion temperature, which at $t=0$ s is
equal to the neutral atom temperature, we show the temperature maps for three
different instants of time $t=30,60,90$ s, corresponding to Run#1 in Fig. 10.
On the left panel (corresponding to $t=30$ s), we see the two different parts
of temperature: (i) the top part between the heights 3-4 Mm is made up of hot
plasma of about $10^{5}$ K, which corresponds to the upward propagating first
shock; (ii) the lower part below the height of 3 Mm shows colder ($\sim
10^{4}$ K) inverted-Y shape plasma, which corresponds to the chromospheric
jet. In the middle panel (corresponding to $t=60$ s), the hot pulse propagates
upwards and reaches the height of $y=9$ Mm, i.e., the mean velocity is around
170 km s-1 (similar to local sound speed). At the same time, the cold jet
reaches up to $y=4$ Mm, i.e., the mean velocity is 30 km s-1 (similar to the
speed of type I spicules). Finally, in the right panel ($t=90$ s), the hot
pulse passed the considered heights, while the cold inverted-Y shape jet
reached the height of $y=6$ Mm. It is also important to note that the plasma
above the jet (green light color) is hotter ($\sim 10^{5}$ K) than inside (red
and yellow colors). The temperature increase above the jet could be caused by
the heating associated with ion-neutral collisions. This increase in
temperature might be a novel result, which is connected to the two-fluid
approximation. Indeed, Fig. 5 (bottom-left) shows that $|V_{iy}-V_{in}|$ is
large there, which results in deposition of the thermal energy.

From Fig. 10, one can see that the jet shows inverted-Y shape after $t=60$ s.
This inverted-Y shape is related to the magnetic null point as it is absent in
similar simulations for a straight magnetic field.

Figure 10: Ion temperature, $T_{i}$, maps in Kelvin, at time $t=30$ s (left),
$t=60$ s (middle), and $t=90$ s (right), corresponding to Run#1.

The inverted-Y shape jets have been observed in various layers of the solar
atmosphere and can be related to observational signatures of magnetic
reconnection. For instance, Shibata et al. (2007) found that Hinode
observations show the ubiquitous presence of chromospheric anemone jets
outside sunspots in active regions; these small jets had an inverted-Y shape,
similar to the shape of x-ray anemone jets in the corona. Thus, they concluded
that magnetic reconnection could occur at small spatial scales throughout the
chromosphere. Additionally, Tian et al. (2018) used the 1.6 m Goode Solar
Telescope to detect prevalent reconnection through frequently occurring fine-
scale jets in the H$\alpha$ line wings at light bridges of sunspots. They
found that many jets had an inverted-Y shape, shown by models as typical of
reconnection in a unipolar field environment. Moreover, Nelson et al. (2019)
linked the inverted-Y shaped jets to numerous apparent signatures of magnetic
reconnection in the solar photosphere. They used $H\alpha$ imaging data
sampled by the Swedish Solar Telescope to investigate whether bidirectional
flows can be detected within inverted-Y-shaped jets near the solar limb. They
concluded that magnetic reconnection might cause bidirectional flows within
inverted-Y-shaped jets and could be the driver of surges. Furthermore, in the
context of solar jets and reconnection, Sterling et al. (2015) reported high-
resolution X-ray and extreme-ultraviolet observations of jets that form in
coronal holes at the Sun’s poles. They stated that contrary to the emerging-
flux model, a miniature version of the filament eruptions that initiate
coronal mass ejections drives the jet-producing reconnection.

###  3.6 Dynamics of magnetic field lines

In order to show the dynamics of the magnetic field during the simulations, on
Fig. 11 we display initial and final maps of horizontal ($B_{x}$) and vertical
($B_{y}$) components corresponding to Run#1. The upper panels of the figure
illustrate $B_{x}$ at the initial ($t=0$ s) and final ($t=600$ s) times of
simulation, which indicate that the configuration remains almost unchanged,
with only a slight reconstruction of the field lines near the location of the
null point ($y=1.6$ Mm). The lower panels show that the configuration of
$B_{y}$ also did not change during simulation. This is because the generated
pulse and jet propagate almost along the magnetic field lines and do not
affect the lines themselves. On the other hand, the inverted-Y shape of the
exciting jet is connected with the magnetic field configuration around the
null point.

Figure 11: Maps of the horizontal, $B_{x}$ (top panels) and vertical, $B_{y}$,
(bottom panels) components of the magnetic field, expressed in Gauss and
overlaid by the magnetic field lines at times $t=0$ s (left panels) and
$t=600$ s (right panels).

##  4 Final comments and conclusions

We studied the formation and evolution of jets excited by the pulses in ion
and neutral gas pressures. These pulses were launched from the null point of a
potential magnetic arcade in a partially ionized and gravitationally
stratified plasma. We used the two-fluid approximation in the numerical
simulations, where electrons+ions and neutral atoms are recognized as two
different fluids which are coupled by ion-neutral collisions. The two-fluid
effects become significant near the time scales of ion-neutral collisions,
much higher than $1$ Hz in the chromosphere. Nevertheless, for the processes
with longer time scales in partially ionized plasma, two-fluid and single
fluid (with Cowling diffusion) approximations can be applied without
restriction. Both approximations give essentially the same results for the
longer time-scales (see for example, (Zaqarashvili et al., 2011a)). We chose
the two-fluid approach because it is more appropriate for numerical
simulations.

For the magnetic field configuration, we used the 2D potential arcade model
with a null point, which may result from the possible magnetic reconnection
events in the solar atmosphere. First, we studied two different locations of
the magnetic null points at the chromospheric ($Y=1.6$ Mm) and the coronal
($Y=2.2$ Mm) heights. Then, we launched Gaussian pulses in ion and neutral gas
pressures from the null points and followed their evolution in the stratified
atmosphere. The evolution of the pulses and the whole system crucially depend
on the initial location of the pulses (hence on the location of the null
point). When the pulses were launched from $Y=1.6$ Mm height, the following
process was observed in numerical simulations: the pulses steepened into
shocks, penetrating the corona and propagating with the mean velocity of
coronal sound speed. The chromospheric cold and dense plasma followed the
shock as a collimated jet and started to rise, reaching $8-9$ Mm height with
the mean velocity of $20-30$ km s-1. After a few hundreds of seconds, the
dense and cold plasma of the jet started to flow down due to the gravitational
attraction. On the other hand, the pulses generated the quasi-periodic wake
with secondary shocks. These shocks formed quasi-periodic upward motions of
chromospheric partially ionized plasma in the form of collimated jets.
However, for the pulses were launched at $Y=2.2$ Mm height, no jets were
observed as the excited waves did not steepen into shocks, and hence the
chromospheric plasma did not rise into the corona. The important contribution
of the magnetic structure is that the null points lead to the particular
inverted-Y shape of the jets near the footpoints. Therefore, the observed
inverted-Y shape of jets is probably connected to the reconnection events in
the chromosphere rather than the corona. This point needs further study by
observations and simulations.

Our simulations show that, according to our expectations, ions and neutrals
evolve similarly from the high ion-neutral collision frequency. This fact
means collisions couple the charged and neutral fluids on longer time scales
than the collision time scales. On the other hand, the temperature around the
jets significantly grows, which could be connected to the plasma heating due
to the collisions. This temperature increase may explain the observed rapid
disappearance of some jets from cold chromospheric spectral lines and their
consequent appearance in hot transition region lines.

Another interesting consequence of the simulations is that the period of
quasi-periodic upward motions of jets is longer than the acoustic cut-off
period, which for the C7 model, ranges from about $250$ s at the photosphere
($y=0$ Mm) to $350$ s at the upper chromosphere ($y=2$ Mm). In particular, for
near $y=y_{0}=1.6$ Mm, the cut-off period is approximately $350$ s, but the
periodicity of plasma oscillations over the TR is almost $450$ s, as shown on
the left panel of Fig. 9. This is associated with the fact that the
periodicity of quasi-periodic shocks depends on the initial amplitude of the
pulse: stronger initial pulses lead to the longer periods of secondary shocks
in the stratified atmosphere (Zaqarashvili et al., 2011b).

In conclusion, the pressure pulses excited (presumably by magnetic
reconnection) near the magnetic null points in the chromosphere lead to the
formation of quasi-periodic collimated jets, which rise into the hot solar
corona. The lower part of these jets has an inverted-Y shape, which is a
typical feature of some observed jets. In addition, the upper parts of these
jets are heated, which could be connected to ion-neutral collisions. Future
simulations, including variable ionization/recombination, may reveal more
detailed features of the jet heating due to ion-neutral effects in partially
ionized plasma.

## Acknowledgements

We thank the referee for the constructive comments that helped us to improve
the manuscript. JJGA thanks to Investigadores por México-CONACYT (CONACYT
Fellow), CONACYT LN 315829 (2021) and CONACYT-AEM 2017-01-292684 grants for
partially support this work. The program "investigadores for México", project
1045 sponsor space Weather Service Mexico (SCIESMEX). We carry out the
simulations in the facilities of “Centro de Supercćomputo de Clima Espacial
(CESCOM)" part of the “Laboratorio Nacional de Clima Espacial (LANCE)" at
IGUM, Morelia, Michoacan, Mexico. The JOANNA code was developed by Darek
Wójcik with some contribution from Piotr Wołoszkiewicz and Luis Kadowaki. KM’s
work was done within the framework of the National Science Centre (NCN) grant
projects nos. 2017/25/B/ST9/00506 and 2020/37/B/ST9/00184. TVZ was supported
by the Austrian Fonds zur Förderung der Wissenschaftlichen Forschung (FWF)
project P30695-N27 and and by Shota Rustaveli National Science Foundation of
Georgia (project FR-21-467). We visualize the simulation data using the VisIt
software package (Childs et al., 2012).

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Al Shidi et al. (2019) Al Shidi Q., Cohen O., Song P., Tu J., 2019, Sol. Phys., 294, 124
  * Alharbi et al. (2021) Alharbi A., Ballai I., Fedun V., Verth G., 2021, Monthly Notices of the Royal Astronomical Society, 501, 1940
  * Avrett & Loeser (2008) Avrett E. H., Loeser R., 2008, ApJS, 175, 229
  * Bagwell & Ma (2020) Bagwell T. E., Ma Z. W., 2020, The Astrophysical Journal, 898, 167
  * Ballester et al. (2018) Ballester J. L., et al., 2018, Space Sci. Rev., 214, 58
  * Braginskii (1965) Braginskii S. I., 1965, Reviews of Plasma Physics, 1, 205
  * Chapman & Cowling (1970) Chapman S., Cowling T. G., 1970, The mathematical theory of non-uniform gases. an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. Cambridge University Press
  * Childs et al. (2012) Childs H., et al., 2012, VisIt: An End-User Tool For Visualizing and Analyzing Very Large Data. Chapman and Hall/CRC
  * Choe & Lee (1996) Choe G. S., Lee L. C., 1996, The Astrophysical Journal, 472, 360
  * Close et al. (2004) Close R. M., Parnell C. E., Longcope D. W., Priest E. R., 2004, The Astrophysical Journal, 612, L81
  * De Pontieu et al. (2004) De Pontieu B., Erdélyi R., James S. P., 2004, Nature, 430, 536
  * Dedner et al. (2002) Dedner A., Kemm F., Kröner D., Munz C. D., Schnitzer T., Wesenberg M., 2002, Journal of Computational Physics, 175, 645
  * Deubner & Gough (1984) Deubner F.-L., Gough D., 1984, Annual Review of Astronomy and Astrophysics, 22, 593
  * Draine et al. (1983) Draine B. T., Roberge W. G., Dalgarno A., 1983, ApJ, 264, 485
  * Durran (2010) Durran D. R., 2010, Numerical Methods for Fluid Dynamics. Springer, New York, NY
  * Freed et al. (2015) Freed M. S., Longcope D. W., McKenzie D. E., 2015, Sol. Phys., 290, 467
  * Galsgaard & Nordlund (1997) Galsgaard K., Nordlund Å., 1997, J. Geophys. Res., 102, 231
  * Hollweg (1982) Hollweg J. V., 1982, ApJ, 257, 345
  * Judge (2006) Judge P., 2006, in Leibacher J., Stein R. F., Uitenbroek H., eds, Astronomical Society of the Pacific Conference Series Vol. 354, Solar MHD Theory and Observations: A High Spatial Resolution Perspective. p. 259
  * Khomenko et al. (2014) Khomenko E., Collados M., Díaz A., Vitas N., 2014, Physics of Plasmas, 21, 092901
  * Kiss et al. (2017) Kiss T. S., Gyenge N., Erdélyi R., 2017, The Astrophysical Journal, 835, 47
  * Kuridze et al. (2009) Kuridze D., Zaqarashvili T. V., Shergelashvili B. M., Poedts S., 2009, A&A, 505, 763
  * Kuźma et al. (2015) Kuźma B., Murawski K., Solov´ev A., 2015, A&A, 577, A138
  * Kuźma et al. (2017) Kuźma B., Murawski K., Kayshap P., Wójcik D., Srivastava A. K., Dwivedi B. N., 2017, The Astrophysical Journal, 849, 78
  * Kuźma et al. (2021) Kuźma B., Murawski K., Poedts S., 2021, Monthly Notices of the Royal Astronomical Society, 506, 989
  * Leake & Linton (2013) Leake J. E., Linton M. G., 2013, The Astrophysical Journal, 764, 54
  * Leake et al. (2014) Leake J. E., et al., 2014, Space Sci. Rev., 184, 107
  * Lipartito et al. (2014) Lipartito I., Judge P. G., Reardon K., Cauzzi G., 2014, The Astrophysical Journal, 785, 109
  * Longcope (2005) Longcope D. W., 2005, Living Reviews in Solar Physics, 2, 7
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, Journal of Computational Physics, 208, 315
  * Morita et al. (2010) Morita S., Shibata K., UeNo S., Ichimoto K., Kitai R., Otsuji K.-i., 2010, Publications of the Astronomical Society of Japan, 62, 901
  * Murawski & Zaqarashvili (2010) Murawski K., Zaqarashvili T. V., 2010, A&A, 519, A8
  * Murawski et al. (2011a) Murawski K., Zaqarashvili T. V., Nakariakov V. M., 2011a, A&A, 533, A18
  * Murawski et al. (2011b) Murawski K., Srivastava A. K., Zaqarashvili T. V., 2011b, A&A, 535, A58
  * Murawski et al. (2020) Murawski K., Musielak Z. E., Wójcik D., 2020, The Astrophysical Journal, 896, L1
  * Nelson et al. (2019) Nelson C. J., Freij N., Bennett S., Erdélyi R., Mathioudakis M., 2019, The Astrophysical Journal, 883, 115
  * Nishizuka et al. (2008) Nishizuka N., Shimizu M., Nakamura T., Otsuji K., Okamoto T. J., Katsukawa Y., Shibata K., 2008, The Astrophysical Journal, 683, L83
  * Nishizuka et al. (2011) Nishizuka N., Nakamura T., Kawate T., Singh K. A. P., Shibata K., 2011, The Astrophysical Journal, 731, 43
  * Oliver et al. (2016) Oliver R., Soler R., Terradas J., Zaqarashvili T. V., 2016, The Astrophysical Journal, 818, 128
  * Pennicott & Cally (2021) Pennicott J. D., Cally P. S., 2021, Sol. Phys., 296, 26
  * Pylaev et al. (2017) Pylaev O. S., Zaqarashvili T. V., Brazhenko A. I., Melnik V. N., Hanslmeier A., Panchenko M., 2017, A&A, 601, A42
  * Régnier et al. (2008) Régnier S., Parnell C. E., Haynes A. L., 2008, A&A, 484, L47
  * Sabri et al. (2020) Sabri S., Vasheghani Farahani S., Ebadi H., Poedts S., 2020, Sci Rep, 10, 14
  * Secchi (1877) Secchi S. J., 1877, Le Soleil. Paris: Gauthier-Villars
  * Shibata et al. (2007) Shibata K., et al., 2007, Science, 318, 1591
  * Singh et al. (2012) Singh K. A. P., Isobe H., Nishizuka N., Nishida K., Shibata K., 2012, The Astrophysical Journal, 759, 33
  * Smirnova et al. (2016) Smirnova V., Konkol P. M., Solov’ev A. A., Murawski K., 2016, Sol. Phys., 291, 3207
  * Smith & Sakai (2008) Smith P. D., Sakai J. I., 2008, A&A, 486, 569
  * Soler et al. (2013) Soler R., Carbonell M., Ballester J. L., 2013, The Astrophysical Journal Supplement Series, 209, 16
  * Srivastava et al. (2017) Srivastava A. K., et al., 2017, Scientific Reports, 7, 43147
  * Srivastava et al. (2018) Srivastava A. K., et al., 2018, Nature Astronomy, 2, 951
  * Sterling (2000) Sterling A. C., 2000, Sol. Phys., 196, 79
  * Sterling & Hollweg (1988) Sterling A. C., Hollweg J. V., 1988, ApJ, 327, 950
  * Sterling et al. (2015) Sterling A. C., Moore R. L., Falconer D. A., Adams M., 2015, Nature, 523, 437
  * Tian et al. (2018) Tian H., et al., 2018, The Astrophysical Journal, 854, 92
  * Vögler et al. (2005) Vögler A., Shelyag S., Schüssler M., Cattaneo F., Emonet T., Linde T., 2005, A&A, 429, 335
  * Vranjes & Krstic (2013) Vranjes J., Krstic P. S., 2013, A&A, 554, A22
  * Wójcik et al. (2019) Wójcik D., Murawski K., Musielak Z. E., 2019, The Astrophysical Journal, 882, 32
  * Wójcik (2017) Wójcik D., 2017, MSc Thesis. UMCS, Lublin
  * Wójcik et al. (2018) Wójcik D., Murawski K., Musielak Z. E., 2018, Monthly Notices of the Royal Astronomical Society, 481, 262
  * Yang et al. (2013) Yang L., He J., Peter H., Tu C., Zhang L., Feng X., Zhang S., 2013, The Astrophysical Journal, 777, 16
  * Yeates & Hornig (2011) Yeates A. R., Hornig G., 2011, Physics of Plasmas, 18, 102118
  * Zaqarashvili et al. (2011a) Zaqarashvili T. V., Khodachenko M. L., Rucker H. O., 2011a, A&A, 529, A82
  * Zaqarashvili et al. (2011b) Zaqarashvili T. V., Murawski K., Khodachenko M. L., Lee D., 2011b, A&A, 529, A85