# Generation and Life Cycle of Solar Spicules

Hamid Saleem  Theoretical Research Institute, Pakistan Academy of Sciences,
3-Constitution Avenue, G-5/2, Islamabad 44000, Pakistan.  Department of
Physics, School of Natural Sciences (SNS), National University of Sciences and
Technology (NUST), H-12, Islamabad 44000, Pakistan  Space and Astrophysical
Research Laboratory (SARL), National Center of GIS and Space Applications,
Islamabad 44000, Pakistan  Zain H. Saleem  Argonne National Laboratory, 9700
S. Cass Ave, Lemont, IL 60439, USA.

###### Abstract

Physical mechanism for the creation of solar spicules is proposed with three
stages of their life cycle. It is assumed that at stage-I, the density hump is
formed locally in the xy-plane in lower chromosphere in the presence of
temperature gradients of electrons and ions along z-axis (the vertical
direction). In this region, the density structure of quasi-neutral
$(n_{i}\simeq n_{e}=n)$ plasma after taking birth is accelerated in the
vertical direction due to the thermal force ${\bf F}_{th}\propto\nabla
n(x,y,t)\times(\nabla T_{e}+\nabla T_{i})$. The exact time-dependent
analytical solution of two fluid plasma equations is presented assuming that
density is maximum at the center and decays away from it gradually. The two
dimensional (2D) density structure is created as a step function $H(t)$ in
time at bottom of the chromosphere and consequently the vertical plasma
velocity turns out to be the ramp function of time $R(t)=tH(t)$ whereas the
source term $S(x,y,t)$ for the density follows the delta function $\delta(t)$
form. The upward acceleration ${\bf a}=a(x,y)\hat{z}$ produced in this density
structure is greater than the downward constant solar acceleration $-{\bf
g}_{\odot}$ in the chromosphere. In the transition region (TR), the
temperature gradients are steeper; therefore, the upward acceleration
increases in magnitude $g_{\odot}<<a$ and the density hump spends lesser time
there. This is stage-II of its life cycle. In stage-III, the density structure
enters into corona where the gradients of temperatures vanish and the
structure decelerates to zero velocity under the action of the solar
gravitational force.

Generation of Solar Spicules, Life Cycle of Solar Spicules, Astrophysical
Jets, Two Fluid Plasma, Time-dependent Baroclinic Vectors

##  1 Introduction

Material ejection in the vertical direction from the surfaces of astronomical
objects is a very common phenomenon. Small-scale plasma jets, called spicules,
were first observed long ago Secchi (1877) in the solar atmosphere. The
physical mechanism for the generation of large-scale classical plasma jets has
similarities with the small-scale outflows of plasma observed in different
regions of the solar atmosphere Shibata & Uchida (1986); Bellan (2018).
Vertical plasma jets emerging in the perpendicular direction to the solar disk
can always be observed, which have lengths of the order of thousands of
kilometers, diameters of the order of hundreds of kilometers and life time of
the order of (5-10) minutes. These density structures move in the upward
direction with speeds of tens of kilometers per second Priest (1982). About
half a million Lang (2001) or at least about hundred thousand $100,000$ De
Pontieu et al. (2004) spicules can always be seen at the solar disk with
variations in size. So far, there has been no consensus on the fundamental
mechanism for the generation of astrophysical jets at different spatial
scales. Several physical mechanisms can be involved in producing these jet-
like flows in different astronomical environments.

Sun is the closest astronomical object to our home, Earth, and a great variety
of small-scale localized jet-like plasma flows and mass ejections have been
detected in the three regions of its atmosphere; the surface, the chromosphere
and the corona Priest (1982, 2014); Woo (1996); De Pontieu et al. (2007a);
Aschwanden & Peter (2017); Klimchuk (2006). It includes spicules, coronal
loops, coronal mass ejections (CMEs), and the much smaller cylindrical plasma
structures; the threads and strands observed within the coronal loops Goddard
et al. (2017).

During the past several years, a lot of interest has been invoked in the study
of solar spicules, particularly because of the advancement in observational
techniques which led De Pontieu et al. (2007a) to divide the spicules into two
categories; one has been named as Type-I and the other Type-II. Type-I
spicules exhibit slower velocities $(10-40)kms^{-1}$ and longer life times
$(3-10)$ minutes, while Type-II spicules have larger velocities
$(80-300)kms^{-1}$ and shorter life times $(1-3)$ minutes De Pontieu et al.
(2007a). In addition to upward plasma motions, downward flows have also been
observed in the transition region (TR) and lower chromosphere. High-speed
downward flows with velocities $(60-200)kms^{-1}$ occur near the active
regions as well as in the quiet Sun (Bose et al. 2021). Thermal energy flows
from the chromosphere to the corona, but only a small fraction of energy
escapes to the corona Carlsson et al. (2019) and the rest remains trapped
within the chromosphere. After their life time, spicules either fade out or
fall back into the chromosphere Bose et al. (2021); Scalisi et al. (2021a).

In order to gain a deeper understanding of the morphology and characteristics
of solar jets, a few physical mechanisms have been proposed. These mechanisms
include the shock waves generated by pressure/velocity pulses Hollweg (1982);
Shibata et al. (1982); Kuźma et al. (2017a), shock waves driven by nonlinear
Alfven waves in magnetic flux tubes Matsumoto & Shibata (2010) as well as the
creation of jets by the magnetic reconnection Yokoyama & Shibata (1995);
González-Avilés et al. (2017). The reflection of Alfven waves has also been
considered as one of the possible mechanisms for the creation of spicules
Scalisi et al. (2021b).

Numerical simulations can consider a multitude of physical effects by
including the terms representing radiation, conduction, resistivity, and
pressure in the set of plasma equations to investigate the evolution of jet-
like flows. Series of numerical experiments have been performed using single
fluid magnetohydrodynamics (MHD). The role of several parameters in
determining the dynamics of spicules has been studied Martínez-Sykora et al.
(2017, 2020). The 2D simulations provide a clue to the existence of pressure
force within the spicule structure, and the plasma within the boundaries is
found to be in the state of non-equilibrium. The non-uniform internal pressure
force and magnetic field tension counterbalance each other and cause cross-
sectional deformation. Regions of high density and pressure have been given
the names of knots by some authors Dover et al. (2021). A new phenomenon of
nanojets has also been observed recently Antolin et al. (2021).

The creation and dynamical evolution of the spicules have also been
investigated through numerical simulations by launching the vertical velocity
from the upper chromosphere using MHD equations Kuźma et al. (2017a). The set
of equations used in this simulation reduces to MHD if neutrals are ignored.
The authors have considered several physical effects, such as gravitational
field, ionization, recombination, collisions, and pressure gradient forces.
The current-free plasma is considered and in equilibrium the pressure gradient
balances the gravity. Then a time-dependent perturbation is introduced at the
bottom of the plasma boundary to investigate the temporal evolution of
densities and pressures. The results indicate a rapid increase in the plasma
density with height which is not in complete agreement with the observations.
However, several other features indicated by the simulations match with the
observations. The MHD simulations Dover et al. (2021); Kuźma et al. (2017a)
predict the role of pressure gradients in the formation and evolution of
spicules, and they do not rely only on the magnetic reconnection mechanism.
Numerical simulation has also been performed for the development of spicules
with velocity pulse generated in the chromosphere using ideal and non-ideal
MHD models including pressure force Kuźma et al. (2017b).

Long ago Biermann (1950), it was proposed that the electron baroclinic vector
$(\nabla n_{e}\times\nabla T_{e})$ generates the seed magnetic fields in
stars. Ions were assumed to be static $(m_{i}\rightarrow\infty)$ and electrons
were considered to be inertialess $(m_{e}\rightarrow 0)$ in this theoretical
model. The spatial variations of density and temperature were assumed to be
constant in time and the electron Lorentz force term was also ignored. The
large magnetic fields observed in classical laser plasma experiments of the
order of megagauss ($10^{6}$) $G$ Brueckner & Jorna (1974) were explained on
the basis of the Biermann battery mechanism. Keeping in view the intermediate
temporal and spatial scales for the case $m_{i}\rightarrow\infty$ and
$m_{e}\rightarrow 0$, a theoretical model named electron magnetohydrodynamics
(EMHD) was presented Kingsep et al. (1990); Bol’shov et al. (1991). This
mechanism was also used to estimate the seed fields in galaxies of the order
of ($3\times 10^{-17}$ $G$) Lazarian (1992); Widrow (2002). It may be
mentioned here that the baroclinic flow, baroclinic waves and instabilities
are well-known in atmospheric physics. The equations used in different models
are solved under appropriate boundary conditions applicable in the atmosphere
De-Leon et al. (2020); Baran et al. (2022); Agaoglou et al. (2024).

It was shown that the assumption of stationary ions in Biermann mechanism and
in EMHD model was not convincing for the study of the problem of magnetic
field generation Saleem (1996, 1999). Therefore, the Biermann mechanism was
modified by including ion dynamics and plasma flows to generate the coupled
seed magnetic field ${\bf B}$ and plasma vorticity $(\nabla\times{\bf
v}_{i})$. The 2D exact solutions of two fluid plasma equations including flows
of electrons and ions were found assuming spatial variation of density and
temperatures to be constant with respect to time in Cartesian geometry Saleem
(2007, 2010). In these investigations, it was pointed out that the plasma
vorticity is coupled with the magnetic field and both physical quantities are
created simultaneously by the baroclinic vectors of electrons and ions. Hence
the mechanism of magnetic field generation should not be considered
independent of vorticity (or flow) generation. Later on Saleem (2021), the
above theoretical model was extended to explore the possibility of generation
of jet-like flows by thermal forces of electrons and ions using cylindrical
coordinates. It was shown that plasma flow in the perpendicular direction to
the surface is generated if the density attains locally the spatial profile
like the Bessel function of order one in the radial direction and the
temperatures vary linearly along the axial direction In this theoretical
model, the jet-like flows are produced in both directions; upward and downward
because the Bessel function is an oscillatory function.

Recently, the 3D exact solutions of two fluid plasma, MHD and neutral fluid
equations have been presented in Cartesian geometry Saleem & Saleem (2022)
assuming the gradients of density and temperatures to be constant with respect
to time. Following the 2D cylindrical solution of plasma equations for the
creation of jet-like flows mentioned above, a model for the creation of plasma
jets in Cartesian coordinates has also been presented. In that work, it has
been shown that plasma flow in the upward direction is produced when the
density has a spatial dependence on the surface coordinates $n=n(x,y)$ in a
local region and the electron and ion temperatures have positive gradients
along the vertical direction, that is, $\nabla
T_{j}=\mid\frac{dT_{j}}{dz}\mid\hat{z}$ where $j=(e,i)$. This exact analytical
solution of plasma equations has been applied to explain the generation of
spicules by plasma baroclinic vectors $\nabla n\times\nabla T_{j}$. The
creation of the plasma slab at the bottom of the spicule can be divided into
four quadrants in the xy plane. The density is assumed to be maximum at the
center. Considering only one out of four quadrants, it has been shown that a
slab of a small height at the bottom of the structure will be lifted upward
with the acceleration greater than the solar gravitational acceleration acting
downward. When the slab enters into upper regions, either the density gradient
or temperature gradient vanishes and its speed becomes constant in agreement
with the observations. The vorticity and magnetic field turn out to be
linearly growing with time by the time-independent baroclinic vectors of
electrons and ions. But this model cannot consider the time evolution of the
plasma jet at the initial stage.

Physically a density hump is created in space during a finite time ’t’ however
small it may be. If constant temperature gradients are given as is
approximately true in different regions of the solar atmosphere Slemzin et al.
(2014); Priest (2014), then wherever a density hump or dip is created locally
as a function of $(x,y)$, it must be time-dependent at least at initial
evolution stage of the ordered plasma structure; the plasma jet. A more
realistic solution of plasma equations must have a time-dependent density
$n=n(x,y,t)$ to produce upward flows locally in the presence of constant
temperature gradients along $z$ axis.

Our task is two fold; first we find the exact analytical solution of two fluid
plasma equations with time-dependent density $n=n(x,y,t)$ and constant
gradients of temperatures along the positive z-axis $T_{j}=T_{j}(z)$. This
solution can also be used to investigate the generation of large-scale plasma
jets emerging from YSOs and other astronomical objects in the classical limit.
Second, we apply the above mentioned 2D solution of two fluid plasma equations
to explain the evolution of spicules and their life cycle. The birth and death
of spicules are explained by considering three stages in three different
regions; chromosphere, transition region (TR), and lower corona. The 2D time-
dependent density structure presented in this work will be used in the
evolution stage of the spicule in the chromosphere. Then the 2D solution
Saleem & Saleem (2022) obtained for time-independent density will be employed
to elaborate the dynamics of the spicule in TR and corona using basic
classical kinematics.

##  2 Theoretical Model for Plasma Jet Formation

Let us consider two fluid ideal plasma equations of momentum conservation for
electrons and ions in the presence of constant gravitational acceleration
${\bf g}$ which are written, respectively, as Chen et al. (1984),

$m_{e}n_{e}(\partial_{t}+{\bf v}_{e}\cdot\nabla){\bf v}_{e}=-en_{e}({\bf
E}+\frac{1}{c}{\bf v}_{e}\times{\bf B})-\nabla p_{e}+m_{e}n_{e}{\bf g}$ (1)

and

$m_{i}n_{i}(\partial_{t}+{\bf v}_{i}\cdot\nabla){\bf v}_{i}=en_{i}({\bf
E}+\frac{1}{c}{\bf v}_{i}\times{\bf B})-\nabla p_{i}+m_{i}n_{i}{\bf g}$ (2)

Since spicules have an ambient magnetic field created by solar plasma
dynamics, therefore we add a unidirectional constant magnetic field ${\bf
B}_{0}=B_{0}\hat{z}$ to the weak magnetic field generated by the baroclinic
vectors ${\bf B}_{(}x,y,t)=\tilde{\bf B}$ and write the total field within the
structures ${\bf B}_{T}$ as,

${\bf B}_{T}={\bf B}(x,y,t)+{\bf B}_{0}$ (3)

Both electrons and ions are assumed to obey the ideal gas law
i.e.$p_{j}=n_{j}T_{j}$ where subscript $j=(e,i)$ denotes electrons and ions.
It will be shown that when a density hump is created within a local region of
the chromosphere in xy-plane in the presence of constant gradients of the
temperatures along ${\bf B}_{0}$, it will be accelerated upward by the force
of baroclinic vectors $\nabla n\times\nabla T_{j}$. The quasi-neutral plasma
$n_{e}\simeq n_{i}=n(x,y,t)$ is in a state of non-equilibrium $T_{e}\neq
T_{i}$. The continuity equations require a source term $S(x,y,t)$ and they can
be written for jth species in the following form,

$\partial_{t}n+\nabla\cdot(n{\bf v}_{j})=S(x,y,t)$ (4)

Since the current density is ${\bf j}=en({\bf v}_{i}-{\bf v}_{e})$, therefore
Amperes’ law gives electron velocity in terms of ions velocity,

${\bf v}_{e}={\bf v}_{i}-\frac{c}{4\pi e}(\frac{\nabla\times{\bf B}}{n})$ (5)

The equations of motion of electrons and ions for longitudinally uniform flow
$\nabla\cdot{\bf v}_{j}=0$ take the form of following two coupled equations
Saleem (2010), respectively,

$\partial_{t}{\bf B}=\nabla\times({\bf v}_{i}\times{\bf B})+\nabla\times({\bf
v}_{i}\times{\bf B}_{0})-(\frac{c}{4\pi ne})\\{\nabla\times[(\nabla\times{\bf
B})\times{\bf B}]\\}$ (6) $-(\frac{c}{4\pi
ne})\\{\nabla\times[(\nabla\times{\bf B})\times{\bf B}_{0}]\\}+\frac{c}{4\pi
ne}\\{\nabla\psi\times[(\nabla\times{\bf B})\times{\bf B}]\\}$ $+\frac{c}{4\pi
ne}\\{\nabla\psi\times[(\nabla\times{\bf B})\times{\bf
B}_{0}]\\}-\frac{c}{e}(\nabla\psi\times\nabla T_{e})$

and

$a\partial_{t}{\bf B}+\partial_{t}(\nabla\times{\bf v}_{i})=\nabla\times[{\bf
v}_{i}\times(\nabla\times{\bf v}_{i})]+a\nabla\times({\bf v}_{i}\times{\bf
B})$ (7) $+a\nabla\times({\bf v}_{i}\times{\bf
B}_{0})+\frac{1}{m_{i}}(\nabla\psi\times\nabla T_{i})$

where $a=\frac{e}{m_{i}c}$, $\psi=\ln\bar{n}$, $\bar{n}=\frac{n}{N_{0}}$, and
$N_{0}$ is arbitrary constant density. If following two conditions also hold,

$\nabla\psi\cdot{\bf v}_{i}=0$ (8)

and

$\nabla\psi\cdot(\nabla\times{\bf B})=0$ (9)

then all nonlinear and complicated terms of Eqs. (6) and (7) vanish and they
reduce to two simpler equations Saleem (2010, 2021); Saleem & Saleem (2022),
respectively, given as,

$\partial_{t}{\bf B}=-\frac{c}{e}(\nabla\psi\times\nabla T_{e})$ (10)

and

$\frac{e}{m_{i}c}\partial_{t}{\bf B}+\partial_{t}(\nabla\times{\bf
v}_{i})=\frac{1}{m_{i}}(\nabla\psi\times\nabla T_{i})$ (11)

We assume $\nabla\psi$ is not parallel to $\nabla T_{j}$. Then (10) and (11)
give an expression for the generation of ions vorticity by baroclinic vectors,

$\partial_{t}(\nabla\times{\bf v}_{i})=\frac{1}{m_{i}}\nabla\psi\times(\nabla
T_{e}+\nabla T_{i})$ (12)

Equation (12) implies that the thermal force of plasma is the source for the
generation of ions vorticity which is coupled with the electrons vorticity
through Eq. (5).

###  2.1 Generation of Plasma Jet

If $\psi=\psi(x,y,t)$ and $T_{j}=T_{j}(z)$, then the right hand side of Eq.
(12) can be expressed as,

$\frac{1}{m_{i}}\nabla\psi\times(\nabla T_{e}+\nabla
T_{i})=a_{0}(\partial_{y}\psi,-\partial_{x}\psi,0)$ (13)

where $a_{0}=(\frac{T_{e}^{\prime}+T_{i}^{\prime}}{m_{i}})$ and
$T_{j}^{\prime}=\frac{dT_{j}}{dz}$. We have $\nabla\times{\bf
v}_{i}=[(\partial_{y}v_{iz}-\partial_{z}v_{iy}),(\partial_{z}v_{ix}-\partial_{x}v_{iz}),(\partial_{x}v_{iy}-\partial_{y}v_{ix})]$
and relation (13) demands that $\partial_{x}v_{iy}-\partial_{y}v_{ix}=0$, in
general, which is possible if $v_{ix}=v_{iy}=0$. In this case, Eq. (12)
becomes,

$\partial_{t}(\partial_{y}v_{iz},-\partial_{x}v_{iz},0)=a_{0}(\partial_{y}\psi,-\partial_{x}\psi,0)$
(14)

The above relation shows that the ions velocity has non-zero component only in
the direction perpendicular to the plasma surface that is along z-axis. The
temperatures in solar atmosphere vary along z-axis almost linearly in
different layers Slemzin et al. (2014); Priest (2014), therefore we assume,

$T_{j}=T_{j}^{\prime}z+T_{0j}$ (15)

where $T_{j}^{\prime}$ and $T_{0j}$ are constants. Density is assumed to be
the spatial function of the coordinates $(x,y)$ and time ”t”, therefore
$\psi=\psi(x,y,t)$. The ions velocity is determined by the Eqs. (12) and (13)
in the following form,

${\bf v}_{i}=(0,0,v_{iz})=v_{i}(x,y,t)\hat{z}$ (16)

Equations (10) and (11) imply that the baroclinic vectors of ions and
electrons are parallel to each other and Eq. (14) indicates that the spatial
profile of $v_{i}$ is similar to $\psi$, therefore, ${\bf B}$ is also parallel
to the ion vorticity $\nabla\times{\bf v}_{i}$. Hence, we assume,

${\bf B}=b_{0}(\nabla\times{\bf v}_{i})$ (17)

where $b_{0}$ is constant. The form of $v_{i}$ given by Eq. (16) yields,

$\nabla\times[{\bf v}_{i}\times(\nabla\times{\bf v}_{i})]=0$ (18)

Equations (17) and (18) imply,

$\nabla\times({\bf v}_{i}\times{\bf B})=0$ (19)

Equation (17) yields,

$\nabla\times{\bf B}=-b_{0}\nabla^{2}{\bf v}_{i}$ (20)

If

$\nabla^{2}{\bf v}_{i}=\eta{\bf v}_{i}$ (21)

where $\eta$ is constant, then both generated fields ${\bf B}$ and ${\bf
v}_{i}$ turn out to be curls of each other. Equations (14) and (21) demand,

$\nabla^{2}\psi=\eta\psi$ (22)

##  3 Analytical Solutions

The general solutions of Eqs. (10), (11) and (12) depend upon the form of
$\psi$ in the presence of spatial gradients of temperatures along the z
direction defined by Eq. (15). We look for an exact analytical solution of two
fluid equations with a time-dependent density function $\psi=\psi(x,y,t)$.

###  3.1 The $\psi$ as a Step Function in Time

Let us consider the physical phenomenon of the the creation of a density hump
in chromosphere at initial time $t=0$ with two dimensional spatial form
$F(x,y)$,

$F(x,y)=(A_{1}e^{-(\mu x+\nu y)}+A_{2}e^{(\mu x-\nu y)})$ (23)

where $A_{1}$, $\mu$, $\nu$, $A_{2}$ are constants. Let the source term in
Eqs. (4) be,

$S(x,y,t)=nF(x,y)\delta(t)$ (24)

Then continuity equations (4) reduce to,

$\partial_{t}\psi=F(x,y)\delta(t)$ (25)

The above relation defines $\psi$ as follows,

$\psi(x,y,t)=F(x,y)H(t)$ (26)

where $H(t)$ is the step function and hence $H(t)=1$ for $0\leq t$ which
gives,

$\psi(x,y,0)=F(x,y)$ (27)

The Eq. (23) implies,

$\nabla^{2}F=\eta F$ (28)

where $\eta=(\mu^{2}+\nu^{2})$ and it is in agreement with Eq. (22). Equation
(14) can be written in terms of $F$ and $H(t)$ using the definition of
function $\psi$,

$\partial_{t}(\partial_{y}v_{i},-\partial_{x}v_{i},0)=a_{0}(\partial_{y}\psi,-\partial_{x}\psi,0)=a_{0}(\partial_{y}F,-\partial_{x}F,0)H(t)$
(29)

The above equation implies that $v_{i}$ should have the spatial profile
similar to $F$.

$\partial_{t}v_{i}(x,y,t)=a_{0}F(x,y)H(t)=a(x,y)H(t)$ (30)

Integration of above equation with the initial condition $v_{i}=0$ at $t=0$,
yields the ions vertical flow as,

$v_{i}(x,y,t)=a_{0}F(x,y)R(t)=a(x,y)t$ (31)

where $R(t)=tH(t)$ is the ramp function which obeys the following relation,

$\frac{d}{dt}R(t)=H(t)$ (32)

Equations (16) and (31) determine the acceleration in vertical direction,

${\bf a}(x,y)=a(x,y)\hat{z}=a_{0}F(x,y)\hat{z}$ (33)

which is also a function of $(x,y)$ coordinates. The density profile will move
upward if the condition $g_{\odot}<a_{0}F(x,y)$ holds. Using Eq. (10) along
with Eq. (26), the magnetic field generated by baroclinic vectors can be
expressed in terms of functions $F$ and $H(t)$ as follows,

$\partial_{t}{\bf
B}(x,y,t)=-\frac{cT_{e}^{\prime}}{e}(\partial_{y}\psi,-\partial_{x}\psi,0)=-\frac{cT_{e}^{\prime}}{e}(\partial_{y}F,-\partial_{x}F,0)H(t)$
(34)

Integration of above equation with the condition $B(x,y,0)=0$ at $t=0$ gives
the magnetic field created by baroclinic vectors as a ramp function of time,

${\bf
B}_{(}x,y,t)=-\frac{cT_{e}^{\prime}}{e}(\partial_{y}F,-\partial_{x}F,0)R(t)$
(35)

Figure 1: The function used to plot the above figure is
$H(x+y)e^{\psi_{0}e^{(-\mu x-\nu y)}}+H(-x+y)e^{\psi_{0}e^{(\mu x-\nu
y)}}+H(x-y)e^{\psi_{0}e^{(-\mu x+\nu y)}}+H(-x-y)e^{\psi_{0}e^{(\mu x+\nu
y)}}$ and here $H(x,y)$ is the Heaviside function on 2D space. The above form
of the function allows us to plot the solution in all four quadrants. It
represents the density profile in xy-plane in the simpler case corresponding
to $A_{2}=0$ in Eq.(23).

###  3.2 The $\psi$ as an Exponential Function in Time

It may be mentioned here that we can also get a mathematically consistent
solution by considering,

$\psi(x,y,t)=F(x,y)e^{\gamma t}$ (36)

where $\gamma$ is a constant.

In this case, we can write Eq. (10) as,

$\partial_{t}{\bf
B}(x,y,t)=-\frac{cT_{e}^{\prime}}{e}(\partial_{y}\psi,-\partial_{x}\psi,0)=-\frac{cT_{e}^{\prime}}{e}(\partial_{y}F,-\partial_{x}F,0)e^{\gamma
t}$ (37)

Integrating above equation with the initial condition $B(x,y,0)=0$ at $t=0$,
we obtain,

${\bf
B}(x,y,t)=-\frac{cT_{e}^{\prime}}{e\gamma}(\partial_{y}F,-\partial_{x}F,0)(e^{\gamma
t}-1)$ (38)

The ions vorticity equation (12) becomes,

$\partial_{t}(\nabla\times{\bf
v}_{i})=\partial_{t}(\partial_{y}v_{i},-\partial_{x}v_{i},0)=a_{0}(\partial_{y}F,-\partial_{x}F,0)e^{\gamma
t}$ (39)

Integrating above equation with the initial condition $v_{iz}(x,y,0)=0$, we
obtain the jet-like flow ${\bf v}_{i}$ to be generated in the following form,

${\bf v}_{i}=\frac{a_{0}}{\gamma}F(x,y)(e^{\gamma t}-1)\hat{z}$ (40)

In this case, the density source can be defined as,

$S(x,y,t)==\gamma n(x,y,t)(e^{\gamma t}\psi_{0}F(x,y))=\gamma n\psi$ (41)

The continuity equations reduce to,

$\partial_{t}\psi=\gamma\psi$ (42)

At $t=0$, Eq.(36) becomes Eq. (27). Since density becomes double exponential
function of time, therefore the parameter $\gamma$ must be chosen very
carefully. Note that in this case, the density $n(x,y,t)$ increases very
rapidly with time and this solution can be valid only over an extremely short
duration of time. Thus, the choice of step function form of $\psi$ during
evolution time seems to be physical and we will use it to investigate the
birth and death of spicules in the next section.

##  4 Creation and Life Cycle of Spicules

Let us suppose that a density hump is created in the form of a rectangle
within the chromosphere with area $\sigma=2x_{m}\times 2y_{m}$ where
$-x_{m}\leq x\leq x_{m}$ and $-y_{m}\leq x\leq y_{m}$. Keeping in view the
observed order of the spatial dimensions of the spicules Priest (1982, 2014),
we choose $x_{m}=(200)km$ and $y_{m}=(300)km$ . The density function $F$ is
maximum at the point $(0,0)$ and decreases smoothly within the rectangular
area away from the centre for $A_{2}<<A_{1}$.

Since, we define density as a function of two coordinates $(x,y)$, the density
profile does not follow the exact exponential form contrary to one dimensional
(1-D) case used for the study of electrostatic drift waves in magnetized
plasmas. It may be mentioned here that the equilibrium density gradient
$\nabla n_{0}$ (where $n_{0}$ denotes equilibrium density) in magnetized
plasmas is assumed to be one-dimensional (1-D), in general, to study the
linear and nonlinear waves. For example, the electrostatic drift wave is an
important low frequency wave which gives rise to instabilities, energy
transport and nonlinear structures in plasmas. In most of the research work on
that topic, the equilibrium density variation is assumed to be of the
exponential form i.e. $n(x)=n_{0}e^{-\frac{x}{L_{n}}}$ where
$\mid\frac{1}{n_{0}}\frac{dn_{0}}{dx}\mid=\kappa_{n}=\frac{1}{L_{n}}=$
constant. The equilibrium density decays exponentially away from the point
$x=0$ where $L_{n}$ is the density gradient scale length Weiland (1999). These
waves were first studied theoretically long time ago Kadomtsev (1963) and the
same year they were observed in laboratory D’Angelo & Motley (1963). Later,
the electrostatic turbulence in magnetized plasma was explained by the
nonlinear drift waves Hasegawa & Mima (1978). If ambient magnetic field is
along the z-direction ${\bf B}_{0}=B_{0}\hat{z}$, then this longitudinal wave
propagates predominantly along positive y-axis. However, in the case of
spicules, the density must have at least two-dimensional (2-D) variation.

The spicules originate from different regions of active and quiet Sun and have
been divided into two main types. The type-I spicules have typical life times
of the order of $3$ to $10$ minutes, diameters of $120$ to $700$ $km$, heights
of $3$ to $5$ $Mm$, temperatures $10,000$ to $15,000$ $K$ electron densities
$10^{17}$ $m^{-3}$ and reach the speeds $10$ to $50$ $kms^{-1}$. Type-II
spicules reach speeds of $30$ to $150$ $kms^{-1}$ and have longer heights
compared to Type-I and are shorter in active regions ($2$ to $5$ $Mm$) and
taller in coronal holes $5$ to $10$ $Mm$ Priest (2014); Beckers (1972); De
Pontieu et al. (2007a). On the other hand, the macrospicules are bigger
structures with heights $4$ to $40$ $Mm$, diameters $4$ to $11$ $Mm$ and have
life time of $8$ to $45$ minutes Priest (2014). Multitude of theoretical
models have been presented for understanding the dynamics of spicules and
causes of their creation Hollweg (1982); De Pontieu et al. (2014); Iijima &
Yokoyama (2015); Shimojo et al. (2020). There is no consensus on any single
complete theory for the creation of spicules and complete comparisons between
simulation results and observations have not been presented. Furthermore, the
observations do not give the detailed changes of density of spicules at
different altitudes. The present theoretical model also does not take into
account the density dependence on z-coordinate. We assume that the spicules
are created in chromosphere, pass through transition region and rise to
heights of thousands of kilometers in corona. At every point within the
spicule, the upward acceleration ${\bf a}(x,y)$ and velocity ${\bf
v}_{i}(x,y)$ have different magnitudes. These three stages are discussed using
the exact solutions of two fluid plasma equations.

In the solar atmosphere, the temperature first decreases from the surface
value $T\simeq 6600$ $K$ to $T\simeq(4300)$ $K$ at an altitude of about $500$
$km$ and then gradually increases from $T\simeq(4300)$ $K$ to $(2.5)\times
10^{4}$ $K$ in the upper chromosphere Priest (2014) and after passing through
the thin transition region (TR), the temperature increases rapidly to
$T\simeq(2-3)\times 10^{5}$ $K$ Slemzin et al. (2014). Above transition region
lies the lower corona where the temperature reaches the value of the order of
$T\simeq 10^{6}$ $K$ and its gradient becomes negligibly small.

The plasma in chromosphere is not in thermal equilibrium and several physical
phenomena such as ionization, recombination, radiation emission and collisions
among charged-charged, charged-neutral and neutral-neutral particles are
continuously taking place Vranjes & Krstic (2013); Kuźma et al. (2019);
Martínez-Sykora et al. (2017). The data of the temperature of electrons
$T_{e}$ and protons $T_{i}$ is not available separately within the region of
chromosphere. However, the values of temperature $T$ at different altitudes
are mentioned in literature assuming the plasma as a single fluid. We are
considering only the hydrogen plasma ignoring the contribution of alpha
particles, neutrals and collisions. Theoretical models from mid chromosphere
to the corona and the solar wind predict that the proton temperature is higher
than the electron temperature $T_{e}<T_{i}$ Hansteen et al. (1997). In the
light of these models, we can expect the condition $T_{e}<T_{i}$ also holds in
chromosphere. Thus we choose $T_{i}=2T_{e}$ and interestingly this relation of
temperatures yields reasonably good estimates of the magnitudes of upward
accelerations, flows and the times spent by the density hump in chromosphere,
transition region and corona. It is important to mention that in the present
theoretical model, the acceleration ${\bf a}$, upward plasma flow ${\bf v}$
and times $\tau$ (spent by the spicule in different regions) are functions of
$F(x,y)$ and hence depend upon $(x,y)$ coordinates. The values of these
quantities have been estimated at two points; the centre point $(0,0)$ and a
point near the edge $(x_{m},y_{m})$ of spicule. It may be noted that we get
$\psi=0$ at the points where $n(x,y,t)=N_{0}$ and the phenomenon of flow
generation disappears. The vertical velocities in each region depend upon the
form of $F$ which is a function of $(x,y)$ coordinates, therefore the density
structure has a non-uniform vertical flow.

###  4.1 Life Cycle of Spicule

Keeping in view the above mentioned temperatures in different regions of solar
atmosphere and their thicknesses, we divide the spicule generation process
into three stages. The first stage is its birth in the lower chromosphere
where the density hump in $(x,y)$ coordinates is created due to the plasma
dynamics in the presence of constant vertical temperature gradients. The
thermal force is proportional to $\nabla\psi\times(\nabla T_{e}+\nabla T_{i})$
and it produces the ion vorticity that is associated with their vertical flow.
It is assumed that the density structure formed at $t=0$ is accelerated in
upward direction and spends time $\tau_{ch}$ in the chromosphere before
leaving it after traveling the distance $1500$ $km$. It then crosses the thin
transition region in a short time $\tau_{tr}$ as the acceleration produced in
TR is larger due to the larger temperature gradients whereas thickness of the
layer is small. After gaining a larger velocity in TR, it enters the corona
where $\nabla T_{j}\rightarrow 0$ and there is no upward force acting on the
structure. The solar gravitational acceleration $(-g_{\odot})$ then slowly
brings the velocity of the structure to zero after a time $\tau_{cr}$. The
total life time of the spicule is $\tau_{l}=(\tau_{ch}+\tau_{tr}+\tau_{cr})$.

###  4.2 Creation of Spicule in Chromosphere

Let us estimate the magnitudes of the acceleration and velocity of the density
hump in this area during time $t=0\rightarrow\tau_{ch}$. Four corner points of
this part of the area are $(0,0)$, $(x_{m},0)$, $(0,y_{m})$ and
$(x_{m},y_{m})$ and we choose dimensionless parameters $\mu x_{m}=(0.3)$ and
$\nu y_{m}=(0.4)$ along with $\frac{n(0,0,0)}{N_{0}}=5$. We find
$\mu=(1.5)\times 10^{-6}$ $m^{-1}$ and $\nu=(1.33)\times 10^{-6}$ $m^{-1}$ for
$x_{m}=(200\times 10^{3})$ $m$ and $y_{m}=(300\times 10^{3})$ $m$. The
function $F(x,y)$ is flexible and the constants $A_{1},A_{2},\mu,\nu$ can have
different values in different cases. Furthermore, the ratios of the densities
at points $(0,0)$ and at $(x_{m},y_{m})$ can be varied and hence the present
theoretical model suggests that vertical structures of different sizes can be
created. The above choice of the values of the constants produces values of
the upward acceleration at all four points
$[(0,0),(x_{m},0),(0,y_{m}),(x_{m},y_{m})]$ greater than the solar gravity
i.e. $g_{\odot}<a(x,y)$. Therefore, the density hump starts moving upward and
attains observable velocity. The magnitude of the upward acceleration $a(x,y)$
created by the baroclinic vectors depends also upon the magnitude of the
temperature gradients. The formation of spicules in chromosphere under the
framework of present theoretical model is based on the following
considerations.

I. Let $n_{0}$ be the backgound density in chromosphere. Then $n_{0}\leq
N_{0}$ while we have chosen $\frac{n(0,0)}{N_{0}}=5$, therefore the density
within the spicule can be even larger than five times the background density.
But this ratio ensures that the whole part of the area in first quadrant is
uplifted because $g_{\odot}<a(x,y)$ at all four corner points. The function
$F$ defined by Eq. (23) with $A_{2}=0.1A_{1}$, gives $F(0,0)=(1.1)A_{1}$ while
we assume $\psi(0,0,0)=\ln 5=\simeq(1.61)$ and comparing these two values, we
obtain $A_{1}\simeq(1.46)$ because $\psi(x,y,0)=F(x,y)$.

II. The profile of density structure is proportional to $F(x,y)$ and it does
not change with $z$ coordinate. This assumption is a limitation of the exact
analytical solutions of two fluid equations.

III. The continuity equations have the source term with $\delta(t)$ and it
forces the $\psi$ function to take the form of a step function $H(t)$. The
density variation is not unidirectional and it does not have exact exponential
form as it is assumed for the investigation of electrostatic drift waves. The
upward acceleration ${\bf a}$, ion velocity ${\bf v}_{i}$, and times $\tau$
spent by the density structure in different regions are functions of $(x,y)$
coordinates.

IV. The magnetic field generated by thermal forces turns out to be much
smaller than the constant ambient magnetic field within the spicule i.e.
$\mid{\bf B}(x,y,t)\mid<<B_{0}$.

V. Density structures are created one after the other in the chromosphere
until the final velocity of the first such structure in the corona becomes
zero. At this stage, the formation of the density lump at the bottom of the
spicule stops.

The plasma temperatures $T$ at about $500$ $km$ and $2000$ $km$ altitudes
above the solar surface are given, respectively, as $T_{e}(z_{1})=4.3\times
10^{3}\;$K and $T_{e}(z_{2})=(2.5)\times 10^{4}\;$ $K$ Priest (2014) and we
assume the plasma temperature to be equal to the electron temperature i.e.
$T=T_{e}$. Accordingly, the electron temperature gradient in chromosphere can
be estimated approximately as follows,

$[T_{e}^{\prime}]_{ch}\simeq\frac{T_{e}(z_{2})-T_{e}(z_{1})}{(z_{2}-z_{1})}=(19.04)\times
10^{-26}\;\hbox{J/m}.$ (43)

where $(z_{2}-z_{1})=(1.5)\times 10^{3}$ $km$. Assuming $T_{i}\simeq 2\;T_{e}$
for hydrogen plasma $m_{i}=1.67\times 10^{-27}$ $kg$, we obtain,

$[a_{0}]_{ch}=\frac{T_{i}^{\prime}+T_{0e}^{\prime}}{m_{i}}\simeq(342.10)m/s^{2}$
(44)

Using Eq. (33), the net upward acceleration in chromosphere becomes,

$[a(x,y)]_{ch}=(a_{0})_{ch}F(x,y)-g_{\odot}$ (45)

where $g_{\odot}=(274)$ $m/s^{2}$ Priest (1982) is magnitude of the
gravitational acceleration in the downward direction at the solar surface. At
$t=\tau_{ch}(x,y)$, we get the vertical velocity of the density hump in
chromosphere before entering the TR,

$[v_{i}(x,y)]_{ch}=[a(x,y)]_{ch}\tau_{ch}$ (46)

where the initial velocity at $t=0$ has been assumed to be zero. The vertical
distance to be travelled by the density is $S_{ch}=(z_{2}-z_{1})=1500$ $km$
and hence the following relation can be used to estimate the time spent by the
density hump in chromospeher,

$S_{ch}=(v_{i})\tau_{ch}+\frac{1}{2}a(x,y)\tau_{ch}^{2}$ (47)

We estimate the magnitudes of velocities and the spent times of the structure
only at two points $(0,0)$ and $(x_{m},y_{m})$. From Eq. (45), we find
$[a(0,0)]_{ch}=275.42$ $m/s^{2}$ and $[a(x_{m},y_{m})]_{ch}=19.22$ $m/s^{2}$.
The spent time $[\tau(x,y)]_{ch}$ as well as initial and final velocities at
points $(0,0)$ and $(x_{m},y_{m})$ are mentioned in Table-1.

Figure 2: Schematic diagram of density hump and its velocity in three
different regions of solar atmosphere. The time axis is in the upward
direction.

###  4.3 Velocity Gain in the Transition Region

The value of the temperature at top of transition region, say at point $z_{3}$
is given $T_{e}(z_{3})=(2-3)\times 10^{5}\;$ $K$ (Slemzin 2014) while the
thickness of this layer is about $100\;$ $km$. We take $T_{e}(z_{3})=3\times
10^{5}$ $K$ and estimate,

$T_{e}^{\prime}\simeq\frac{T_{e}(z_{3})-T_{e}(z_{2})}{(z_{3}-z_{2})}=(3.79)\times
10^{-23}\;\hbox{J/m}.$ (48)

where $(z_{3}-z_{2})=100$ $km$. It gives,

$(a_{0})_{tr}=\frac{T_{i}^{\prime}+T_{e}^{\prime}}{m_{i}}=(68)km/s^{2}$ (49)

The net upward acceleration is,

$[a(x,y)]_{tr}=[a_{0}]_{tr}F(x,y)-g_{\odot}\simeq[a_{0}]_{tr}F(x,y)$ (50)

because in this region $g_{\odot}<<[a_{0}]_{tr}F(x,y)$. We obtain
$[a(0,0)]_{tr}\simeq(109.21)$ $km/s^{2}$ and
$[a(x_{m},y_{m})]_{tr}\simeq(58.15)$ $km/s^{2}$. Considering the thickness of
TR to be $(100)$ $km$, the time $\tau_{tr}$ spent in TR by these two points is
estimated using the following relation,

$[a(x,y)]_{tr}\tau_{tr}^{2}+2[v_{i}(x,y)]_{ch}\tau_{tr}-2\times 10^{5}=0$ (51)

The velocity gained in this region by ions is estimated using following
equation,

$[v_{i}(x,y)]_{tr}=[v_{i}(x,y)]_{ch}+[a(x,y)]_{tr}\tau_{tr}(x,y)$ (52)

The estimated values of velocities and spent time at two points $(0,0)$ and
$(x_{m},y_{m})$ are mentioned in Table-1.

###  4.4 Rise and Fall of Spicule in the Corona

In corona, the upward acceleration vanishes because $\nabla T_{j}\rightarrow
0$ and the density structure enters into this region with the velocity gained
in TR i.e. $[v_{i}(x,y)]_{tr}$. The final ion velocity in corona ,say
$[v_{i}]_{cr}$, becomes zero after the time $\tau_{cr}$ in corona which is
longer than $\tau_{ch}$ and $\tau_{tr}$. Following relation is used to
estimate the time spent in corona.

$[v_{i}(x,y)]_{cr}=[v_{i}(x,y)]_{tr}-g_{\odot}\tau_{cr}(x,y)=0.$ (53)

. In above equation $[v_{i}(x,y)]_{cr}$ is the final velocity of ions in
corona. In this region, the downward solar acceleration does not allow the
spicule to move continuously in the vertical direction. In fact, its velocity
decreases slowly here because $g_{\odot}$ is small and finally the velocity
becomes zero. The corresponding heights $\textit{H}_{cr}(x,y)$ attained in
corona can be estimated using the relation
$\textit{H}_{cr}(x,y)=\frac{1}{2}\frac{[v_{i}(x,y)]_{cr}^{2}}{g_{\odot}}$ and
it gives $\textit{H}_{cr}(0,0)=(21.28)\times 10^{3}$ $km$ and
$\textit{H}_{cr}(x_{m},y_{m})=(10.54)\times 10^{3}$ $km$. Total height of the
spicules is $H_{T}=(1500)+(100)+H_{cr}$ $km$. Therefore, we obtain
$H_{T}(0,0)\simeq(2.28)\times 10^{7}$ $m$ and
$H_{T}(x_{m},y_{m})\simeq(1.21)\times 10^{7}$ $km$ Similarly, the life time
can be estimated
$\tau_{l}(0,0)=[\tau]_{ch}(0,0)+\tau]_{tr}(0,0)+\tau]_{cr}(0,0)\simeq(8.32)$
minutes and $\tau_{l}(x_{m},y_{m})\simeq(11.25)$ minutes.

Figure 3: Schematic diagram of how the upward velocity of density hump
increases slowly in stage-I, gets larger speed in stage-II and finally
decreases to zero in corona after time $\tau_{cr}$. The $\tau_{ch}$ is the
time for which the spicule evolves in the chromosphere and $\tau_{tr}$ is the
time the spicule spends in the transition region.

###  4.5 Magnetic Field Generation in Spicule

Magnitude of the generated magnetic field by electron baroclinic term can be
estimated by using the electron equation. Let us evaluate it at two points
chosen above in the chromosphere during $t=\tau_{ch}$. Equation (35) can be
expressed in the following form,

${\bf
B}(x,y,t)=(\frac{cT_{e}^{\prime}}{e}\nu)t[F(x,y),\frac{\mu}{\nu}A_{1}(-e^{-(\mu
x+\nu y)}+(.1)(e^{(\mu x+\nu y)})),0]$ (54)

The generated magnetic field is estimated at two points in the chromosphere as
follows,

${\bf B}(0,0,\tau_{ch})\simeq(1.64)\times 10^{-6}(1.61,-1.47,0)G$ (55)

and

${\bf B}(x_{m},y_{m},\tau_{ch})\simeq(6.24)\times 10^{-6}(0.84,-0.65,0)G$ (56)

In transition region,

${\bf B}(0,0,\tau_{tr})\simeq(2.26)\times 10^{-6}(1.61,-1.47,0)G$ (57)

and

${\bf B}(x_{m},y_{m},\tau_{ch})\simeq(3.71)\times 10^{-6}(0.84,-0.65,0)G$ (58)

The order of magnitude of the created magnetic field is very small compared to
the ambient magnetic field $B_{0}$ which varies in solar atmosphere from a few
Gauss to a few hundred Gauss Priest (1982). Since $\mu$ and $\nu$ are very
small on astrophysical scales, including the solar atmosphere, the magnetic
fields generated by the baroclinic vectors turn out to be very small.

It may be noted that if ions are assumed to be stationary, we get only
electrons Eq. (10) for the magnetic field generation. In this limit, Biermann
battery mechanism is generalized by including electron density dependence on
$(x,y)$ coordinates as well the time ’t’. The ions vorticity equation (12)
disappears and flow generation mechanism is killed. Equation (10) with $\psi$
defined by Eq. (26) modifies the EMHD result in two ways. One is that the
density gradient scale length is not constant and the other is that the
generated magnetic field is not unidirectional and constant in space. In the
original Biermann battery effect, the electric field is assumed to be created
by the electron pressure gradient, viz,

$0=en_{e}{\bf E}+\nabla p_{e}\\\ $

in the background of the stationary ions. The ${\bf v}_{e}\times{\bf B}$ term
in Lorentz force is also neglected. The curl of above equation yields,

$\partial_{t}{\bf B}=-\frac{c}{e}(\frac{\nabla n_{e}\times\nabla
T_{e}}{n_{e}})=-\frac{cT_{e}}{e}(\frac{\nabla n_{e}}{n_{e}}\times\frac{\nabla
T_{e}}{T_{e}})\\\ $

The density gradient scale length $L_{n}=\mid\frac{\nabla n_{e}}{n_{e}}\mid$
and temperature gradient scale length $L_{T}=\mid\frac{\nabla
T_{e}}{T_{e}}\mid$ are assumed to be unidirectional and constant. Thus the
generated magnetic field also turns out to be unidirectional perpendicular to
both the vectors $\nabla n_{e}$ and $\nabla T_{e}$. The simple Biermann
battery mechanism or EMHD cannot be applied to the plasma of solar spicules
where ions motion is confirmed by observations. To find out an exact two
dimensional (2D) analytical solution of the highly nonlinear set of coupled
partial differential equations of two fluid plasma including the flows of both
species ${\bf v}_{j}\neq 0$, we have ignored the collisions. Dissipative
effects in these physical mechanisms can be included, but in that case the
numerical methods are needed to solve the relevant set of equations for
galactic plasma and solar atmosphere. It has been shown that the baroclinic
vectors can generate flows and magnetic fields simultaneously (Saleem2010,
Saleem2021) and ion dynamics should not be ignored.

In the present physical model we have not used the force free condition
$(\nabla\times{\bf B})\times{\bf B}=0$ because Eq. (35) implies that,

$(\nabla\times{\bf B})\times{\bf
B}=-(\frac{cT_{e}^{\prime}}{e}R(t))^{2}\eta(\partial_{x}F,\partial_{y}F,0)\neq
0\\\ $

where

$\nabla\times{\bf
B}=\frac{cT_{e}^{\prime}}{e}R(t)(0,0,\nabla^{2}F)=\frac{cT_{e0}^{\prime}}{e}R(t)(0,0,\eta
F)$ (59)

Figure 4: Spatial profile of the upward velocity of ions $[v_{i}(x,y)]_{ch}$ as function of $(x,y)$ coordinates in the chromosphere.  |  Time spent by density hump $\tau$(s) |  Ion’s initial velocity $v_{i}(km/s)$ |  Ion’s final velocity $v_{f}(km/s)$ |  Temperature gradient $T_{e}^{\prime}$ $(10^{-26}J/m)$
---|---|---|---|---
Chromosphere (ch) |  $\tau_{ch}(0,0)=104$ $\tau_{ch}(x_{m},y_{m})=395$ |  $[v_{i}(0,0)]_{ch}=0$ $[v_{i}(x_{m},y_{m})]_{ch}=0$ |  $[v_{f}(0,0)]_{ch}=29$ $[v_{f}(x_{m},y_{m})]_{ch}=8$ | 19
Transition Region (tr) |  $\tau_{tr}(0,0)=0.72$ $\tau_{tr}(x_{m},y_{m})=1.18$ |  $[v_{i}(0,0)]_{tr}=29$ $[v_{i}(x_{m},y_{m})]_{tr}=8$ |  $[v_{f}(0,0)]_{tr}=108$ $[v_{f}(x_{m},y_{m})]_{tr}=76$ | 3790
Corona (cr) |  $\tau_{cr}(0,0)=395$ $\tau_{cr}(x_{m},y_{m})=279$ |  $[v_{i}(0,0)]_{cr}=108$ $[v_{i}(x_{m},y_{m})]_{cr}=76$ |  $[v_{f}(0,0)]_{cr}=0$ $[v_{f}(x_{m},y_{m})]_{cr}=0$ | 0
Table 1: The time spent by the density hump $\tau$ and the initial(final)
velocities of ions at two points $(0,0)$ and $(x_{m},y_{m})$ are mentioned in
three different regions to show the life cycle of the spicules from birth in
the chromosphere with initial velocity zero to death in the corona with final
velocity zero. The upward velocities are created by the thermal force
propotional to $\nabla n\times(\nabla T_{i}\times\nabla T_{e}).$

##  5 Discussion

A physical mechanism has been proposed for the generation of spicules and
their three stage life cycle. There are two interesting questions related to
this phenomenon; one is how the density structure is created in chromosphere
with non-zero vertical velocity and the second is why this density structure
decays in corona after getting a certain height.

The process of taking birth of spicules at the bottom of chromosphere has been
considered as a sudden rise in local plasma density within a finite area as a
function of $(x,y)$ coordinates at time $t=0$ in the presence of gradients of
the temperatures of electrons and ions along vertical diretcion i.e. along the
z-axis. Since temperatures of charged particles have positive gradients along
z-direction in chromosphere $\nabla T_{j}=T_{j0}^{\prime}\hat{z}$, therefore a
net upward acceleration is produced by the thermal force which is proportional
to the cross product of gradients of density and temperatures $\nabla
n(x,y,t)\times(\nabla T_{e}+\nabla T_{i})$. The density structure attains a
final non-zero velocity at top of chromosphere after spending time $\tau_{ch}$
here and enters into transition region (TR). The steeper gradients of
temperatures along vertical direction in TR give larger upward push to this
density structure. The density hump passes through this thin layer of about
$100$ $km$ thick in a very short time $\tau_{tr}$. After gaining large
velocity it enters into corona where the upward force vanishes due to $\nabla
T_{j}=0$ and the solar gravitational force acting in the negative direction
reduces its final upward velocity to zero after time $\tau_{cr}$. The process
of the creation of density structures at the bottom of chromosphere continues
till the first such structure reaches the point of death in corona with zero
final velocity. At this moment, the process of creation of spicule stops. The
life time of the spicule is the sum of its time durations spent in the three
regions $\tau_{l}=\tau_{ch}+\tau_{tr}+\tau_{cr}$.

The above model for the genertaion of spicules and their life cycle is based
on the two fluid plasma theory. The two fluid plasma equations have been
extensivley used to study the dynamics of astrophysical and laboratory
plasmas. But it is not straightforward to obtain the exact solutions of this
coupled set of highly nonlinear partial differential equations. However,
during the past two decades, the two dimensional (2D) Saleem (2007, 2010,
2021) and three dimensional (3D) Saleem & Saleem (2022) exact analytical
solutions of these equations have been obtained for a restricted type of flow
which is longitudinally uniform with $\nabla\cdot{\bf v}_{j}=0$ assuming the
gradients of density and temperatures to be time-independent
$\partial_{t}n=0$. Recently Saleem & Saleem (2022), it has been shown that the
jet-like vertical flows can be generated at the surface of neutral fluid or
plasma by the thermal force proportional to $\nabla n(x,y)\times(\nabla
T_{e}+\nabla T_{i})$ where $\nabla T_{j}=T_{j0}^{\prime}\hat{z}$. This
solution can be applied to the transition region and corona during the life
cycle of the spicule. But to explain the generation of the spicule at bottom
of the chromosphere, we have to assume that density structure is created
during a finite time however small it may be. For this purpose, we assume that
the local density hump is created with a delta function in time $\delta(t)$ in
the chromosphere. Luckily, the exact solution of two fluid plasma equations
turns out to be similar to the one found for the case of constant density with
respect to time obtained previously. This new 2D analytical solution with
time-dependent density is very interesting because it can explain the process
of creation of plasma structures of different sizes and different scales in
astrophysical systems. The density function becomes a step function in time
and continuity equations have a source term $S(x,y,t)$ for $\partial_{t}n\neq
0$. It justifies the generation mechanism with consistent mathematical
formalism. At present, we have applied this solution to explain the creation
of spicules in chromosphere.

The spatial profile of density hump created in the chromosphere is similar to
the function $F(x,y)$ as shown in Fig. (1). The $N_{0}$ is a free parameter
and its value is assumed to be larger than the background density i.e.
$n_{0}<N_{0}$. Finite region of chromosphere has been chosen in the first
quadrant to show that in the presence of observed vertical gradients of
temperatures in choromosphere and transition region, all points within this
finite density structure i.e $0<x_{m},y_{m}$ get a net upward acceleration and
move in vertical direction. However, the accelerations and velocities of
different points $(x,y)$ have different values. The life time of the spicule
turns out to be a few minutes which is of the order of observed values. But
the estimated heights turn out to be larger than the observed ones. The
observed heights mentioned in literature are of the order of $10^{6}$ $m$ in
general while we obtain the values of the order of $10^{7}$ $m$. If the
dissipative effects are taken into account, then these values may reduce.
However, it should be noted that the combined thickness of chromosphere and
transition region is $(1.6)\times 10^{6}$ $m$ in our model and we have
included the total length of the structure from foot of the spicule in
chromosphere to the head in corona. It may happen that by the time, it reaches
its final destiny in corona, the lower parts had vanish already. On the other
hand, since the upward flow is not uniform, therefore the different layers of
flow within the structure give rise to resisitive forces which have not been
taken into account in the present analytical model based on exact solutions of
two fluid equations. The variation of density with respect to altitude may
also cause modifications in the shape and acceleration etc. Such
considerations are beyond the scope of the present work.

The thermal force mainly creates the ions vorticity $(\nabla\times{\bf
v}_{i})$ which is coupled with the generation of magnetic field. The generated
vorticity requires the ions to flow in the vertical direction, viz. ${\bf
v}_{i}=v_{i}(x,y,t)\hat{z}$. The ion velocity is related to the flow of
electrons through Ampere’s law, and hence the plasma jet-like flow along the
z-axis is created. The magnitude of the generated magnetic field ${\bf
B}(x,y,t)$ turns out to be very small compared to the constant ambient
magnetic field ${\bf B}_{0}$ because the scale lengths of gradients are large.
The Figs. (2) and (3) represent the behavior of physical quantities, but are
not according to physical scales. The Fig. (4) shows the variation of vertical
velocity within the spicule as a function of $(x,y)$ coordinates. This is a
general model which can explain the fundamental physical mechanism for the
generation of vertical plasma motions of different spatial scales in the solar
atmosphere. It can be applied to the generation of several different type of
structures such as threads and strands within coronal loops, solar mass
ejection and it also seems to add plasma in the solar wind at least partially.

Acknowledgement
Authors would like to thank the anonymous referee for his very valuable
suggestions and comments which have helped the authors in making reasonable
numerical estimates and improving the presentation of this research work.

## References

  * Agaoglou et al. (2024) Agaoglou, M., García-Garrido, V., Harlander, U., & Mancho, A. 2024, Physics of Fluids, 36
  * Antolin et al. (2021) Antolin, P., Pagano, P., Testa, P., Petralia, A., & Reale, F. 2021, Nature Astronomy, 5, 54
  * Aschwanden & Peter (2017) Aschwanden, M. J., & Peter, H. 2017, The Astrophysical Journal, 840, 4
  * Baran et al. (2022) Baran, B., Machaj, K., Malecha, Z., & Tomczuk, K. 2022, Energies, 15, 854
  * Beckers (1972) Beckers, J. M. 1972, Annual Review of Astronomy and Astrophysics, 10, 73
  * Bellan (2018) Bellan, P. 2018, Journal of Plasma Physics, 84, 755840501
  * Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Bol’shov et al. (1991) Bol’shov, L., Dykhne, A., Kowalski, N., et al. 1991, Handbook of Plasma Physics, North-Holland, New York
  * Bose et al. (2021) Bose, S., Joshi, J., Henriques, V. M., & van der Voort, L. R. 2021, Astronomy & Astrophysics, 647, A147
  * Brueckner & Jorna (1974) Brueckner, K. A., & Jorna, S. 1974, Reviews of modern physics, 46, 325
  * Carlsson et al. (2019) Carlsson, M., De Pontieu, B., & Hansteen, V. H. 2019, Annual Review of Astronomy and Astrophysics, 57, 189
  * Chen et al. (1984) Chen, F. F., et al. 1984, Introduction to plasma physics and controlled fusion, Vol. 1 (Springer)
  * D’Angelo & Motley (1963) D’Angelo, N., & Motley, R. 1963, The Physics of Fluids, 6, 422
  * De-Leon et al. (2020) De-Leon, Y., Paldor, N., & Garfinkel, C. I. 2020, Quarterly Journal of the Royal Meteorological Society, 146, 2096
  * De Pontieu et al. (2004) De Pontieu, B., Erdélyi, R., & James, S. P. 2004, Nature, 430, 536
  * De Pontieu et al. (2007a) De Pontieu, B., McIntosh, S., Hansteen, V. H., et al. 2007a, Publications of the Astronomical Society of Japan, 59, S655
  * De Pontieu et al. (2014) De Pontieu, B., Title, A., Lemen, J., et al. 2014, Solar Physics, 289, 2733
  * Dover et al. (2021) Dover, F. M., Sharma, R., & Erdélyi, R. 2021, The Astrophysical Journal, 913, 19
  * Goddard et al. (2017) Goddard, C. R., Pascoe, D. J., Anfinogentov, S., & Nakariakov, V. 2017, Astronomy & Astrophysics, 605, A65
  * González-Avilés et al. (2017) González-Avilés, J., Guzmán, F., & Fedun, V. 2017, The Astrophysical Journal, 836, 24
  * Hansteen et al. (1997) Hansteen, V., Leer, E., & Holzer, T. 1997, The Astrophysical Journal, 482, 498
  * Hasegawa & Mima (1978) Hasegawa, A., & Mima, K. 1978, The Physics of Fluids, 21, 87
  * Hollweg (1982) Hollweg, J. V. 1982, The Astrophysical Journal, 257, 345
  * Iijima & Yokoyama (2015) Iijima, H., & Yokoyama, T. 2015, The Astrophysical Journal Letters, 812, L30
  * Kadomtsev (1963) Kadomtsev, B. 1963, SOVIET PHYSICS JETP, 16
  * Kingsep et al. (1990) Kingsep, A., Chukbar, K., Yankov, V., & Kadomtsev, B. 1990, Electron Magnetohydrodynamics, 243
  * Klimchuk (2006) Klimchuk, J. A. 2006, Solar Physics, 234, 41
  * Kuźma et al. (2017a) Kuźma, B., Murawski, K., Kayshap, P., et al. 2017a, The Astrophysical Journal, 849, 78
  * Kuźma et al. (2017b) Kuźma, B., Murawski, K., Zaqarashvili, T., Konkol, P., & Mignone, A. 2017b, Astronomy & Astrophysics, 597, A133
  * Kuźma et al. (2019) Kuźma, B., Wójcik, D., & Murawski, K. 2019, The Astrophysical Journal, 878, 81
  * Lang (2001) Lang, K. R. 2001, The Cambridge Encyclopedia of the Sun (Cambridge University Press)
  * Lazarian (1992) Lazarian, A. 1992, Astronomy and Astrophysics, 264, 326
  * Martínez-Sykora et al. (2017) Martínez-Sykora, J., De Pontieu, B., Hansteen, V. H., et al. 2017, Science, 356, 1269
  * Martínez-Sykora et al. (2020) Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020, The Astrophysical Journal, 889, 95
  * Matsumoto & Shibata (2010) Matsumoto, T., & Shibata, K. 2010, The Astrophysical Journal, 710, 1857
  * Priest (1982) Priest, E. 1982, Solar Magnetohydrodynamics (Dordrecht: D. Reidel Publ. Co.), Ch
  * Priest (2014) —. 2014, Magnetohydrodynamics of the Sun (Cambridge University Press, p-31)
  * Saleem (1996) Saleem, H. 1996, Physical Review E, 54, 4469
  * Saleem (1999) —. 1999, Physical Review E, 59, 6196
  * Saleem (2007) —. 2007, Physics of plasmas, 14
  * Saleem (2010) —. 2010, Physics of Plasmas, 17, 092102
  * Saleem (2021) —. 2021, Physics of Plasmas, 28
  * Saleem & Saleem (2022) Saleem, H., & Saleem, Z. H. 2022, The Astrophysical Journal, 927, 72
  * Scalisi et al. (2021a) Scalisi, J., Oxley, W., Ruderman, M. S., & Erdélyi, R. 2021a, The Astrophysical Journal, 911, 39
  * Scalisi et al. (2021b) Scalisi, J., Ruderman, M. S., & Erdélyi, R. 2021b, The Astrophysical Journal, 922, 118
  * Secchi (1877) Secchi, A. 1877, Le soleil (Gauthier-Villars)
  * Shibata et al. (1982) Shibata, K., Nishikawa, T., Kitai, R., & Suematsu, Y. 1982, Solar Physics, 77, 121
  * Shibata & Uchida (1986) Shibata, K., & Uchida, Y. 1986, Astronomical Society of Japan, Publications (ISSN 0004-6264), vol. 38, no. 5, 1986, p. 631-660., 38, 631
  * Shimojo et al. (2020) Shimojo, M., Kawate, T., Okamoto, T. J., et al. 2020, The Astrophysical Journal Letters, 888, L28
  * Slemzin et al. (2014) Slemzin, V., Goryaev, F., & Kuzin, S. 2014, Plasma Physics Reports, 40, 855
  * Vranjes & Krstic (2013) Vranjes, J., & Krstic, P. 2013, Astronomy & Astrophysics, 554, A22
  * Weiland (1999) Weiland, J. 1999, Collective Modes in Inhomogeneous Plasmas: Kinetic and Advanced Fluid Theory (CRC Press)
  * Widrow (2002) Widrow, L. M. 2002, Reviews of Modern Physics, 74, 775
  * Woo (1996) Woo, R. 1996, Nature, 379, 321
  * Yokoyama & Shibata (1995) Yokoyama, T., & Shibata, K. 1995, Nature, 375, 42