Preprint
Review

This version is not peer-reviewed.

The Fluid Ionosphere

A peer-reviewed article of this preprint also exists.

Submitted:

30 December 2024

Posted:

31 December 2024

You are already at the latest version

Abstract
In this review paper the equations of motion describing the fluid dynamics of the ionosphere are constructed step by step, so that any master or post-graduate student may get familiar with the general theory of the “traditional” approach to Ionospheric Physics, in which chemicals forming the Earth’s upper atmosphere are represented as fluids in mutual interaction. The hypotheses on which the smooth-field fluid representation is based are discussed in terms of microscopic dynamics of the gas particles: this discussion is oriented to prepare the reader for the post-fluid approaches to the physics of turbulence, to be treated in forthcoming manuscripts. The fluid dynamical picture of the ionosphere is the most classical and conceptually simple one, and this makes it extremely widespread in terms of theoretical models and applications.
Keywords: 
;  ;  ;  

1. Introduction

The aim of this review is to present the general differential equations describing the Earth’s ionosphere system as composed by a certain number of classical fluids interacting one with the others in mechanical, thermodynamical and chamical way. Each of these fluids represents one of the various components of the ionosphere (ions, free electrons and neutrals). Mechanical interactions lead to energy and momentum exchanges. Thermodynamical interactions yield thermal exchanges of energy. Finally, the presence of chemical reactions allows for transformation of one fluid into another one via production and loss processes. Last but not least, elecromagnetic fields are generated by, and influence the dynamics of all those chemical species, many of which are electrically charged. The name we give to this system is generically ionospheric local medium (LIM), because its properties may depend on the position r and on time t. Here the equations of evolution of the LIM are constructed step by step in great detail, taking into account for all the aforementioned aspects (mechanical, thermodynamical, electrodynamical and chemical) within a classical fluid approximation: through the concept of fluid parcel the dynamics is introduced keeping a macroscopic classical point of view, that does need suitable hypotheses of Statistical Mechanics, but does not stem from a particular microscopic dynamics and calculations on it. Great attention is payed to carefully expaining what all these approximations mean and what their implications are.
The review is organized as follows.
Section 2 is devoted to the rigorous definition of the system at hand and to the explanation of the approximations used.
In Section 3 the equations of motion are written down fully, including mass, momentum and internal energy balance for each constituent, plus Maxwell Equations for electromagnetic fields. Also, constitutive hypotheses are presented, to make the system closed mathematically. The dynamical system is formed by coupled partial derivative equations (PDE) plus constitutive hypotheses as contraints.
Section 4 uses the results obtained in Section 3 to develop some examples and applications to concrete ionospheric phenomenology: in particular, one treats the LIM as a gaseous conductor, and writes down its electrological relationships.
Finally Section 5 underlines again the main points explained in this review, and traces its further developments, most of all in terms of what happens when the conditions for fluid approximation are removed, then introducing idea of post-fluid formalism [1,2].

2. The System

The system of interest here is the Earth’s upper atmosphere, considered above the minimum quote at which the ionization due to Sun’s radiation or cosmic rays is so important to affect the propagation of electromagnetic waves [3]: this might appear a technology-centreed definition, but actually this “lower limit” coincides with the quote at which the atmospheric medium is active from an electromagnetic point of view, that is a rather intrinsic physical concept. All in all, our reasoning holds from approximately 60 km height up to the magnetopause of the Earth, in many different declinations. Actually, some authors refer to as “ionosphere” only to the electrically charged chamicals, i.e. ions and free electrons. Yet, as the equations of motion of these charged species do involve neutrals, and as the equations of neutrals can be expressed formally as those of the charged species, we are considering them as part of our system.
The reference frame used here is mainly the one co-rotating with the solid Earth, referred to as Solid Earth Reference Frame (SERF); sometimes it will be of use to introduce another auxiliary reference frame, namely the Neutral Wind Reference Frame (NeWReF), see below. Last but not least, the unit system used will always be the MKSA system.
The system of interest is neither a closed nor an isolated one [4]: it interacts mechanically and radiatively with what lies below it (lower atmosphere, lithosphere [5] and hydrosphere) and with the space environment surrounding the Earth (i.e., the solar wind, the Sun’s radiation and cosmic rays). The system does also exchange matter with what is below it and above it.

2.1. Smooth and Deterministic: The Continua

At a certain point of their history, physicists accepted that the world is made of very tiny systems, as molecules, atmos, electrons, and after all elementary particles, that render matter granular. However, since a couple of centuries at that moment, they had already developed a picture of macroscopic material systems as continua, i.e. extended distributions of matter occupying finite portions of the physical space and experiencing and applying forces all over their volume, and surface. From the invention of Calculus, physicists had given a precise mathematical meaning to the Latin proverb Natura non facit saltus, taken by Leibniz as a manifesto for modelling physical systems: any physical quantity should have a C -dependence on the variables it is function of. This “phylosophy of non-roughness” applies to the description of continua stating that if a macroscopic continuum C occupies a volume D t C at a certain time t, any positional quantity ψ r , t describing its state (e.g., mass density, velocity, pressure...) must be a smooth field all over D t C :
ψ C D t C × R , K ,
being the factor “ × R ” encoding the time dependence of ψ and K the mathematical set of values of ψ r , t (numbers, vectors, tensors...). In Figure 1 the assessment (1) is cartoonised: a greenish continuum C evolves from an initial configuration C t 0 to a later one C t through a supposed smooth evolution; this means that any local quantity ψ describing it undergoes precisely the condition (1). Actually, considering the border of the sets D C the region where the local proxies of the continuum stop being defined, or may be understood to go from a finite value to zero “within no distance”, one should replace (1) with something more precise like ψ C D t C D t C × R , K (this does not, however, lighten the request (1)).
Despite its high requests on the mathematical description of matter reality, this picture has been able to give very satisfactory theories, as hydrodynamics, magneto-hydrodynamics, elasticity theory and many more [6]. Note that the picture satisfying (1) is also a deterministic picture, as the smooth time-behaviour of ψ  excludes noise from its time evolution, which is another very strong requirement.
Even if the continuous picture of macroscopic systems is very beautiful, yet at a certain point the granular nature of matter was experimentally proved, and in order to reconcile the latter with the smooth picture of macroscopic bodies some reasoning had to be introduced, namely the so called continuous limit of pointlike particle systems. The fluid dynamics (FD) of the ionosphere is one important product of this “ideology”, that has an elegant and effective resoning, but relies on precise and not-that-universal properties of the microscopic systems at hand.
The central concept of the continuous limit contruction is considering “small” portions of the continuum d C , referred to as parcels, and make the following reasoning (the term “small” has to be specified later). In such a d C there exist N d C elementary components of the macroscopic matter system, each of which is moving with a certain velocity. In general, these N d C components may or may not have all the same chemical nature, hence the same mass: the local representation of the neighbourhood d C within the continuum is obtained considering the centre-of-mass P of those N d C components, to which the velocity of the centre-of-mass of the N d C particles is attributed as V P . The local representation of the continuum at P is constructed starting by saying that the parcel of continuum d C P moves with velocity V P . Figure 2 represents a cartoon of this idea. Recalling the definition of centre-of-mass, we may say that the point P is attributed the mass-average velocity of the N d C P particles within d C P , i.e. the total linear momentum of the matter in d C P .
The bridge between the micro-physics of particles and macro-physics of continua was traced by Statistical Mechanics, in which the macroscopic appearance of the world is conceived as the statistical sum of all the microscopic, fluctuating contributions of particles. Boltzmann’s Kinetic Theory paved the way for a quantitative study of those cumulative effects of what the micro-components are doing: this makes it possible to recover the local variables ψ in (1) as suitable momenta of the statistical distributions of the N d C P particles in d C P [7]: this is how mass density, bulk velocity, pressure, tension ... are defined in FD. This sounds quite intuitive, as we humans simply accept that our senses, and our classical instruments, simply sense the world of elementary particles “defocussing” their tiny details, and saving only a coarse-grained, average picture of it. In this reasoning, when one looks at the N d C P particles within the parcel d C P of centre-of-mass P, the position r mentioned in the FD fields ψ r , t is the position vector of P: in what follows, r and P are used with the same meaning.
All in all, the macroscopic picture ψ would be obtained via a “statistical” way:
ψ P , t = ψ micro d C P , t ,
being ψ micro some physical quantity depending on the single particle variables, and . . . d C P , t a sign of suitable statistical averaging(whatever this may mean at this early point) performed within d C P at time t. In order to rely on a statistical procedure as (2), one needs a first important ingredient to the very existence of the ψ s: the number N d C P has to be so huge to point towards the Thermodynamic limit
N d C P = O N A .
This is the first tough requirement to go from particle physics to FD, that had become a necessity once molecules, atoms and all their tinier fellows were discovered. In any statistical theory, due to results as the Central Limit Theorem [8], there is the need of having enough samples to be able to neglect stochastic fluctuations, that is precisely what one needs to do in order for ψ to be a deterministic (hence, time-smooth) process. Note: the condition (3), that may sound definitely sensible considering a small glass of water to contain N A s of H 2 O molecules, is not at all obvious when one would like to represent the physics of galaxies (made of “few” solar systems) or even a finite portion of very dilute space plasmas (made of “few” p and e ) with FD, in which the constituents are all but N A s of bodies. In those cases, the theory of small thermodynamic systems, large fluctuations and extreme events have to be invoked [9].
Once we have defined the continuum quantities ψ as non-fluctuating fields on D t C , one has to have them smooth, as prescribed by (1). This means that, if one has two “nearby” points P and P within D t C separated by some “small” displacement δ r , as in Figure 3, one has to be able to calculate the non-fluctuating value ψ P , t as ψ P , t = ψ micro d C P , t , the non-fluctuating value ψ P , t = ψ micro d C P , t , and then find the incremental limit between them with respect to the P P = δ r distance to exist:
lim δ r 0 ψ P , t ψ P , t δ r < + δ r .
Equally, there must exist also the limit
lim δ t 0 ψ P , t + δ t ψ P , t δ t < +
defining the instantaneous rate of change with time of ψ P , t = ψ micro d C P , t . Actually, the two requirements (4) and (5) demand for differentiability of ψ in space and time, while the requirement (1) is extremely stronger because it means ψ to have infinite space and time derivatives: clearly (1) imply both (4) and (5). More cleanly, one could replace the two requirements (4) and (5) getting the whole (1) by assessing
m , n N / ψ P , t ψ P , t = O δ r m , ψ P , t + δ t ψ P , t = O δ t n :
the differences between ψ at two (macroscopically) nearby points or at two nearby times must scale to zero as positive integer powers of the distance or of the time lapse considered. Clearly, this (6) rules out any turbulent regime of classical fluids [10]. In (4), (5) and (6) the symbol . . . refers to a length in the set of the valus of ψ P , t , for instance the absolute value for real numbers, or Pitagora Theorem norm for three-dimensional vectors.
The construction of the continuum representation of a macroscopic system of matter is all pivoted on the definitions (2) satisfying the conditions (6), starting with the definition of a velocity field  V P , t as described before and sketched in Figure 2, and a mass density field  ρ P , t , defined for the matter in d C P as the ratio M d C P d V , being M d C P the total mass of particles forming d C P , and d V the parcel volume. Defining these ρ P , t and V P , t helps us explaining how “small” the “parcel” must be taken around P: reference is made to the idea of macroscopic infinitesimal portion [11], for which one starts to ideally shrink the volume in C around P and applies the definitions to smaller and smaller volumes d V ; at a certain size of d V one stops seeing ρ P , t and V P , t fluctuating, while go on with shrinking too much, fluctuations re-appear. At the beginnig of this “plateau” of the ψ values with decreasing d V one stops including in the parcel regions of the continuum in different local conditions; when fluctuations re-appear at too small values of d V , the number N d C P is not satisfying the condition (3) any more, and the granularity of matter is unwantedly appreciated in terms of fluctuations.
Concluding this rapid introduction to continua, that is the necessary foreword to FD, we refer again to Figure 2: the set of particles in d C P have their 3 N d C P degrees of freedom, three of which have been taken into account considering the position of P and the velocity V P : these are the centre-of-mass variables [12]. What about the other 3 N d C P 1 degrees of freedom, namely the variables relative-to-P, e.g. the difference between the position of each of the N d C P particles of d C P and that of P, and their canonically conjugated momenta?
Treating also this huge number of relative variables is necessary to completely describe che phase space of the N d C P particles in the parcel, but despite it is sensible to use the per se deterministic variables r P , p P for the centre-of-mass P of the parcel, being p P = M d C P V P the total momentum of d C P , the 3 N d C P 1 = O N A relative variables of the parcel will be treated in a statistical way. That is: while V P is already a field ψ r , t as those involved in equations (1) and (2), the positions q i rel = def r i r P of the i = 1 , . . . , N d C P particles in the parcel relative to P and their conjugate momenta become the subject of the d C P -local statistical mechanics. As the dynamics in this relative phase space
Γ rel = q 1 rel , p 1 rel , . . . , q N d C P 1 rel , p N d C P 1 rel
is treated statistically, becasue dim Γ rel = O N A , the statistical distributions on Γ rel become the completion of the description of the exactly N d C P particles in d C P .
On a theoretical ground, describing statistically the relative variables within the parcel vmay consist of the dynamics of the Gibbs distribution in Γ rel [7,13]; in practice, one always assumes the particles in d C P to be represented by some kinetic theory that will cut the BBGKY hierarchy in Γ rel to 1-particle distributions. Moreover, as the only well established statistical mechanics is equilibrium thermodynamics, the d C P -local kinetic theory will converge to the definition of local thermodynamic variables attributed to the point P D t C of position vector r , as new smooth fields ψ r , t . All in all, the 3 N d C P microscopic degrees of freedom of the particles in the parcel will turn into the 3 ones of the centre-of-mass plus an equilibrium thermodynamics for the 3 N d C P 1 relative coordinates, representing the inner structure of d C P , converging however to some other local smooth continuum system variables ψ r , t (e.g., the pressure, or density, and temperature, see SubSection 3.4 and Section 3.5).
With this, the FD philosophy of the LIM is established; before moving on focussing on the system composition, let us resume it in items:
  • the system is formed by parcels with thermodynamic numbers of particles (3);
  • the statistical mechanics of these particles converges to space- and time-differentiable emerging variables (6);
  • the system is locally in thermal equilibrium, hence described via a kinetic theory converging to a local equilibrium thermodynamics.
The requirement 1 makes it possible to have well defined non-fluctuating statistics tout court; the requirement 2 leads to the feasibility of deterministic smooth FD at all. Last but not least, the requirement 3 is the path leading to FD as we know it today. The thorny debate about physical regimes when those requirements do not appear to be met deserves future manuscripts.

2.2. The System: Matter and Interactions

Let us now introduce how all this theoretical background is applied to the Earth’s ionosphere.
The system is considered to be formed by S chemical species, indicated as X I , with I = 1 , 2 , . . . , S . For example, molecular species as O2, N2, CO2, NH4, but also atomic species as O and N, are in this list. If the neutral species considered are
O 2 , N 2 , CO 2 , NH 4 , O , N
and the ion species are
NO + , O 2 ,
then the example composition reads:
X 1 = O 2 , X 2 = N 2 , X 3 = CO 2 , X 4 = NH 4 , X 5 = O , X 6 = N , X 7 = NO + , X 8 = O 2 , X 9 = e .
In the ionosphere, these S chemicals are in the state of gases. Each of the different gases is an extended system of matter represented as a continuum C I in the way discussed before, i.e. with parcel dynamics: at each point P of the space a parcel d C I P is considered for the I-th chemical, with a centre-of-mass velocity V I P and the I-th relative phase space Γ rel I P described via statistical mechanics, in general, more practically via P-local equilibrium thermodynamics.
This is only the matter part; there is a field part too, made of those fundamental force fields acting on the atmospheric gases: gravity and electromagnetic fields. Due to the average speeds of the subsystems in the SERF, and to the energy of the processes treated, here we represent everything in terms of classical physics, involving only small sips of Special Relativity, here and there, when the electromagnetic fields are to be transformed from a first observer to another one. Gravity will not be a dynamical variable of the system!
Every component C I experiences the action of the gravitational field defined around the Earth, which includes the field of properly massive origin, the inertial positional forces and the time varying field due to the tidal action of other celestial bodies. From a General relativistic point of view one should include Coriolis force too in the gravitational field, while here the Coriolis force will be kept as separated from the whole “effective” gravity. In general one can write:
g = g 0 + g Ω + g tide .
Coriolis forces are thus due to Coriolis acceleration field
g C I = 2 Ω × V I ,
being Ω the rotation vector of the Earth as seen in a frame in which this Coriolis force is not experienced, and V I the fluid velocity of the I-th species, constructed as explained before.
Electromagnetic fields are experienced only by charged continua (in the example (7) the ions C 7 and C 8 and the electrons C 9 ), and for these particles their strength is much more important than that of gravity. In general there is the geomagnetic field B 0 , that is considered as a rigid external field, and some non geomagnetic, perturbing magnetic field B temp , due to temporary phenomena (e.g., the geomagnetic storms), so that one has
B = B 0 + B temp .
Many electric fields of different origin are also present, and this will be clarified later.
Quite the opposite to what happens to gravity, the electromagnetic fields acting on the C I s undergo important back reactions due to charge displacements, net currents and dynamoes produced by the motion of these fluids. The “complete” dynamical system we will deal with will then include all these fluids C I plus the fields E and B influenced by the atmospheric motions (it should be stated that the mechanics of the C I s are only a part of the source of the atmospheric E and B , since B is generated, in its component B 0 in equation (9), by some conducting lavic movement inside the solid Earth [14]).

2.3. Fluid Matter Variables

Every single chemical component C I is described as made of particles of mass m I and electrical charge q I , we will use ρ I to indicate mass density, N I for numerical density and δ I for charge density. These quantities are related as:
ρ I = m I N I , δ I = q I N I .
Of course, ρ I r , t , N I r , t and δ I r , t are proper fluid fields as the ψ I r , t discussed in SubSection 2.1.
The components involved can be neutral species
q I = 0
or charged ones
q I 0 .
In ionospheric studies the use of N I is often preferred instead of ρ I , because the numerical densities of particles enter the typical expressions of cross sections, determining the production and loss rates of the C I s in chemical processes, and also because the expressions of many measurable quantities are entered by the numerical density N e of free electrons. For what concerns δ I one must say that the functional relationship (10) renders it practically equivalent to N I or ρ I for the charged species, while when q I = 0 it simply vanishes.
About these δ I s, in general a constraint of neutrality of the LIM
I δ I = 0
may exist, besides the global one:
R 3 I δ I r , t d 3 x = 0 .
The constraint (13) is a standard in many applications of ionospheric physics, even if particular regions, as the dawn and dusk terminator, may violate it.
Local thermal equilibrium is assumed for each C I : then, temperature fields T I r , t will come into the play; it is useful to remember that the different chemicals have different temperatures, they are not in thermal equilibrium with the others; for example, electron temperature has been measured as about 900 K , while positive ions were found around 2500 K in the polar E layer at about 150 km height, when intense electric field are present, as reported in [15]. This should not be surprising, despite the different C I s are in thermal contact with each other: as they are moving, floating and interacting, these gases are in a well different condition from that of fluids in a calorimeter. This also keeps the system far from a global thermodynamic equilibrium.
The continuous approach, which is largely valid for studying plasma and neutrals of the “standard” ionosphere, can not be used to describe the rain of high energy particles entering the Earth’s atmosphere from the space, because of the hypothesis of local thermal equilibrium which is lacking for them.

2.4. Dynamical Variables for Matter

A fluid occupying a certain volume in the space is fully described when point by point and time by time its velocity field, mass density and tensions are given [16]; this means that for each constituent C I one should use the S O 3 -vector
V I = V I x r , t , V I y r , t , V I z r , t
(bulk velocity given by the centre-of-mass velocity of the parcel with the centre-of-mass in r ), the S O 3 -scalar
ρ I = ρ I r , t
(bulk mass density) and the distinct elements of the stress tensor
T I i j = T I i j r , t
with i , j = x , y , z , a two rank S O 3 -tensor. Usual fluids as those at hand have symmetric stress tensor, so that one has 6 independent components T I i j r , t . This ρ I , plus V I , plus T I i j give all the fields to describe each continuum C I composing the atmosphere in the Eulerian picture. However, tensions are not further variables to be added to V I and ρ I : the T I i j s are in general functionally dependent on the velocity field V I , and “fluids” are classified into different species according to the mathematical form of the dependence
T I i j = T I i j V I , . . . .
The tension dependence on velocities is restricted in general to a particular subdivision of T I i j , that can be written as the sum of a pure positional addendum and on another term depending on the V I s too:
T I i j = p I r δ i j + τ I i j V I , r , . . . .
The continuum C I is defined as fluid if
V I = 0 τ I i j = 0
holds. In the very particular case in which
τ I i j V I = 0 V I
one speaks about ideal fluids. The definition (20) means that the system is not able to exert any off diagonal stress when it is at rest. The statement (21) indicates as ideal fluids those fluids unable to exert off diagonal stresses also when they are moving.
The S O 3 -scalar p I r in (19) is the hydrostatic pressure. The case of perfect fluids well represents the C I s in very dilute regions [17]. Note that since p I is a scalar under the rotation group it must be related to the trace of T I i j , so that τ I i j is traceless.
In order to describe energy exchanges the thermodynamical nature of the system must be taken into account; then another quantity is introduced, the local temperature T I r , t , which is also an S O 3 -scalar. The full set of dynamical variables we have introduced to describe each C I is then given by the density ρ I , the velocity V I , the pressure p I , the τ I i j tensor and the temperature T I , giving a set of 12 unknowns in the equations of motion. In SubSection 3.4 the roles of T I and p I will be clarified from a fully mechanical point of view, as coarse grained proxies of the microscopic degrees of freedom of the matter of the fluid at hand, contained in the volume d V P , in terms of variables relative-to-the-centre-of-mass of the “parcel”.
The various fluids C I are generally classified and studied as follows [18]:
  • Neutral particles. Various molecules or atoms (e.g. O, N, O2, N2, CO2, NH4 and so on) with zero electric charge: these ones are the major part, both numerically as well as massively, of all the high atmosphere, up to the level of about 1000 km height, starting point of a region where practically all the matter is ionized. From that point on, one speaks about plasmasphere.
    If each neutral species is indicated with the index h = 1 , 2 , . . . , the speed of their centre of mass is point by point some field U , defined as follows in terms of partial matter densities ρ h and velocities U h :
    U = 1 ρ h ρ h U h , ρ = h ρ h .
    This velocity undergoes a proper dynamics that may be attributed to the neutral matter as a unique fuid C 0 .
  • Negative ions. These are the species which have captured one or more electrons. These particles are very rare and their presence is practically negligible in general. Nevertheless, one should consider them when dealing with the lower part of the ionosphere, as the region D and the lower E [19]. Negative ion abundance is reasonably non negligible only under 95 km, even if in particular conditions this can be false. Sometimes a coefficient λ Y is defined for each negative ion Y as
    λ Y = N Y N e ,
    giving the measure of the contribution of the species Y to the ionospheric negative charge, relative to the electron one.
  • Positive ions. These are the neutrals that have lost one or more electrons: actually, the 2- or 3-valent positive ions are very rare, so we can restrict to the study of 1-valent positive ions. When the constraint (13) holds and no negative ions exist, positive ions have a numerical density equal to that of free electrons.
  • Electrons. Electrons are the most important element for ionospheric phenomenology as far as electromagnetic disturbances are concerned, because of their very small mass that gives them a great mobility [3], making them very effective in producing local and travelling electromagnetic fields.
In general one could represent in a Lagrangian way the fluid motion as the evolution of a manifold, a domain continuously filled in with matter, and describe its motion as the transformations undergone by this manifold. The manifold can stretch, roll, it can be deformed by tensions, its centre-of-mass moves under the action of external forces. Basically one physical portion of it can be identified as the evolution of a fixed portion of the initial continuum [16]. This way of representing the evolution of a continuous system has led to a Hamiltonian description of fluids alternative to the local approach (for reference see [20,21,22] and [23]).
In the absence of chemical reactions the mass of the parcel d C I remains constant:
NO CHEMISTRY d d t m d C I = 0 ,
being d d t the operator making the time derivative along the motion [24] . What happens if chemical reactions take place? When chemical reactions enter the fluid dynamics, it is convenient to treat the systems as non constant mass systems, writing:
CHEMISTRY d d t m d C I = G I d C I .
This G I d C I is the variation rate of the mass of the parcel as following the parcel motion, and it can be either positive or negative, according to production or loss processes.
The description of the fluid d C I is referred to as Lagrangian or substantial, because the quantities used refer to the degrees of freedom of a material portion of the system C I ; on the contrary, fluids are often described in the Eulerian or local way, in which one describes the situation via quantities thought of as properties of the spatial point at hand [16]. When some quantity G is attributed to a fluid, the relationship between the Lagrangian version G L and the Eulerian version G E of it is given in terms of time derivatives, and reads
d G L d t = G E t + V · G E ,
being V the local speed of the fluid to which G is attributed (for ionospheric continua, the G s and V in (26) have a chemical species index I)
It has already been stated that the extended system d C I r has centre-of-mass (CoM) variables and relative-to-CoM (R2CoM) variables as explained in [12] and references therein. The evolution of the CoM variables will be governed by Newton’s principle as explained in SubSection 3.3 below. For what concerns the R2CoM variables,a statistical approach is used for them in ordinary fluid dynamics, see better SubSection 3.4.

2.5. Electromagnetic Dynamical Variables

The mathematical tools describing the system are completed when the variables used for the electromagnetic fields are given. At each point r in the space around the Earth an electric vector E and a magnetic vector B can be given, both time as well as position dependent:
E = E r , t , B = B r , t .
These electric and magnetic fields are obviously governed by Maxwell equations. While there exists a Lagrangian description of fluids with a straightforward interpretation, electromagnetic fields do have only an immediate Eulerian dynamics: this is why all the d C I -dynamics is reformulated in local terms, through the translation (26).
In general, the electromagnetic field in a point of interest will be the superposition of many fields of different origin.
The magnetic field B is given by the superposition of a “permanent” geomagnetic component B 0 , that is generated by sources in the interior of the solid Earth, and other “temporary” terms, indicated collectively as B temp in (9). B 0 is not constant on the long term, because it has been proved to vary: roughly speaking, it is some dipole magnetic field that slowly rotates and changes its shape. At Italian latitudes and longitudes, for instance, it is representable as a rigid dipole migrating of about one prime per year. Actually, this dipole B 0 is not rigid, and its motion should be studied in a complex theoretical way (see [14] and references therein). However, the variability of this B 0 takes place on times which are too long compared to the processes usually discussed in LIM modelling.
The inclusion of a field in the system at hand as a dynamical variable depends, by the way, on the inclusion of its sources: since we decided not to include the currents inside the solid Earth in our system, we will regard the geomagnetic field generated by them as an external constant field  B 0 . On the contrary, the part B temp can be affected by ionospheric or magnetospheric sources.
The electric field E is generally produced by charge separations related to the motions of the ionospheric plasma. In many cases this E can be considered as slowly varyingwith time, and treated as electrostatic, and so a potential ϕ giving E = ϕ is introduced.

3. Equations of Motion

The equations ruling the motion of the atmospheric system are essentially those of classical Newtonian physics, plus thermodynamics, plus classical electrodynamics.
For the electromagnetic fields the four Maxwell equations will be used directly, while for what concerns the matter part, conservation principles applied to the parcel d C I r , t dynamics will drive the reasoning.

3.1. Mass Balance

The mass balance is written for the I-th species by considering that each of the many species present in the atmosphere undergoes creation and destruction due to chemical and photochemical processes. Neutrals are “destroyed” in ionization and “reconstructed” due to recombination of electrons and positive ions, or detachment of electrons away from the negative ions. In the same way free electrons are created in ionization together with positive ions, in the detachment processes of negative ions, and destroyed in general in attachments and recombinations. Finally, positive ions are created together with electrons in ionization processes, and destroyed in recombinations.
In general the mass balance can be obtained by elaborating equation (25): if the parcel of the I-th element is some d C I which occupies a volume d V , and the mass density is ρ I , one can write:
m d C I = ρ I d V d d t ρ I d V = G I d C I .
The total derivative is applied as
d ρ I d t d V + ρ I d d t d V = G I d C I ,
and all we have to do is to write d ρ I d t d V + ρ I d d t d V as a function of the matter degrees of freedom presented in SubSection 2.4: the physical quantities we will involve in our discussion are the density ρ I and the velocity V I of the parcel.
The density is already involved explicitly in (28); one has to evaluate d d t d V in terms of the parcel velocity. The variation of the measure d V can be calculated considering that the measure of d V changes due to the motion of the points on its surface d V , so that this time derivative is the outgoing flux of this speed:
d d t d V = d V V I · n ^ d 2 S .
When Stokes theorem is applied to this integral one has
d d t d V = d V · V I d 3 x .
Considering the smallness of d V , it is possible to write:
d d t d V = · V I d V .
Placing (30) into equation (28), we can write:
d ρ I d t d V + ρ I · V I d V = G I d C I
The mass source term G I d C I can be imagined as proportional to the total mass of the parcel as G I d C I m d C I = ρ I d V = m I N I d V , because each particle in the parcel may undergo the same chemical reactions with the same probability distributions, so that G I d C I m I d V may be written as:
G I d C I = m I Γ I d V ,
being m I the mass of the single particle of the I-th species and the source field Γ I suitably defined in (32). This expresses equation (31) involving the Lagrangian derivative of ρ I
d ρ I d t + ρ I · V I = m I Γ I ,
and going from the Lagrangian to the local expression of d ρ I d t , one has:
d ρ I d t = t ρ I + V I · ρ I t ρ I + V I · ρ I + ρ I · V I = m I Γ I .
The total space derivative
V I · ρ I + ρ I · V I = · ρ I V I
is recognized, then equation (33) turns out to read:
t ρ I + · ρ I V I = m I Γ I .
It is useful to write the source term Γ I as the difference
Γ I = P I L I
between the local rate of production P I of the I-th species, and the corresponding rate of loss L I . Equation (34) is then re-written as:
t ρ I + · ρ I V I = m I P I L I ,
Each species has its own loss and production term in (36), but these are not all independent; their forms encode the chemical processes taking place in the ionosphere, that have the elementary constituents of the continua C I as reagents or products. These reactions undergo some conservation relations, and the source terms in (36) must satisfy those conservation laws. The matter of the total system
C = I C I
must have constant mass (under the hypothesis that no nuclear reaction takes place, and considering not to discuss the loss or gain of chemicals due to some matter exchange with the near-Earth space). This means that there must exist a total mass balance equation with zero source,
t ρ + · ρ V = 0 ,
where ρ is the total mass density of the multi-fluid system and V is “the right vector” to take part to equation (38) (being the “right vector” V to be defined few lines below). The way of deducing the total equation (38) from all the particular equations (36) is straightforward, it suffices to sum them all and define V accordingly, that will turn out to be rather intuitive. Summing all the equations (36) over the chemical species, one gets
I t ρ I + · ρ I V I = I m I P I L I .
Total quantities are recognized straigthforwardly: firs of all, the total density
I ρ I = ρ ,
while transport terms give the right addendum in (38) when V is defined so to satisfy:
I ρ I V I = ρ V .
This is a well known quantity in mass kinematics: V is as the “speed of the centre of mass” for all the fluids. Due to these definitions, we have
t ρ + · ρ V = I m I P I L I .
If equation (41) has to coincide with equation (38), enforcing the mass conservation, we must have
I m I P I L I = 0 ,
that is the mass conservation condition on the production and loss terms P I and L I . The equation (42) is a basic condition for aeronomical chemistry, claiming that chemicals are not created or destroyed “out of the blue” or “into nothing”, but simply turned the ones into the others. The only external contribution to this apparently closed system balance is the supply of energy from solar photons, that are destroyed when photochemical processes take place. Yet, these photons will not be recognised as a futher species of chemicals, essentially as long as chemical, and not nuclear reactions take place, so the energy of these massless photons will not be translated into mass.
The constraint (42) is not the only condition one has to fulfill when writing the production and loss terms. Conservation of the electrical charge must be taken into account, as well. In order to get the charge conservation condition, let us go back to (36) and divide it by m I , that is clearly a constant parameter. What one gets is
t ρ I m I + · ρ I m I V I = P I L I .
Conventions (10) give
t N I + · N I V I = P I L I ,
that is the balance equation for number densities. This equation (43) is the key continuity equation for any additive quantity held by the single fluid species: an example of such quantities is indeed the electric charge. Equation (43) is multiplied by the I-th electric charge
t q I N I + · q I N I V I = q I P I L I
and definitions (10) are turned to again, so that:
t δ I + · δ I V I = q I P I L I .
If δ is the total electric charge density of the whole system of fluids, then one would like to get
t δ + · J q = 0 ,
and this is obtained by summing all the equations (44) as I t δ I + · δ I V I = I q I P I L I and recognizing
I δ I = δ , I δ I V I = J q .
This gives
t δ + · J q = I q I P I L I ,
so that if equation (47) has to be equivalent to electric charge conservation, the following condition on the P I s and L I s
I q I P I L I = 0
must be realized.
In general, if an extensive quantity H exists and has density η I for the I-th continuous species, and total density
I η I = η ,
one can define its total current as
J H = I η I V I .
The single species balance equation for H reads
t η I + · η I V I = h I P I L I ,
being h I the amount of H carried by the single particle of the I-th species. This yields
t η + · J H = I h I P I L I ,
giving finally
t η + · J H = 0 I h I P I L I = 0
as a “total H conservation condition”.

3.2. Chemical Space Formalism

A particularly suitable symbolism for fluid with dynamics, loss and production processes, is referred to as vector chemical space formalism.
There is one equation (36) for each chemical species, let S be the total number of chemical species. Let us rearrange the mass conservation to involve the numerical concentration instead
t N I + i J I i = P I L I , I = 1 , . . . , S , J I = N I V I ,
then let us think that these S equations are the components of the following system:
t N ̲ + i J ̲ i = P ̲ L ̲
(where N ̲ , J ̲ i , P ̲ and L ̲ are S-component real vector in the “space of chemicals”: in equations (54) and (55) the summation over i = x , y , z is intended in the terms i J I i and i J ̲ i . Summations take place over repeated indices in contravariant positions). The local conservation constraint (53) for any conserved quantity H is read as an orthogonality relationship in R S between the collection of the values h ̲ and the collection of the source terms:
h ̲ · Γ ̲ = 0 , Γ ̲ = P ̲ L ̲ .
The equation (55) is referred to as species-vector continuity, and its components are coupled via the source term Γ ̲ that satisfies (56). In fact, this source term depends on the densities
Γ ̲ = P ̲ N ̲ ; r , t L ̲ N ̲ ; r , t ,
and maybe on the time and position explicitly, so that equation (55) all in all reads:
t N ̲ + i J ̲ i = Γ ̲ N ̲ ; r , t .
In order to understand equations (57) in a general way, we must classify explicitly the properties of the function Γ ̲ N ̲ , . . . , its possible symmetries and the simplifications. Turning back to the components t N I + i J I i = Γ I , we can appreciate the various levels of coupling between different species by expanding Γ ̲ in the densities, via a simple Taylor series approach, e.g. up to the third order:
Γ I N ̲ = Γ I 0 + Γ I N J 0 N J + 1 2 2 Γ I N J N H 0 N J N H + + 1 6 3 Γ I N J N H N K 0 N J N H N K + O N 4
(the following reasoning makes an exact sense as long as Γ ̲ N ̲ has a convergent Taylor series, otherwise it is just an approximation). When such expressions are put into (54)
t N I + i J I i = Γ I 0 + Γ I N J 0 N J + + 1 2 2 Γ I N J N H 0 N J N H + 1 6 3 Γ I N J N H N K 0 N J N H N K + O N 4
one discovers that the various terms involve zero, one, two, three, ... species of particles, and thus correspond to chemical processes with zero, one, two, three, ... chemical involved. It is not said that in the n-th order term n different species are involved: in fact, there are terms in Γ ̲ in which a given N I can appear with a degree higher than one. These different momenta of Γ I in the space of species must form a function which satisfies the constraint (56): putting (58) into it, one has
I q I Γ I 0 + Γ I N J 0 N J + 1 2 2 Γ I N J N H 0 N J N H + + 1 6 3 Γ I N J N H N K 0 N J N H N K + . . . = 0 .
The various terms must be zero each, because we are setting to zero a polynomial in the components of N ̲ . For the zeroth order one has
I q I Γ I 0 = 0 ,
which essentially means the impossibility of creating charge in the empty space, with N ̲ = 0 . One can see in Chapman ionization theory in [3,17] an explicit application of this. For example, if our model atmosphere is represented by the three species X, X + and e , equation (61) reads
0 · Γ X 0 + e Γ X + 0 e Γ e 0 = 0 ,
that is
Γ X + 0 = Γ e 0 .

3.3. Momentum Balance (Newton’s F = m a )

The fluids C I composing the system described at the beginning of this Section are made out of point particles that essentially obey to Newton’s principles. In the approximation described in SubSection 2.3 the initial condition Cauchy problem
F = m a , r 0 = r 0 , v 0 = v 0
for a particle of mass m can be transported to the continuous system of the I-th constituent as the Navier-Stokes equation [7], that is exactly Newton’s Second Principle, or equally the linear momentum balance, written for the parcel.
The linear momentum of the parcel is some p d C I , so we have
d d t p d C I = d F I .
If its mass is referred to as m I d V and we use the symbol d F I to indicate the force experienced by this d C I , we can write:
d d t m I d V V I = d F I ,
stressing again that V I is the velocity of the centre of mass of the parcel d C I at hand. From equation (25), this can be written as:
G I d C I V I + m I d V d V I d t = d F I .
The mass of the parcel reads
m I d V = ρ I d V ,
and considering also the definition (32), equation (65) can be written as:
m I Γ I V I + ρ I d V I d t d V = d F I .
The volume element ρ I d V can be factorized on both sides of (67) by defining
d F I = f I d V ,
so that all in all one has
m I Γ I V I + ρ I d V I d t = f I .
This (69) will be written in local terms by applying (26) to the fluid velocity itself, so to obtain
ρ I t V I + ρ I V I · V I = f I m I Γ I V I .
In order to enumerate the forces composing d F I we have to list those systems acting on the portion of fluid. The parcel d C I interacts with the following systems:
  • the electromagnetic fields E and B defined at its position;
  • the gravitational field g defined at its position;
  • the inertial forces measured in the reference frame co-rotating with the Earth (the SERF);
  • the remaining part of the system C I which interacts with d C I at the border d V through the nearby parcels;
  • the other systems C J with J I ;
  • the mechanical waves passing by the position of d C I .
Let us come to the complete expression of the force d F I = f I d V : as announced before, let’s assume that C I is a generic charged fluid and let’s attribute to d C I a total electric charge
d q I = δ I d V .
This charge d q I makes d C I experience a force d F I e m due to the presence of the electromagnetic field, that is written in Newtonian terms as
d F I e m = d q I E + V I × B ,
being the fields those defined in SubSection 2.5.
The second force acting on d C I is the gravitational force, represented as
d F I g r a v = m I d V g e f f :
the field g e f f is obtained in the SERF. It is basically composed by the Kepler-like gravitational potential g 0 given by the rest mass of the Earth, plus the positional part of the inertial forces due to the Earth’s motion. If the angular velocity Ω of the Earth is supposed to be constant within the time scales treated here, we can say that g e f f includes g 0 and the centrifugal force:
g e f f = g 0 + g Ω , g Ω = Ω × Ω × r .
The Coriolis force is calculated for the parcel of mass m I d V and velocity V I as:
d F I C = 2 m I d V Ω × V I = 2 ρ I Ω × V I d V .
The force applied by the rest of C I on the d C I under exam through the border d V of the small volume can be expressed as an integral over the surface d V of the outward flux of an S O 3 -tensor of rank 2, referred to as the tension tensor; this has been indicated in SubSection 2.4 as T I i j , and one can write:
d F I i = d V T I i j n j d 2 S ,
being d 2 S the infinitesimal element of the area of d V and n j the outward pointing versor locally perpendicular to d V . From Stokes theorem one can write
d V T I i j n j d 2 S = d V j T I i j d 3 x ,
that can be approximated as
d V j T I i j d 3 x j T I i j d V
as d V tends to zero. When the symbol T I is used for the tensor variable whose components are the T I i j , the force d F I may be written as
d F I = · T I d V .
In general equation (19) holds, and divides T I into a pressure part and a viscous part. When this is put into equation (76), one has
d F I = · p I 1 + Ø I d V = p I · Ø I d V .
The forces due to the other systems C J I are viscosity forces that can be seen as the friction due to the motion of the parcel d C I through the viscous medium C J I . Let the velocity of C J be V J at the position of the parcel d C I : if the time frequency of collision between an I and a J particle at r  per unit mass is some ν I J , the friction force d F J I is
d F J I ν = ν I J m I d V V I V J ;
then the total friction force due to the systems C J I reads
d F I ν = ρ I d V J I ν I J V I V J .
The last system interacting with d C I can be anymechanical wave field travelling through the atmosphere. When such a wave goes across the space it carries momentum and energy, and can exchange these quantities with whatever it meets. In particular, this may be conveniently as the divergence of a tensor · Ø I w , so the volume force applied by a mechanical wave reads:
d F I w = · Ø I w d V .
We are finally in the position of writing the force f I completely, getting rid of the infinitemal volume d V , as:
ρ I t V I + ρ I V I · V I + m I Γ I V I = δ I E + V I × B + ρ I g e f f + 2 ρ I Ω × V I p I · Ø I ρ I J I ν I J V I V J · Ø I w .
Equation (81) is the Navier-Stokes equation for theI-th fluid involved in the system, and it is clearly coupled with all the other fluids C J I , and with the electromagnetic fields. One may work with the number densities of the chemical species, putting (10) into (81):
N I t V I + N I V I · V I + Γ I V I = q I m I N I E + V I × B + N I g e f f + 2 N I Ω × V I 1 m I p I 1 m I · Ø I + · Ø I w N I J I ν I J V I V J .
Wherever the I-th element is present and N I 0 this can be furtherly written as
t V I + V I · V I + Γ I N I V I = q I m I E + V I × B + g e f f 2 Ω × V I + 1 N I m I p I 1 N I m I · Ø I + · Ø I w J I ν I J V I V J .
The sum of the two tensor terms · Ø I + · Ø I w is often referred as a unique term ϑ I = Ø I + Ø I w , so that
t V I + V I · V I + Γ I N I V I = q I m I E + V I × B + g e f f 2 Ω × V I + 1 N I m I p I 1 N I m I · ϑ I J I ν I J V I V J .
As τ I w is mainly a friction term, the ϑ I in (84) is generally referred to as general viscosity tensor.

3.4. Internal Energy Balance

The quantities as p I , N I or Ø I refer to the “inner structure” of d C I , that is composed by particles. As the parcel d C I does contain order of Avogadro number of particles of the I-th chemical, it is considered a thermodynamical system, its physical variables undergo thermodynamical relations: then p I , N I and Ø I will be involved in an internal energy balance
d U = d Q + d L .
where U is the internal energy, d Q is a portion of energy thermically exchanged and d L is a portion of work. The signs in (85) are chosen so that d Q is the heat acquired by the system, and d L is the quantity of mechanical work exerted on the parcel [25]. It is very important to underline that the energy balance we are speaking about is a thermal balance: we are not discussing about the energy that d C I releases or incorporates due to its centre-of-mass motion throughout the space; here we deal with energy exchanges that concern the internal degrees of freedom of the mechanical system d C I treated in terms of equilibrium thermodynamics of the O N A particles of the parcel (this is why equation (85) does not contain any potential energy term due to mechanical work made by external conservative forces on the centre-of-mass degrees of freedom, i.e. no gravitational or electrostatic energy term). The inner structure of d C I is mimicked again in Figure 4
Let us consider the three elements d U , d Q and d L and write them in terms of the degrees of freedom of C I . These quantities d U , d Q and d L are variations/increments happening while the parcel CoM moves along the trajectory during some time d t of a displacement d r I = V I d t . This means that d U d C I is for example the Lagrangian variation of the internal energy of the parcel d C I , so that:
d U d C I = d t d d t U d C I .
In order to write U d C I correctly we must first of all think that it is a function of volume and temperature of the gas to which it refers:
U d C I = U d V , T I .
One can write its value working on the volume-temperature plane starting at some point O = d V , T 0 with the same volume of the parcel, and integrate along an iso-volumic trajectory
U d C I = T 0 T I C V I T d T + U O d C I ,
being U O d C I the energy at the point O , and T 0 its temperature. Due to this expression, the variation d U d C I on the left hand side of (86) must be carefully calculated as
d U d C I = d t d U d C I d t = d t d d t T 0 T I C V I T d T + U O d C I .
It is possible to make this derivative when all the time dependences are carefully considered: this is done in Appendix A. All in all, the result turns out to read:
d U d C I = m I Γ I T 0 T I c V I T d T + ρ I c V I T I t T I + V I · T I d t d V .
A little argument requires a comparison between this (89) and the expression (86): in equation (86), U d C I is not a scalar field in the physical space, but a function of the parcel degrees of freedom relative to its centre-of-mass. When the expression (87) is introduced one just put in evidence this thing, making U d C I depend on r and t through the temperature T I , so that the mathematical framework to arrive to the local representation is gained. What one must keep in mind is that the thermodynamics of d C I r just depends on the local position through belonging to the parcel of CoM equal to r at time t.
Let us turn to the term d Q I of thermal exchange of energy between the parcel and the environment. Three mechanisms of heat exchange must be considered: thermal conduction d Q I , T , due to temperature differences between the d C I and the nearby parcels; mechanical dissipation d Q I , w due to interaction with the environmental part of C I (and related to the tensor Ø I ); heat internal production and dispersion d Q I , i n t due to other processes, more explicitly discussed below:
d Q I = d Q I , T + d Q I , w + d Q I , i n t .
Let us examine these addenda one by one.
The quantity of heat exchanged by the parcel due to thermal gradients can be expressed as the integral of the heat current all around the surface d V : considering that the heat spontaneously flows from warmer to cooler regions one can write
d Q I , T = d t d V κ I T I · n ^ d 2 S ,
being κ I T I the thermal current, and κ I the I-th thermal conductivity. The integral in (91) is positive when T I points from inside to outside the surface d V , so that the environment being hotter than the parcel will give a positive d Q I , T to d C I . Mathematical symbols in (91) are the same as in (76). Transforming the surface integral in (91) via Stokes theorem
d Q I , T = d t d V · κ I T I d 3 x
and considering the smallness of the material parcel, one has:
d Q I , T = · κ I T I d V d t .
The quantity d Q I , w is indicated in general as
d Q I , w = Φ I d V d t .
In a similar manner we introduce a quantity Q I so that
d Q I , i n t = ρ I Q I d V d t .
The total heat contribution is included in the equation (85) of internal energy as:
m I Γ I T 0 T I c V I T d T + ρ I c V I T I t T I + V I · T I d t d V = = · κ I T I d V d t + Φ I d V d t + ρ I Q I d V d t + d L .
Let us now examine the work term: consider that this infinitesimal work d L d C I is the portion of work made by the environment onto the thermodynamical degrees of freedom of the parcel in the time d t . In particular, d L is the expansion work made by varying the value of d V :
d L = p I d d V = p I d t d d t d V
(this expansion-compression work is different from the acceleration of the centre-of-mass of d C I : the latter would be calculated through the power of the external forces applied to the d C I in the momentum balance). Equation (30) leads to:
d L = p I · V I d V d t ,
and the internal energy equation is finally added to (81) as a proper equation of motion
m I Γ I T o T I c V I T d T + ρ I c V I T I t T I + V I · T I = = · κ I T I + Φ + ρ I Q I p I · V I .
The whole set of “fundamental” equations involving the matter degrees of freedom can be displayed as
( 34 ) , ( 81 ) , ( 97 ) ,
which give us 1 + 3 + 1 = 5 PDEs. The unknowns involved in them are more than five, anyway, let us remind that they are
ρ I , V I i , p I , τ I i j , T I
for each fluid, giving 1 + 3 + 1 + 6 + 1 = 12 quantities for each C I . Some more relationships must be added, in order to have a well posed problem.

3.5. Some More Relationships

The next relationship useful to get as many equations as the dynamical variables is the angular momentum balance for the parcel d C I , that actually gives some useful aid. Indeed, the fact that there must be angular momentum conservation in fluid dynamics leads instead to some constraints on the stress tensors, that should be symmetric for the matter systems [16]: actually, this reduces the components of τ I i j from 6 to 3. We are then left with ρ I , V I , p I , τ I i j , T I being 9 unknowns, while the system (98) contains just 1 + 3 + 1 = 5 equations.
The fundamental equations for fluids are all in (98): writing them we are finished with exploiting the general laws concerning every kind of fluid. This means that in order to obtain a number of equations equal to the number of matter unknowns, one has to start with particular hypotheses defining the special kind of fluid one is dealing with.
First of all, since we are treating the parcel d C I of the fluid just like a thermodynamical system, we can turn to thermodynamics to find other useful equations. When the hypotheses described in SubSection 2.3 for the classical continuous description hold, we can say that the system is locally near its thermodynamical equilibrium, it has a well defined state (represented by p I and T I ), and evolves by undergoing quasi-static transformations. As long as local thermodynamical equilibrium is supposed, some equation of state will hold and the tension p I , the density ρ I and the temperature T I will be constrained to fulfill some relationship of the form:
f I p I , ρ I , T I = 0
(as discussed in [1], this may not be always the case; in general the local equilibrium hypothesis is consistent with smooth dynamics). In general the ionospheric fluids and the thermosphere are assumed to be ideal gases [17]. For a volume V of n I moles of ideal gas C I one has
p I V = n I R T I ;
when the system is the parcel d C I one has V = d V and n I is the ratio between the mass of d C I and the molar mass of the I-th gas μ I :
p I d V = m d C I μ I R T I , p I d V = ρ I d V μ I R T I , p I = ρ I μ I R T I .
When the Avogadro number N A is involved the equation of state turns into the following form
p I = ρ I m I k B T I ,
being k B = 1.3811 × 10 23 Joule / mol · K the Boltzmann factor. Then (99) reads
m I p I = ρ I k B T I .
Equation (99) adds one relationship to the five equations in (98): then we have 6 conditions, 5 PDEs and 1 constraint, one the 9 variables ρ I , V I , p I , τ I i j , T I , 3 more equations are still needed.
Other important relationships rising the number of useful equations up to 9, and thus eliminating the underdetermination question of the boundary problem given by (98) plus (99), are provided by the dependence (18) giving the tensor T I in terms of the velocity V I , in particular of the non ideal part Ø I . These relationships will be of the form
τ I i j = τ I i j V I , . . .
and will be as many as the components of Ø I , i.e. 3 more conditions, as these off-diagonal terms are symmetric as τ I i j = τ I j i due to the angular momentum conservation.
As explained in [6], the components τ I i j are expected to be composed by linear functions of the space-gradients of V I : fluid viscosity tries to keep the nearby strata of the fluid in solidal motion, so the strength of its force must be (at least) linear in these infinitesimal differences; moreover, the viscous forces must be zero when the fluids move with rigid rotation
V I i = ϵ j k i Ω I j x k , i Ω I j = 0 τ I i j V I , . . . = 0 ,
because in this case no relative velocity exists between two points of the system. Thanks to this reasoning, the form suggested in [6] reads
τ I i j V I , η I , ζ I = η I j V I i + i V I j 2 3 δ i j k V I k + ζ I δ i j k V I k .
η I and ζ I are two constant scalar coefficients characterizing the nature of the fluid C I . This form of τ I is referred to as Newtonian viscosity, and it turns out to be the best fitting solution for the atmospheric fluids.
When the divergence of this τ I is evaluated, one gets
· τ I = η I 2 V I + η I 3 + ζ I · V I ,
and introducing
η I = η I 3 + ζ I
one immediately writes
· τ I = η I 2 V I + η I · V I .
As the 3 assumptions (101) are made and the τ I i j are not independent any more but rather expressed as functions of V I , the equations of motion form the following partial differential equation-plus contraint system:
t ρ I + · ρ I V I = m I Γ I , ρ I t V I + ρ I V I · V I + m I Γ I V I = δ I E + V I × B + ρ I g e f f + 2 ρ I Ω × V I p I · τ I V I , . . . ρ I J I ν I J V I V J · τ I w , m I Γ I T o T I c V I T d T + ρ I c V I t T I + ρ I c V I V I · T I = = · κ I T I + Φ I + ρ I Q I p I · V I , f I p I , ρ I , T I = 0 .
Through the assumption (105) one finally gets:
t ρ I + · ρ I V I = m I Γ I , ρ I t V I + ρ I V I · V I + m I Γ I V I = δ I E + V I × B + ρ I g e f f + 2 ρ I Ω × V I p I η 2 V I η · V I ρ I J I ν I J V I V J · τ I w , m I Γ I T o T I c V I T d T + ρ I c V I t T I + ρ I c V I V I · T I = = · κ I T I + Φ I + ρ I Q I p I · V I , f I p I , ρ I , T I = 0 .
Possible solutions, i.e. motions for the fluid, are found when one writes explicitly the quantities Γ I , Φ I and Q I in terms of the other variables, making hypotheses as (99) and (105). These assumptions will be referred to as constitutive hypotheses. The system of PDEs plus the equation of state (107) is namely all we need to describe the ionospheric gases as “smooth” fluids.
The term · τ I w in Navier-Stokes equation can be calculated when a field of mechanical waves is assigned: otherwise, one should add another dynamical variables representing the field of elastic oscillation throughout the medium, say u , and its own equation of motion, then evaluating the force · τ I w as a function of u . The elastic oscillations are usually thought of as a dissipative perturbation of the continuum that propagates throughout the system, superimposed on the motions described by the field V I ; we can say that the medium with velocity V I can be regarded as some background along which the perturbation u travels.
Note that the study of large structures of the Earth’s atmosphere as a collection of shells of matter is generally dealt with by assuming the velocity to be divergenceless, i.e. the fluid to be incompressible
· V I = 0 :
this (108) means there is no source and no sink in the flow of the C I fluid, as equation (30) says. When this incompressibility hypothesis is done, equations (107) are re-written as
t ρ I + V I · ρ I = m I Γ I , ρ I t V I + ρ I V I · V I + m I Γ I V I = δ I E + V I × B + ρ I g e f f + 2 ρ I Ω × V I p I η 2 V I ρ I J I ν I J V I V J · τ I w , m I Γ I T o T I c V I T d T + ρ I c V I t T I + ρ I c V I V I · T I = = · κ I T I + Φ I + ρ I Q I , f I p I , ρ I , T I = 0 .
It is important to underline that (108) and (109) are essentially approximations: in general ionospheric continua are compressible fluids. Nevertheless this · V I = 0 works fine for neutrals; the charged part of the upper atmosphere shows compressibility along the magnetic field B and incompressibility perpendicular to it.

3.6. Maxwell Equations

The system (107), or its incompressible equivalent (109), contains six more variables that may either be assigned as “external” rigid fields, or integrated consistently as dynamical variables: the electromagnetic fields E and B , subject to Maxwell Equations that couple them to charged matter.
The equations of motion for E and B are the well known Maxwell equations [18], that can be written in S O 3 -vector formalism as
· E = δ ε 0 , · B = 0 , × E = t B , × B = μ 0 J q + ε 0 t E .
The sources of E and B are related to C I local respective quantities. In particular, the distribution δ reads:
δ = I δ I ,
or equivalently:
δ = I q I N I , δ = I q I m I ρ I .
The vector J q is calculated as follows:
J q = I q I N I V I = def I q I J I .
These two source fields couple the (charged) ionospheric matter with the electromagnetic fields. The densities δ and J q are related to each other by a charge conservation balance that reads:
t δ + · J q = 0 .
There is no source term anyway, and this is a very severe conservation: on the being zero of t δ + · J q the whole gauge nature of the electromagnetic field does depend [26].
Equations (110) may be expressed as follows:
· E = I q I ε 0 N I , · B = 0 , × E = t B , × B = μ 0 I q I N I V I + ε 0 t E .
Maxwell Equations (115) are conveniently specialised to describe the LIM excluding extremely fast phenomena, making the following approximations:
  • the displacement current is considered much smaller than the matter currents
    ε 0 t E J q ,
    and is usually neglected:
    ε 0 t E 0 ;
  • the magnetic field produced by ionospheric currents is much smaller than the geomagnetic one
    δ B B 0 ,
    and it is assumed as zero too: this leads to the electrostatic nature of E , considering that t B 0 = 0 and that the approximation δ B 0 ends up yielding t B = 0 , so that:
    × E 0 ;
  • electric charge densities producing important electric fields are very small, and so are their time derivatives
    t δ 0 ,
    so that the corresponding currents must be thought of as divergenceless:
    · J q 0 ,
    and this is practically true almost everywhere in the ionosphere.
When the assumptions (116), (117) and (118) are considered, equations (110) become:
· E = δ ε 0 , · B = 0 , · J q = 0 , × E = 0 , × B = μ 0 J q ,
or
· E = I q I ε 0 N I , · B = 0 , · I q I N I V I = 0 , × E = 0 , × B = μ 0 I q I N I V I .
As the electric neutrality is assumed, one has a divergenceless electric field too
· E = 0 .
However there are very important exceptions to this case, as in the case of terminators and the consequent “fountain effect”, see [18,27].

3.7. Neutral Components

In this Subsection we will consider the application of equations (106) to the case of neutral constituents, those C I for which one has:
q I = 0 .
Putting (121) into (106) one finds
t ρ h + · ρ h V h = m h Γ h , ρ h t V h + ρ h V h · V h + m h Γ h V h = ρ h g e f f 2 ρ h Ω × V h + p h · τ h V h , . . . ρ h J h ν h J V h V J · τ h w , m h Γ h T O T h c V h T d T + ρ h c V h t T h + ρ h c V h V h · T h = = · κ h T h + Φ h + ρ h Q h p h · V h , f h p h , ρ h , T h = 0 :
of course,one can eliminate the electromagnetic part of the forces.
Equations (122) are usable in their own right as field equations describing the single neutral species interacting with all the other parts of the physical system. However, as one is often interested in treating the ionosphere as the charged part of the atmosphere embedded in the system of neutrals as a whole, it is useful to define the two quantities ρ and U as in (22) to describe the neutrals all together, and consider to obtain equations for these collective variables starting from (122). The equations are summed over the neutral species as
t h ρ h + · h ρ h V h = h m h Γ h , h ρ h t V h + h ρ h V h · V h + h m h Γ h V h = h ρ h g e f f 2 h ρ h Ω × V h + h p h · h τ h V h , . . . h ρ h J h ν h J V h V J h · τ h w , h m h Γ h T O T h c V h T d T + h ρ h c V h t T h + h ρ h c V h V h · T h = = · h κ h T h + h Φ h + h ρ h Q h h p h · V h .
Summing over h the equation of state would not make any sense. Instead one can add an equation of state for the whole neutral matter if suitable assumptions can be made.
Consider the mass balance equations summed: the left hand side of the first PDE in (123) can be transformed by using the definitions (22); moreover, giving the following definition
G 0 = h m h Γ h ,
one ends up with
t ρ + · ρ U = G 0 .
There are regions of the atmosphere in which the rate G 0 can be neglected with respect to the other rate of production-loss, essentially those where the presence of neutrals prevails over that of charged species. There is also a region in which the ionization is so small that it can be neglected at all: in this case one should set
G 0 0
and then (125) becomes the total mass conservation paradigm tout court, as neutrals are basically the only form of matter:
t ρ + · ρ U = 0 .
When (127) holds, the equations of motion of the whole system reduce to (122), and one has no need to refer to the equations of charged species, or to what comes out from (123) when neutrals are treated as a unique fluid component.
On the other hand, there are regions in which neutrals almost do not exist, such as the plasmasphere, where the atmospheric matter is fully ionized. In these cases one can neglect the density ρ comparing to the ρ I s of charged species. Then equation (125) becomes unimportant, the dynamical effect of neutrals becomes simply some perturbative drag affecting the motion of plasma: in such situations, the perturbative drag is simply considered as an assigned external term.
As usual, the most difficult case is that in which G 0 must be considered different from zero: if one defines some mass of the “average neutral species”
m 0 = h m h N h h N h
and also
Γ 0 = 1 m 0 h m h Γ h ,
the mass balance equation for neutrals turns out to be:
t ρ + · ρ U = m 0 Γ 0 .
Equation (130) is then put in (123) replacing the first one:
t ρ + · ρ U = m 0 Γ 0 , h ρ h t V h + h ρ h V h · V h + h m h Γ h V h = h ρ h g e f f 2 h ρ h Ω × V h + h p h · h τ h V h , . . . h ρ h J h ν h J V h V J h · τ h w , h m h Γ h T O T h c V h T d T + h ρ h c V h t T h + h ρ h c V h V h · T h = = · h κ h T h + h Φ h + h ρ h Q h h p h · V h .
Let us now consider the second equation in (131): one would like to obtain a Navier-Stokes Equation for all the neutrals as a whole, and some manipulations are necessary in order to do this, precisely those exposed in Appendix B. One finally gets an equations involving ρ and U only:
d d t ρ U = ρ g e f f 2 Ω × ρ U p · τ 0 · τ 0 w Y ν 0 Y ρ U V Y .
For applied studies the divergence · τ 0 must be written more explicitly · τ 0 = η 0 2 U , as the neutral component of the LIM is considered a Newtonian fluid (all the neutrals are), so that:
d d t ρ U = ρ g e f f 2 ρ Ω × U p η 0 2 U + η 0 · U · τ 0 w Y ν 0 Y ρ U V Y .
In order to study the large scale phenomena, it is usually assumed
· U = 0 ,
so that the equation reads:
d d t ρ U = ρ g e f f 2 ρ Ω × U p η 0 2 U · τ 0 w Y ν 0 Y ρ U V Y .
In general the incompressibility (134) is an approximation that works fine in ”quiet” conditions, i.e. when small scale mixing and turbulent phenomena are neglected. It is better said that under the turbopause, for describing the neutral matter evolution at macroscopic scales, one think of a gas of extended structures, the eddies, to which the velocity U is attributed. The viscosity η 0 is then the coefficient evaluated in the eddy interaction theory. At the same large scales, the fluid over the turbopause is treated as a perfect gas of molecules, and the factor η is the viscosity due to particle collisions.
As neutral atmosphere incompressiblity is postulated, one can re-write (131) as:
t ρ + · ρ U = m 0 Γ 0 , d d t ρ U = ρ g e f f 2 ρ Ω × U p · τ 0 · τ 0 w Y ν 0 Y ρ U V Y , h m h Γ h T O T h c V h T d T + h ρ h c V h t T h + h ρ h c V h V h · T h = = · h κ h T h + h Φ h + h ρ h Q h h p h · V h .
Let us deal with the last of the equations (136), namely thermodynamics. First of all, one has to assume that all the neutral species have the same temperature:
T h = T .
Then
h m h Γ h T O T c V h T d T + t T h ρ h c V h + T · h ρ h c V h V h = = · T h κ h + h Φ h + h ρ h Q h h p h · V h .
The extensive quantities for the neutral fluid as a whole can be obtained simply by summing the partial ones: then one has
h κ h = κ 0 , h Φ h = Φ 0 , h ρ h Q h = ρ Q 0 ,
and is able to write the equation of internal energy balance as:
h m h Γ h T O T c V h T d T + t T h ρ h c V h + T · h ρ h c V h V h = = · κ 0 T + Φ 0 + ρ Q 0 h p h · V h .
The last term on the right hand is the sum over all the contributions of mechanical work passed to thermodynamical degrees of freedom for the neutral components, so it is the total work applied to the neutrals as a whole. One has:
h p h · V h d V d t = d L 0
and since the neutrals as a whole are considered a unique fluid, still one has
d L 0 = p 0 · U d V d t ,
and then
h m h Γ h T O T c V h T d T + t T h ρ h c V h + T · h ρ h c V h V h = = · κ 0 T + Φ 0 + ρ Q 0 p 0 · U .
With similar reasonings one can also write
h m h Γ h T O T c V h T d T = m 0 Γ 0 T O T c V 0 T d T , h c V h ρ h t T = ρ c V 0 t T , h ρ h c V h V h · T = ρ c V 0 U · T .
With all these positions, one concludes:
m 0 Γ 0 T O T c V 0 T d T + ρ c V 0 t T + ρ c V 0 U · T = = · κ 0 T + Φ 0 + ρ Q 0 p 0 · U ,
so that the full system of equations for the neutrals as a whole can be read as
t ρ + · ρ U = m 0 Γ 0 , d d t ρ U = ρ g e f f 2 ρ Ω × U p · τ 0 · τ 0 w Y ν 0 Y ρ U V Y , m 0 Γ 0 T O T c V 0 T d T + ρ c V 0 t T + ρ c V 0 U · T = · κ 0 T 0 + Φ 0 + ρ Q 0 p 0 · U , f 0 ρ , p 0 , T = 0 :
clearly, in (140) an equation of state f 0 ρ , p 0 , T = 0 is assumed to exist for the neutral part of the atmosphere as a whole, and after assuming that all the T h s are equal, this sounds rather reasonable: this is just a strong mixing statement.
A last remark is deserved here, before going on: the construction of a “total neutral” equivalent fluid makes it possible to define the neutral wind reference frame (NeWReF) very useful later, namely the one with speed U point by point (see SubSection 4.1)

3.8. Charged Components

Let us finally specialise equations (107) for charged chemicals. When no approximation is made, the equations for charged fluids look like exactly the same way as (107), but realistic specializations are very useful to study the ionosphere.
The mass balance equation is exactly written as
t ρ I + · ρ I V I = m I Γ I ,
as discussed in SubSection 3.1.
The linear momentum conservation for charged species reads:
m I N I t V I + m I N I V I · V I + m I Γ I V I = q I N I E + V I × q I N I B + 2 m I N I Ω + + m I N I g e f f p I m I N I J I ν I J V I V J · τ I · τ I w .
In order to treat (142) for the ionospheric plasma some approximations are assumed.
First of all, since the charged species are very dilute in all the parts of the atmosphere, we can neglect both the self-viscosity term of C I and the transmission of the momentum from the mechanical waves to C I :
q I 0 · τ I 0 , · τ I w 0 .
This approximation is justified in different ways in the different regions of the upper atmosphere. In its lower region, where one can speak about ionosphere properly, the charged matter is extremely less abundant than the neutral one: for q I 0 the fluids C I are very dilute gases merged in the much denser medium C 0 . This renders the forces “applied by the charged matter on itself” extremely weaker than those applied by C 0 , and negligible at all. When the ionospheric-plasmaspheric limit is passed, neutrals are not present any more in practice, then the force density · τ I + · τ I w for the plasma cannot be negligible with respect to those other terms, but it becomes negligible by itself due to the diluteness of the very plasma at those heights, where densities are extremely low. One can say that (143) holds everywhere, and equation (142) becomes:
m I N I t V I + m I N I V I · V I + m I Γ I V I = q I N I E + V I × q I N I B + 2 m I N I Ω + + m I N I g e f f p I m I N I J I ν I J V I V J .
There are other negligible terms.
The velocity-skew-linear force
V I × q I N I B + 2 m I N I Ω
can be approximated: due to the very small mass of the particles forming the plasma, and to their diluteness, one can write
q I N I B 2 m I N I Ω ,
and when this is used, the velocity force can be approximated as
V I × q I N I B + 2 m I N I Ω V I × q I N I B ,
so that Coriolis forces disappear from (144): after a little bit of algebra one has:
t V I + V I · V I + Γ I N I V I = g e f f + q I m I E + q I m I V I × B 1 m I N I p I J I ν I J V I V J .
Equation (146) can be specialized for each charged species. If the chemical reactions are neglected in the Navier-Stokes PDE
Γ I N I V I o
equation (146) reduces to
d V I d t = g e f f + q I m I E + q I m I V I × B 1 m I N I p I J I ν I J V I V J .
Otherwise one must re-start from equation (146) and write
d d t N I V I = q I m I N I E + q I m I N I V I × B + N I g e f f 1 m I p I N I J I ν I J V I V J
(the production-loss term is, if any, contained in the Lagrangian derivative d d t N I V I , involving d N I d t ). An example of how (147) holds true is given in Appendix C.
Equation (148) or (149) can be adapted to any charged species, however in this Section we will use a particularly simplified atmospheric picture, in which only two charged species are present: positive 1-valent ions, with charge + e and mass M, and electrons, with charge e and mass m. In this picture negative ions and multi-valent ions are completely neglected. Moreover, the approximation (147) will be considered as true, so we will start from equation (148). This model of population C I = X , X + , e can be applied, e.g., to the topside of the ionosphere and to the plasmasphere (being X + different chemicals for the different regions). For negative ions one should refer to what is reported for example in [19]. Note also that physical models of the free electron density N e are influenced seriously by the presence of negative ions, as shown in [28].

3.8.1. Ion Motions

The specialization of (148) to the e , M positive ions reads:
d V i d t = g e f f + e M E + e M V i × B 1 M N i p i ν i e V i V e ν i n V i U ,
being U the bulk velocity of neutral matter, and V e that of electrons. These positive ions are considered as a perfect gas in local thermodynamical equilibrium:
p i = N i k B T i .
Last but not least, ions undergo mass balance:
t N i + · N i V i = Γ i .
The equations of motion for the positive ionic species then read:
t N i + · N i V i = Γ i , p i = N i k B T i , d V i d t = g e f f + e M E + e M V i × B 1 M N i p i ν i e V i V e ν i n V i U .
In many cases the electron-ion collision frequency ν i e is negligible with respect to neutral-ion collision frequency ν i n , and so equations (153) read
t N i + · N i V i = Γ i , p i = N i k B T i , d V i d t = g e f f + e M E + e M V i × B 1 M N i p i ν i n V i U :
this is the case of the low ionosphere, where the probability of meeting neutrals is much higher than that of meeting electrons.
On the opposite side, in the high ionosphere, neutrals are almost absent, so that one has ν i n = 0 and (153) must be re-written as:
t N i + · N i V i = Γ i , p i = N i k B T i , d V i d t = g e f f + e M E + e M V i × B 1 M N i p i ν i e V i V e .

3.8.2. Electron Motions

Let us specialize equation (148) to electrons, identified by e , m , so the version of equation (148) for them reads
d V e d t = g e f f e m E e m V e × B 1 m N e p e ν e n V e U ν e i V e U .
If the electrons are considered as a perfect gas without quantum effects, the law
p e = N e k B T e
is added; for what concerns the mass balance, one can write
t N e + · N e V e = Γ e .
Free electrons due to ionization are then described by the set of equations
t N e + · N e V e = Γ e , p e = N e k B T e , d V e d t = g e f f e m E e m V e × B 1 m N e p e ν e n V e U ν e i V e U .
When the ion-electron collisions are neglected with respect to neutral-electron ones
ν e i ν e n
it is possible to write:
t N e + · N e V e = Γ e , p e = N e k B T e , d V e d t = g e f f e m E e m V e × B 1 m N e p e ν e n V e U .
Approximation (160) holds in low-middle ionosphere, while in the highest layers one has ν e n = 0 as no neutrals can be found there:
t N e + · N e V e = Γ e , p e = N e k B T e , d V e d t = g e f f e m E e m V e × B 1 m N e p e ν e n V e U ν e i V e U .
The two rates Γ i and Γ e appearing in (153) and (159) respectively are not independent as stated in SubSection 3.1, and equation (48) holds
e Γ i e Γ e = 0 ,
so that as the statement Γ i = Γ e Γ is made, one has:
t N i + · N i V i = Γ , t N e + · N e V e = Γ .
The rate Γ represents the rate of the events X X + + e , each of which will “create” a positive ion and a free electron, so it is no surprise that \Gamma represents both the e as well as the X + production-loss rate.

4. The Quiet Ionosphere

The presence of plasma in the upper atmosphere renders the matter in this region a gaseuos conductor, because free charges can make macroscopic motions under the action of electromagnetic fields. Obtaining a description of the system as a conductor at the equilibrium, with the geomagnetic field impressed, requires the medium to be in local mechanical equilibrium, namely no bulk acceleration for the continua C I :
d V I d t = 0 .
What one wants to find is a relationship between the charged species motions and the electromagnetic fields present under the static condition (165): in this condition one speaks about quiet ionosphere.
When the conditions (165) is put into (84), the momentum balance turns into a relationship between the V I s and the electromagnetic fields:
g e f f + q I m I E + q I m I V I × B 1 m I N I p I J I ν J I V I V J ν n I V I U = 0
(note that, to get (166) from (84) and (165), some approximations were considered, namely those discussed in Section 3). The two assumptions t V I 0 and V I · V I 0 , implying d V I d t 0 , are both supported, in practice, by the experimental reality of the ionosphere in “geomagnetically quiet conditions”. For this, let us refer as τ to the rensponse time of the medium to the applied forces, and as L to the characteristic length of the velocity field space variability. Then one can say
t V I = O V τ , V I · V I = O V 2 L .
The other terms of the acceleration in (84) involving the velocity satisfy:
J I ν J I V I V J = O V Δ t , q I m I V I × B = O V τ I B ,
where Δ t is the inverse of the collisional frequency and τ I B is the gyration period of the I-th species around the magnetic field line. In the experimental reality one always has τ Δ t , τ I B so that one may set t V I 0 . Moreover, in the quiet atmosphere the length L is usually much larger than gyration radii and free mean paths (length scales associated to Lorentz force and collisional friction respectively), so that also the approximation V I · V I 0 holds. From the just supported (166), one has:
m I N I g e f f + q I N I E + V I × B p I m I N I J I ν J I V I V J m I N I ν n I V I U = 0 ,
and here the equation of state (99) can be used to write p I in terms of N I and the temperature T I
m I N I g e f f + q I N I E + V I × B N I k B T I + m I N I J I ν J I V I V J m I N I ν n I V I U = 0 .
In general the medium is treated as electrically neural point by point, so that the densities satisfy the constraint
I q I N I = 0 .
It is not pleonastic to remember that (169) is held true only around those places where
I q I · V I 0
due to the continuity equation. This excludes places as the neighborhood of the terminators, as already noticed. Due to (169), and because we have assumed the ionosphere to be made out of X, X + and e only, one may safely write:
N i = N e = def N .
Equation (168) can be used to determine the velocities V I as functions of the other quantities. This will depict the stream portrait of the continuum C I when it is under the mechanical equilibrium (165) point by point.
The Navier-Stokes equation for the system of the two charged species “i” and “e” at the mechanical equilibrium read:
g e f f + e M E + e M V i × B 1 M N p i + ν e i V i V e ν n i V i U = 0 , g e f f e m E e m V e × B 1 m N p e + ν e i V e V i ν n e V e U = 0 .
Where one can assume that ion-electron collisions are rare
ν i e 0 ,
the equilibrium momentum balance equations read:
m N g e f f e N E + V e × B k B N T e m N ν e n V e U = 0 , M N g e f f + e N E + V i × B k B N T i M N ν i n V i U = 0 .
On the contrary, if the ion-collision frequency is the only relevant ν I J
ν n e 0 ,
equations (172) will read:
m N g e f f e N E + V e × B k B N T e m N ν e i V e V i = 0 , M N g e f f + e N E + V i × B k B N T i M N ν e i V i V e = 0 .

4.1. Thermally Homogeneity and the NeWReF

These systems (174) and (176) can be resticted to the case of thermally homogeneous regions
T I = 0 .
This (177) should not be misunderstood: the ionosphere has never null temperature gradients, here one simply means that we are going to treat small enough regions so that the same temperature T I exists all over it. When this assumption holds, it is possible to go on as follows for (174):
m N g e f f e N E + V e × B k B T e N m N ν e n V e U = 0 , M N g e f f + e N E + V i × B k B T i N M N ν i n V i U = 0 ,
while for (176) one has:
m N g e f f e N E + V e × B k B T e N m N ν e i V e V i = 0 , M N g e f f + e N E + V i × B k B T i N M N ν e i V i V e = 0 .
The electrology of the quiet ionosphere as a conductor is essentially that of a conducting object dragged through the geomagentic field by the motions of neutral matter. In order to simplify our reasoning, it is convenient to work in the NeWReF, because this simplifies the dynamics observed. This NeWReF sees the ionosphere more similar to some conductor studied in its reference frame, in which traditional electrological quantities, such as resistivity, are easily defined and applied. Equations (178) and (179) are written as the SERF observer sees them, with g e f f containing the inertial position forces, negligible Coriolis acceleration and neutral wind velocity equal to U . One can also describe the same system as seen by a reference frame that has, point by point and time by time, the speed U of the neutral wind with respect to the SERF. The fields E and B are relativistic quantities, and will be measured as
E = 1 1 U 2 c 2 E + U × B , B = 1 1 U 2 c 2 B 1 c 2 U × E
in the NeWReF. In the non relativistic wind limit
U c 1 1 ,
which is definitely met, one can neglect the U 2 c 2 effects putting
1 U 2 c 2 1
and write:
E E + U × B , B B 1 c 2 U × E .
The Earth’s ionosphere is a specific situation in which
B 1 c 2 U × E
so that the transformations giving the electromagnetic fields E , B measured in the neutral wind frame in terms of those E , B measured by the Earth-solidal observer, rather read:
E = E + U × B , B = B .
Note that only the electric fields are changed. Note also that (182) simply means that the Lorentz currents due to the neutral drag produce magnetic inductions much smaller than the magnetic field measured in the SERF.
As observed in [18], the relations (183) seem to indicate that the sources of the electric fields change significantly when a “low-relativistic” transformation is applied, while those of the magnetic field do not. This is true, it can be seen rigorously in the case of a well conducting medium, as explained in the following. Suppose to make a Lorentz transformation from some system S and a second system translating rigidly with respect to it, with speed
U = U z ^ :
the relativistic theory say that the collection of four quantities c δ , J transform as a 4-vector, i.e.
J x = J x , J y = J y , J z = 1 1 β 2 J z β c δ , δ = 1 1 β 2 δ β c J z ,
being β = U c . If we make a Taylor expansion up to O β we get
J x = J x , J y = J y , J z = J z c β δ + . . . , δ = δ 1 c β J z + . . .
so that, if Δ f is the difference between the values in S and in S of a quantity f, one has
Δ J z = c β δ , Δ δ = 1 c β J z .
The relative variations are
Δ J z J z = c β δ J z , Δ δ δ = β J z c δ ,
and if we evaluate their ratio we get:
Δ δ δ Δ J z J z 1 = J z 2 c 2 δ 2 .
From (186), we get that if a very little charge density produces an intense current J z , via intense electric fields in an importantly conducting medium, which is the ionosphere, one has
Δ δ δ Δ J z J z 1 1 ,
that means a dependence of the electric density on the reference frame much higher than the current dependence.
Passing from the SERF to the NeWReF, the matter velocities are transformed as
v = v U
in the non relativistic limit, so that one has to put the expressions
V e = V e + U , V i = V i + U , E = E U × B , B = B
in equation (178), getting:
m N g e f f e N E + V e × B k B T e N m N ν e n V e = 0 , M N g e f f + e N E + V i × B k B T i N M N ν i n V i = 0 ,
or in equation (179), obtaining:
m N g e f f e N E + V e × B k B T e N m N ν e i V e = 0 , M N g e f f + e N E + V i × B k B T i N M N ν e i V i = 0 .
Dividing by N m I ν I and collecting the terms involving the velocity in (189) and in (190), one has
V e + e m ν e n V e × B = 1 ν e n g e f f e m ν e n E k B T e m ν e n ln N , V i e M ν i n V i × B = 1 ν i n g e f f + e M ν i n E k B T i M ν i n ln N
and
V e + e m ν e i V e × B = 1 ν e i g e f f e m ν e i E k B T e m ν e i ln N , V i e M ν e i V i × B = 1 ν e i g e f f + e M ν e i E k B T i M ν e i ln N
respectively. It is important to argument about the fact that no new forces appear in the NeWReF with respect to those in the SERF, and this is in principle wrong. Yet, we are considering the condition in which all the species C I have null acceleration, so that one has d U d t 0 too implying no new linear acceleration forces in NeWReF with respect to what seen in SERF. Moreover, we will consider the x , y , z axes of the NeWReF time by time and point by point parallel to the axes in the SERF, so that no new rotational forces may appear.
At this point, it is useful to define the following coefficients
K I = q I B m I ν I , D I = k B T I m I ν I , b I = q I m I ν I , H I = k B T I m I g e f f ,
because through them equations (191) and (192) are unified as:
V e K e V e × b ^ = D e H e g ^ e f f + b e E D e ln N , V i K i V i × b ^ = D i H i g ^ e f f + b i E D i ln N , b ^ = vers B = vers B , g ^ e f f = vers g e f f .
The ultimate use of defining the K I s is to use them as conductivity parameters: indeed, this K I reads
K I = Ω I ν I ,
being Ω I = q I B m I the gyration frequency of the C I particles within the field B , and ν I is the most relevant collisional frequency affecting their motion. All in all, equation (195) states that the K I measures the relevance of magnetic forces relative to the dissipation ones in the charged particle dynamics, being de facto a conductivity measurement. Then, the K I 0 limit represents the isolating limit, while the hyperconductive limit is K I + .

4.2. An Anisotropic Conductor: role of B

We have already stated that electric charges do move much more easily along the magnetic field B than perpendicular to it. It is then expected that the conductance features of the LIM will be different along B and across it. What is done usually it is to subdivide the local 3D-space of the ionospheric medium into some R along B and some R 2 perpendiculat to it. In Cartesian coordinates of local S O 3 -vectors one defines the B -parallel and B -perpendicular projectors
P i j / / B = B i B j B 2 , P i j B = δ i j B i B j B 2
(these projectors are the same when written in the NeWReF and in the SERF, as B = B ).
The analysis starts using some “generalized velocities”
W e = V e K e V e × b ^ , W i = V i K i V i × b ^ , b ^ = B B ,
that coincide with the real velocities when the magnetic field is suppressed. These W read:
W e = D e H e g ^ e f f + b e E D e ln N , W i = D i H i g ^ e f f + b i E D i ln N .
We will treat the two of (198) as one equation, omitting the species index and putting V I = v and K I = K for simplicity. One has
v K v × B = 1 ν g e f f + b E k B T n M ν n W F
(being the vector W F the W written as a function of the forced applied). Working with the local Cartesian components one has
v j K ε j h k v h B k = W j F .
In order to work with the transverse and longitudinal components we can apply to (199) the two projectors defined in (196): doing so, one performs all the calculations reported in Appendix D, and finally obtains:
v = A 1 K · W F ,
being the matrix A 1 K the one reported in (A20). Note the two limits, easily interpreted via the equation (195):
lim K 0 A 1 K = 1 0 0 0 1 0 0 0 1 , lim K + A 1 K = 0 0 0 0 0 0 0 0 1 .
Since W F lays on the plane x O y , orthogonal to B = B k ^ in the Appendix D, one gets:
W F = w x w y 0
(see also Figure 5), and applies the matrix A 1 K in (A20) to this vector, in order to obtain
A 1 K · W F = 1 1 + K 2 w x + K 1 + K 2 w y K 1 + K 2 w x + 1 1 + K 2 w y 0 .
Let us turn again to (201). In the K + limit, the tensor A 1 K projects directly along the magnetic axis: when the conductivity becomes very large, the component v tends to become the component of W along B which is clearly zero. So, this model does not allow for motions orthogonal to B forKexactly infinite. When the ionosphere is “tremendously conducting”, it allows only for motion along B of the charged species.
The vector A 1 K · W F is still orthogonal to B , and can be decomposed into the two directions W F and W F × B ^ spanning the cross- B plane. The direction W F × B ^ is
W F × B ^ = w y w x 0 ,
so the decomposition reads:
V I = 1 1 + K I 2 W I + K I 1 + K I 2 W I × b ^ .
The full solution giving the velocity of the I-th species in the stationary case reads:
V I = W I / / + 1 1 + K I 2 W I + K I 1 + K I 2 W I × b ^ .
This is expanded as:
V I = D I H I g ^ e f f / / + b I E / / D I N N / / + + g e f f ν I 1 + K I 2 + b I 1 + K I 2 E D I 1 + K I 2 N N + + K I 1 + K I 2 1 ν I B g e f f × B K I 1 + K I 2 D I B N N × B + K I 1 + K I 2 b I B E × B .
In order to understand the limit K I 1 in which the gyration prevails over the inter-particle collisions, it is convenient to expand la Taylor the coefficients
1 1 + K I 2 , K I 1 + K I 2
with respect to K I 1 , getting
1 1 + K I 2 = K I 2 K I 4 + O K I 6 K I 1 + K I 2 = K I 1 K I 3 + O K I 5 :
if one has K I 1 it is possible to truncate those expressions as
1 1 + K I 2 = K I 2 + . . . , K I 1 + K I 2 = K I 1 + . . .
When K I 1 tends to zero the surviving terms are:
lim K I 1 V I = D I H I g ^ e f f / / + b I E / / D I N N / / + + 1 ν I K I B g e f f × B D I B K I N N × B + b I B K I E × B .
Considering the definitions (193) one has:
lim K I 1 V I = 1 B 2 E × B k B T I q I B 2 N N × B + m I q I B 2 g e f f × B + + q I m I ν I E / / + 1 ν I g e f f / / k B T I N m I ν I N / / .
In the limit K I 1 the terms depending on the species at hand are
q I m I ν I E / / + 1 ν I g e f f / / + m I q I B 2 g e f f × B + k B T I N m I ν I N / / k B T I q I B 2 N N × B ,
while the term
V E B = 1 B 2 E × B
is the same for every charged species, and is referred to as E-cross-B drift. When gravitational and concentration gradient effects are negligible, the part of V I orthogonal to B is reduced to the E-cross-B drift only for all the charged species. Due to this condition, it is impossible to have any net current perpendicular to the magnetic field.
In the non conducting limit K I 0 the inter-particle collisions prevail on the gyration
lim K I 0 V I = W I / / + W I = W I F :
in this case one experiences a viscous motion for the charged species, with the velocity directed as the forces W i F . In this regimenegative and positive particles have different motions because they undergo different electromagnetic forces.
We are finally in the position of treating the ionospheric medium as an anisotropic conductor.

4.3. Ohm’s Law for the Ionosphere

All the work done down to here has been useful to describe the local ionospheric equilibrium as a conductor, and give its Ohm’s law.
If the local electric field is E we will look for the expression
J q = J q E
relating the electric current to the given field. In general the relation (212) is a smooth dependence of any kind
J k E = J k 0 + n = 1 + n J k E i 1 . . . E i n 0 E i 1 . . . E i n
(the subscript “q” has been omitted for simplicity in the local Cartesian components J h of J q ; the S O 3 -tensor summation convention on indices is again in use) and in our case we will make a linear medium assumption:
J k E = J k 0 + J k E i 0 E i .
Equation (214) may be re-written defining the conducibility tensor
σ k i = J k E i 0 ,
so that one has
J q E = J q 0 + σ · E ;
using fussy primes in the NeWReF:
J q E = J q 0 + σ · E .
σ will be constructed in the NeWReF directly, then it will be transformed to the SERF.
If the medium is “the simplest dirty plasma possible”
C I = X + , e , X ,
then the electric current reads:
J q = N e V i V e .
This vector is calcualted explicitly in all its three components J x , J y , J z in the NeWReF through the calculations in Appendix E. Only some part of the electric currents is determined by the electric field E , since the concentration gradient and gravity force contribute as well. Searching for the form (217) for the part of J determined by the electric field, one has
J x E = N e b i 1 + K i 2 N e b e 1 + K e 2 E x + N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E y , J y E = N e b i 1 + K i 2 N e b e 1 + K e 2 E y N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E x , J z E = N e b i b e E z
(consider the definitions (193)). Then the non vanishing components of σ read:
σ x x = σ y y = N e b i 1 + K i 2 N e b e 1 + K e 2 σ P , σ x y = σ y x = N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 σ H , σ z z = N e b i b e σ 0 .
σ can be written as a matrix:
σ = N e b i 1 + K i 2 N e b e 1 + K e 2 N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 0 N e K e b e 1 + K e 2 N e K i b i 1 + K i 2 N e b i 1 + K i 2 N e b e 1 + K e 2 0 0 0 N e b i b e .
The part J q 0 of the electric current, which is determined by non-electric factors, can be written by introducing two more tensors α and β . If the tensor
α = N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e 0 N e K e 1 + K e 2 1 ν e N e K i 1 + K i 2 1 ν i N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e 0 0 0 N e ν i N e ν e
and also
β = e D e 1 + K e 2 e D i 1 + K i 2 e K e D e 1 + K e 2 e K i D i 1 + K i 2 0 e K i D i 1 + K i 2 e K e D e 1 + K e 2 e D e 1 + K e 2 e D i 1 + K i 2 0 0 0 e D i D e
are defined, the linear relation between J q and its causes takes the form:
J q = σ · E + α · g e f f + β · N .
The electric addendum
J q E = def J q J q 0 = σ · E
can be put into the form
J q E = σ P E σ H E × B ^ + σ 0 E / / , B = B k ^ , E / / = E z k ^ , E = E x i ^ + E y j ^ .
Let us give an order of magnitude for the various contributions to J q . Fixing the height as:
h = 300 km N e = 20 × 10 11 m 3
for a realistic equatorial daytime condition, we have:
σ P , σ H σ 0 , σ 0 39.8 Ohm / m .
With an electric field E of
E 50 mV / m
one obtains:
σ 0 E 1.99 × 10 3 A / m 2 .
According to the values of the involved quantities as given in [18] and [28], the elements of the matrix α are:
α x x = α y y 4 × 10 11 As 2 / m 3 , α x y 4 × 10 9 As 2 / m 3 , α z z 4 × 10 7 As 2 / m 3 ,
so that
α x x g e f f 3.9 × 10 10 A / m 2 , α x y g e f f 4 × 10 8 A / m 2 , α z z g e f f 4 × 10 6 A / m 2 .
The matrix β has the following values instead
β x x = β y y 1.9 × 10 16 Am 2 , β x y 1.6 × 10 18 Am 2 , β z z 0
and since we assume
Δ = 20 km Δ N e = 5 × 10 11 m 3
we can say:
β x x N e 4.75 × 10 9 A / m 2 , β x y N e 4 × 10 11 A / m 2 , β z z N e 0 .
These values, even if roughly calculated, show that the non-electric contributions to J q are negligible with respect to the electric ones.
Equation (223) can be transformed to the SERF, where it becomes
J q = σ · E + U × B + α · g e f f + β · N ,
since J q = J , E = E + U × B and the S O 3 -tensors stay unchanged, because they do not rotate, and the field on which they depend is the magnetic field B , that remains unchanged due to equation (183).
Let’s turn to (224): in the SERF we can write
J q = σ P E σ H E × B ^ + σ 0 E / / + σ · U × B + α · g e f f + β · N ,
so that the part of J q of electric origin reads
J q E = σ P E σ H E × B ^ + σ 0 E / / .
The term
σ H E × B ^ = N e 2 b M ν + e 2 B 2 M ν N e 2 b M ν + e 2 B 2 M ν E × B J E B
is considered asthe part of the electric current due to the E-cross-B drift. This J E B  tends to zero when conducibilitiesKdiverge: in fact the expression
J E B = N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E × B ^
can be Taylor-expanded in K I 1 up to the linear terms. In this way we get
N e K i b i 1 + K i 2 = N e b i K i N e b i K i 3 + O K i 5 = N e b i K i + . . . = N e B + . . . , N e K e b e 1 + K e 2 = N e b e K e N e b e K e 3 + O K e 5 = N e b e K e + . . . = N e B + . . . ,
so that for negligible K i 1 and K e 1 the factor N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 is the difference of two equal quantities, i.e. it vanishes. The following limit holds:
lim o K 0 J E B = 0 .
This was already highlighted at page Section 4.2.
Equation (225) of the upper atmosphere as a conducting medium can be used “directly”, i.e. to calculate the currents observed in presence of a given electric field; it can also be used “inversely”, obtaining the fields E necessary to produce the given currents. The latter is the most important use of equations (225) and (223), where the fields necessary to the observed currents are deduced, so that one tries to understand their possible geophysical/space origin.

4.4. Ionospheric Electrostatic Fields

Ionospheric electric fields do exist because magnetic, gravity and pressure forces produce (even slightly) different effects on ions and electrons: when ions and electrons move differently and get separated, the local neutrality of the medium can be violated, so some δ 0 arises at some position. The local charge separation does produce an electric field so that
· E = δ ε 0 .
As these δ s may vary with time, they give rise to a local current too, since · J q = t δ , while δ passes from zero to a finite value, or viceversa. In the stationary condition described here, the current J q is a divergenceless vector: anyway this is a dynamical equilibrium, since the very complicated forces W F yield some J q that is in general · J q 0 . This will in turn produce a field E re-assessing charges and currents, up to the moment in which t δ = 0 holds again, and then · J q = 0 .
In the case one has some electric charge [18], for instance
δ 9 · 10 17 C m 3 ,
the electric field is constant over a distance of L = 1 km , with a value
E 10 mV m 1 :
in order to see how long · J q 0 will last, one considers this unknown time τ and write:
δ τ = E L σ τ = L δ σ E ,
and sees that the solenoidal nature of J q is restored within some microsecond:
τ 10 6 s .
Processes in which neutral winds (those of speed U in the foregoing pages) produce ionospheric motions that induce electric charge separations, and currents and electric fields, are referred to asdynamo processes. The power spent by the wind dynamo on the system is positive
w wind > 0 ,
as the wind drags actively the ionospheric medium. In the steady regime the atmospheric winds spend some power w wind on the plasma motion balancing exactly the power w E
w E = E · J q
exerted by E on J q :
w wind + w E = 0 w wind = E · J q .
Neutral winds do move the ionosphere across the geomagnetic field, so their power is positive, and this leads to the negativity of w E
E · J q < 0 ,
accordingly to the fact that this E arises in opposition to the motion separating the charges. This dynamo basically works in the same way of working of a battery closed on a conductor: in this way, the electrostatic E due to the wind-induced separation of charges tries to make the charges fall back on their opposite sign fellows, driving δ to zero, as the field of a short circuit nullifying the work of the battery.
In the upper atmosphere there exists also an opposite mechanism in which some external electric field does some positive work on the electric charges of the plasma, i.e. on the currents. These currents in turn spend power on the neutral matter, this time making negative work on them: this is the mechanism of the auroral currents induced by the electrodynamic field, in turn due to the precipitation of the charged Solar Wind. In this case the external field due to the Solar Wind gives some positive power to the currents
w E = E ext · J q > 0 ,
while the current transmission through the thermosphere dissipates power via Joule effect, heating the atmosphere. This mechanism is referred to as ion drag because inside these current regions the ions are the ones that mainly interact with neutrals, due to inertial reasons (see also the approximation (173) in which the electron drag is neglected).

4.5. Electric Field “Propagation”

Due to the anisotropy of the tensor σ , the upper atmosphere can be electrically represented as a conducting cage formed by the geomagnetic field lines along which charges move as freely as along a metallic wire, and across which they cannot move, precisely because σ 0 σ H , σ P . For example, at mid and high latitudes one has
h 3000 km σ P σ 0 10 6 .
Using the definitions (219) and (193), one recognises
σ P σ 0 = M m ν i ν e M m ν i ν e + e 2 B 2 m 2 ν e 2 + e 2 B 2 M 2 ν i 2 + e 2 B 2 = O ν 2
and
σ H σ 0 = e 2 B M m ν i ν e M ν i m ν e m 2 ν e 2 + e 2 B 2 M 2 ν i 2 + e 2 B 2 = O ν 4 ,
and sees clearly that where collision frequencies tend to vanish, those rations do, getting the approximation (233) (see below). This reduces the tensor σ simply to the tensor equal to σ 0 along B and vanishing in the directions perpendicular to it.
In the roughest approximation, for σ 0 , we can think of the field lines of B as of perfectly conducting lines, which are equipotential, and with E that is orthogonal to them point by point, see Figure 6.
In this subsection we want to discuss two examples of how the geomagnetic conducting cage affects the configuration of E r , t . On the one hand, we want to discuss how perturbations δ E propagate along the B -lines; on the other hand, how those perturbations δ E may be produced by the injection of B -aligned currents at high latitudes.
All the calculations are done in the quiet ionosphere approximation, i.e. considering the B -lines as conducting wires in electrostatic equilibrium.

4.5.1. Pole to Equator δ E “Propagation”

Let us consider the equations describing E in the electrostatic regime. Within the approximation
σ H , P σ 0 0
one can see that if some electric field is produced in a point of the B line, this is expected to propagate all along the line, perpendicular to it point by point, in the same order of time in which charges displaced on a conductor weakly out of the equilibrium re-displace and reach the equilibrium again. This is a possible mechanism for the propagation of electric fields from high latitudes at low height to low latitudes at high height, as observed in some electric field correlation study [29], in which some equatorial δ E takes place almost simultaneously to some polar δ E with the same time-behaviour. This is heuristically understandable by looking at the geomagnetic dipole profile, represented in Figure 7: if some field E is produced at the polar oval at some point P of height h (due to the injection of charged particles by the Solar Wind, for example, that causes locally δ 0 ), this will be transported as a vector orthogonal to the geomagnetic line up to equatorial latitudes to a point P of height h , that can be much higher than h, due to the shape of the dipolar geomagnetic line. This is grossly how the precipitation of charged matter at low heights in the polar zone can yield a variation δ E of the ionospheric electric field at low latitudes, at much higher quotes.
All this reasoning has to be put into a quantitative frame.
When the stationary condition · J q = 0 holds, the PDE for E of the ionosphere is written by using (225)
· σ · E = · σ · U × B · α · g e f f · β · N ;
then, going to the NeWReF, we have:
· σ · E = · α · g e f f · β · N
(we are omitting species indices for simplicity). If “non electric” contributions to the electric current are negligible α · g e f f + β · N 0 , the equation for E reads:
E j i σ i j + σ i j i E j = 0 .
Equations (235) state that the source of E can be a discontinuity in the conducibility of the medium, as can be seen by writing it as
σ i j i E j = E j i σ i j ,
as it happens, for example, at the terminators. The stationary condition · J q = 0 is a quasi-equilibrium condition, so E is reasonably “quasi-electrostatic”, and we look for solutions of (235) in the form of gradients of a scalar function ϕ (the electrostatic potential):
ϕ / E = ϕ .
The electrostatic potential is related to J q as
· J q = j ϕ i σ i j σ i j i j ϕ .
The PDE for ϕ reads:
j ϕ i σ i j + σ i j i j ϕ = 0 .
For quasi-homogeneous media
j ϕ i σ i j σ i j i j ϕ ,
meaning that we are working in regions where σ is slowly dependent on space coordinates, the solenoidal current condition reads
σ i j i j ϕ = 0 .
This has some sense, as one sees, for instance, looking at Figure 2.6 of [18] one can redact this table
σ P h 4 × 10 12 ÷ 6 × 10 11 Ohm / m 2 , σ H h 9.9 × 10 11 Ohm / m 2 , σ 0 h 1.13 × 10 4 ÷ 1.2 × 10 4 Ohm / m 2 .
Clearly, conditions (240) excludes σ -discontinuity regions as the terminators.
The more general PDE (239) is processed with some algebraic manipulation, thoroughly illustrated in Appendix F. Its solution is found in terms of ϕ x , y , z , being z the local magnetic direction and x and y the coordinates perpendicular to it (see Figure 5). An assumption is made about how σ 0 and σ P depend on z (see equation (A43)), and the following solution is given, in terms of a Fourier series:
ϕ x , y , z = = e 1 2 c 0 σ P z σ 0 z z k f 0 k e i k x x + i k y y A + k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z + + A k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z .
This solution of the anisotropic Poisson equation (239) gives some indications about the propagation of the electrostatic disturbances along the geomagnetic lines. Of course, the term “propagation” has to be taken with some care, as (239) is an equilibrium, time-independent PDE: it will not describe how electric waves travel, rather how electrostatic fields re-arrange in case of perturbation of an equilibrium, to reach another equilibrium. These solutions are characterized by the value of k 2 : this k is formally the wave vector of the potential along the directions orthogonal to the geomagnetic lines, and the inverse of its strength = 2 π k gives the characteristic (wave)length of the potential perpendicular to B . The k -th addendum in (242) varies orthogonally to B on scales of this order: the inverse of the strength of k gives roughly the transverse extension of the electrostatic potential structures around the geomagnetic line.
Assume to induce the variation
ϕ ϕ + δ ϕ
in the electrostatic potential all over the B -line, passing from a stationary condition to another one, still stationary: in such case the field δ ϕ must solve the electrostatic equation too, because that PDE is linear, and then must have the form (242)
δ ϕ x , y , z = = e 1 2 c 0 σ P z σ 0 z z k F 0 k e i k x x + i k y y B + k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z + + B k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z ,
just as ϕ and ϕ = ϕ + δ ϕ do. In particular one can write:
δ ϕ x , y , 0 = k F 0 k e i k x x + i k y y B + k 2 + B k 2 .
This δ ϕ x , y , 0 is the potential variation produced at some z = 0 , that we consider to be in the Polar region; on the other hand, the difference between δ ϕ x , y , z and δ ϕ x , y , 0
δ ϕ x , y , z δ ϕ x , y , 0 = = k F 0 k e i k x x + i k y y e 1 2 c 0 σ P z σ 0 z z B + k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z + + B k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z B + k 2 B k 2
can be expanded in power series as
δ ϕ k x , y , z δ ϕ k x , y , 0 = = F 0 k e i k x x + i k y y 1 2 c 0 2 + 16 π 2 2 B + k 2 B k 2 σ P z σ 0 z z + 1 2 c 0 B + k 2 + B k 2 σ P z σ 0 z z + O σ P z σ 0 z z 2 :
since the term σ P σ 0 is very small in the upper atmosphere, one can say that the huge part of δ ϕ x , y , z δ ϕ x , y , 0 is already included in (245). Now, it is instructive to note that
lim 0 δ ϕ x , y , z δ ϕ x , y , 0 lim + δ ϕ x , y , z δ ϕ x , y , 0 :
those perturbative structures of the electrostatic potential that have smaller transverse scales show larger differences along the geomagnetic line between z = 0 and some other point on the line. Thus wide structures are propagated much more invariantly along B : structures larger than some kilometers arrive almost unaffected at so low latitudes that they almost reach the plasmasphere (as h > h , so “the same structure” is “lifted up”), dominating its dynamics. In Figure 8 a pictorial comparison between huge and small electrostatic structure propagation is showed.
The behaviour predicted by (246) is even more important on the transverse gradients of ϕ , i.e. on the fields E : from (242) we can write
E x , y , z = e 1 2 c 0 σ P z σ 0 z z k k f 0 k e i k x x + i k y y A + k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z + + A k 2 e 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z ,
while from (245) one may write
δ E k x , y , z δ E k x , y , 0 = = k F 0 k e i k x x + i k y y 1 2 c 0 2 + 4 k 2 B + k 2 B k 2 σ P z σ 0 z z + 1 2 c 0 B + k 2 + B k 2 σ P z σ 0 z z + O σ P z σ 0 z z 2 .
In the upper ionosphere and prodotonosphere the limit σ 0 σ H , P does hold almost perfectly, then the geomagnetic lines are perfect conductors, and the electric field is rigorously perpendicular to them
E / / = 0 ,
for any equilibrium current regime. The condition (249) renders the mathematics of the problem much simpler, but the most interesting phenomena from a ionospheric point of view take place where it does not hold at all. One of such cases is that discussed in [30], in which an abrupt time-change in σ induced by an exceptional photochemical event (a Gamma Ray Burst from the deep space) causes a non-trivial time evolution of E / / .

4.5.2. Current Injection and Perturbation Fields

In the same spirit of reasoning about equilibrium states, one can face the following problem: what happens injecting a current J q at z = 0 , as in Figure 9.
Some current J q is injected at some reference position z = 0 along a geomagnetic line, e.g. by a magnetic storm. This will yield a charge concentration perturbation and cause the presence of an electric field E all along the B -line with a delay of the order of lightspeed propagation time. The re-arrangement of the electric charges will lead very soon to a new electrostatic equilibrium, and we want to know the expression of E . . . , z once J q . . . , 0 is fixed. Note that the currents injected in magnetic storms always have a finite B -aligned component J z because the charged cosmic particles travel along the magnetic lines (at least their guiding centres [4]).
Again, the basic assumption is that the geomagnetic line is perfectly conducting, so that the approximation
E = E
holds. The useful relationships are then:
· J + z J z = 0 , ( ) E z = 0 .
From equation (247) one has, for σ P σ 0 0
lim σ P σ 0 0 E = k k f 0 k e i k x x + i k y y A + k 2 + A k 2 :
it is easily seen from here that in the limit σ P σ 0 0 the whole dependence on the geomagnetic coordinate disappears, so that the solution of (250) can be looked for with the condition:
z E = 0 .
The current orthogonal to the magnetic field is
J = σ · E ,
so that the continuity equation reads
· σ · E = z J z .
When this is integrated along the geomagnetic line from z = 0 up to some z 0 , one has:
0 z 0 · σ · E d z = J z 0 J z z 0 .
If the indices a and b are used to span the plane perpendicular to B we re-write it like this:
0 z 0 a σ a b E b d z = J z 0 J z z 0 .
Under the hypothesis (252) the field E b can be extracted from the integration, so that one has:
a E b 0 z 0 σ a b d z = J z 0 J z z 0 .
With the definition
Σ a b x , y , z 0 = def 0 z 0 σ a b x , y , z d z
and the assumption to arrest the integration (254) at some value z 0 where J z vanishes
J z z 0 = 0 ,
we can write:
a Σ a b . . . , z 0 E b = J z 0 .
Equation (257) is important because it puts in relation the electric field E taken in a point z 0 on a geomagnetic line with the value of the longitudinal current in z = 0 : due to the nature of the geomagnetic field, J z 0 is the field aligned current we can find in the polar oval, the flux of charged particles entering the ionosphere along the geomagnetic lines. It is possible to write:
x Σ P E x x Σ H E y + y Σ H E x + y Σ P E y = J z 0 .
If the ionosphere is assumed to be (locally) homogeneous, this becomes
Σ P . . . , z 0 x E x + y E y + Σ H . . . , z 0 y E x x E y = J z 0
and since E does not depend on z and an orthogonal decomposition of the divergence operator has been done, we have:
Σ P . . . , z 0 · E + Σ H . . . , z 0 y E x x E y = J z 0 .
The term Σ H y E x x E y is proportional to the component of × E along z, but since here we deal with a potential field, we have:
× E = 0 y E x x E y = 0 .
This turns (259) into:
· E = J z 0 Σ P .
From (252) we can say · E = · E so that:
· E = J z 0 Σ P .
Due to Gauss’ law of electrostatic
· E = δ ε 0
we can directly put the electric charge into (259), so we can relate the injected current J z with the electric charge density in the ionosphere
δ = ε 0 Σ P J z 0 .
The magnetic field B enters toward the Earth in the Northern region and comes out of it in the Southern one: than a positive J z in the North means a positive charge flux towards the ground; from equation (263) this means that δ has to be positive. In the Southern region some density δ > 0 in the ionosphere corresponds to currents flowing upward.

5. Conclusions

The aim of this review paper is essentially to write, in a clear way, all the equations concerning the evolution of the ionosphere, with a great care to illustrate their construction, focussing on each mathematical step, and highlighting each physical assumption.
The ionosphere is described as a system of mutually interacting fluids that exchange mechanical energy, momentum and heat; they are also subject to chemical reactions that alter their relative abundance. The possibility of turning one of these fluids into another via reactions is a very important characteristic of the ionospheric dynamics: from the definition of mass change rate G d C I given in (25), the whole sequence of the appearence of the production-loss terms Γ I descends, in the Navier-Stokes equations as well as in the internal energy balance.
In the program developed here it is important to highlight the depth of assuming that the LIM can be represented as a continuum, as well as to disclose its effects as far as possible: this was done in Section 2, where the system is presented.
As the LIM is formed by the continua C I , its dynamics is constructed basing on classical mechanics, equilibrium thermodynamics and classical electrodynamics, arriving to rich but compact self-consistent systems of coupled PDEs. These are the system (107) for matter, plus the quasi-stationary Maxwell equations (120). Usually, some simplifications are done to treat the LIM: a unique neutral continuum is introduced, evolving according to the PDEs in (140); moreover, the charged part of the LIM is represented by a unique positive ion species X + , plus free electrons e : these continua are subject to the dynamics (153) and (153) respectively.
As in many problems it can be of great use to describe the electrology of the LIM as that of a gaseous continuous conductor, the force-free assumption (153) is made, and a proper conductance tensor study is done, in Section 4, becoming able to study the ionospheric electrostatic equilibrium and some implications of it.
The content of this Review is a pre-requisite for studying ionospheric modelization from first principles: armed with the full set of fluid dynamical and thermodynamical equations describing the evolution of the continua C I , one is in the position to use suitable specializations of these differential relations in order to predict different phenomena of the ionospheric system.
It is convenient to conclude this review addressing the reader to further studies in the still partly unexplored land of post-fluid, or beyond-continuum, descriptions. Once again, we want to stress how demanding it is to require that an extended system of particles may be treated as a continuum: beyond the easy “smooth” picture emerging from Section 2, the dirty plasma forming the LIM shpould be treated via a multi-chemical kinetic theory, of through higher order truncations of the BBGKY hierarchy. Another approach of allowing extended singular regions to exist in the C I system, either considering the local proxies as fractal measures [31], or closing the dynamics with noise terms, mimicking local irregularities [32], instead of the smooth equilibrium hypotheses discussed in SubSection 3.5.
All these beyond-continuum pictures could probably enrich our description, and maybe comprehension, of “difficult” regimes, as ionospheric turbulence.

Appendix A Lagrangian time Derivative of the Internal Energy UdC I

Here we make the step-by-step calculation of the derivative of U d C I as indicated in (88). In order to do this, the heat capacity at constant volume of the parcel C V I T will be re-written as the product of the specific heat c V I of the I-th chemical component times the mass of the parcel m d C I
C V I T = c V I T m d C I .
The parcel internal energy then reads:
U d C I = m d C I T 0 T I c V I T d T + U O d C I .
With this expression we can write:
d U d C I = d t d d t m d C I T 0 T I c V I T d T + U O d C I
and consider carefully which quantities in U d C I do really depend on time. The addendum U O d C I surely does not, because the point O is some reference thermodynamical state chosen once and forever on the temperature-volume plane, as a reference. The mass m d C I does depend on time, due to the existence of production-loss processes, while in the integral
T 0 T I c V I T d T
the only quantity depending on t is the temperature T I placed at the end of the iso-volumic integration. Then one has:
d U d C I = d t d d t m d C I T 0 T I c V I T d T = d t d m d C I d t T 0 T I c V I T d T =
+ d t m d C I d d t T 0 T I c V I T d T =
= d t d m d C I d t T 0 T I c V I T d T + d t m d C I d d t T 0 T I c V I T d T = . . .
From (25) one may write
d U d C I = d t G I d C I T 0 T I c V I T d T + d t m d C I d T I d t T I T 0 T I c V I T d T =
= d t G I d C I T 0 T I c V I T d T + d t m d C I d T I d t c V I T I =
= m I Γ I T 0 T I c V I T d T + ρ I c V I T I d T I d t d t d V .
One finally writes down:
d U d C I = m I Γ I T 0 T I c V I T d T + ρ I c V I T I d T I d t d t d V .
Passing from the Lagrangian to the local quantities, this turns into (89).

Appendix B Navier-Stokes Equation for the Neutral Wind

Going from the second PDE in the system (131) to the equation of motion of the momentum ρ U of the neutral matter, is done considering:
h ρ h t V h + h ρ h V h · V h + h m h Γ h V h = h ρ h g e f f 2 h ρ h Ω × V h + h p h · h τ h V h , . . . h ρ h J h ν h J V h V J h · τ h w
and recognizing
h ρ h t V h + h ρ h V h · V h + h m h Γ h V h = = h ρ h t V h + V h · V h + h d ρ h d t V h = h ρ h d V h d t + + h d ρ h d t V h = h d d t ρ h V h = d d t h ρ h V h = d d t ρ U .
Then
d d t ρ U = h ρ h g e f f 2 h ρ h Ω × V h + h p h · h τ h V h , . . . h ρ h J h ν h J V h V J h · τ h w
On the right hand side of this equation some more algebra can be performed:
d d t ρ U = ρ g e f f 2 Ω × ρ U p + · h τ h V h , . . . h J h ρ h ν h J V h V J h · τ h w = . . .
The terms depending on collisions and friction with matter are
· h τ h V h , . . . h J h ρ h ν h J V h V J ,
and they should be re-written by distinguishing the neutrals from charged species in the second addendum
h J h ρ h ν h J V h V J = h k h ρ h ν h k V h V k + h Y ρ h ν h Y V h V Y :
here h and k are indices for neutrals, while I, J are general indices and Y is the general index for charged species. When all the neutral components are represented as a whole, the self-interaction of this whole producing a non vanishing part of its stress tensor at rest will be due partly to interactions within the single species, given by
· h τ h V h , . . . ,
and partly to interactions among different neutral species, given by
h k h ρ h ν h k V h V k .
Then we will define some neutral tensor τ 0 so that
· τ 0 = · h τ h V h , . . . + h k h ρ h ν h k V h V k :
this does represent the total neutral viscosity, allowing for the position:
d d t ρ U = ρ g e f f 2 Ω × ρ U p + · τ 0 h , Y ρ h ν h Y V h V Y h · τ h w .
The collective neutral stress tensor τ 0 can be represented in a form which is analogous to that of those τ h of the various neutral components [18], thus writing
· τ 0 = η 0 2 U + η 0 · U .
A collective stress tensor due to mechanical waves is easily defined as
τ 0 w = h τ h w ,
so that
d d t ρ U = ρ g e f f 2 Ω × ρ U p · τ 0 · τ 0 w h , Y ρ h ν h Y V h V Y .
The drag due to neutral-charged collisions
h , Y ρ h ν h Y V h V Y
should be re-written involving ρ and U : the following way is tried. We give a suitable definition of charged-neutral collision frequency ν 0 Y
h , Y ρ h ν h Y V h V Y Y ν 0 Y ρ U V Y
so that the charged-drag term reads:
h , Y ρ h ν h Y V h V Y = Y ν 0 Y ρ U V Y .
When this form is put into (A9) one finally has the complete Navier-Stokes equation for the neutral atmosphere, where the dynamics of the neutral species C h is represented as that of the fluid with local variables ρ and U , namely equation (132).

Appendix C Irrelevance of Γ e N e V → e in an Experimental Example

An example of how the approximation (147) is reliable, is the mid-latitude (rather quiet) ionosphere. Let us consider some free electron density evaluation obtained via the MIDAS algorithm [33,34], in which, at the point h , λ , φ = 150 km , 13 E , 40 N , on January 10, 2001, the following values were found:
N e 01 : 12 = 0.5 × 10 11 el / m 3 , N e 08 : 40 = 4.5 × 10 11 el / m 3
(the time is given in universal time). At that height and mid latitude the transport terms can be neglected, and a rough estimate of Γ e can be calculated as:
Γ e = Δ N e Δ t = N e 08 : 40 N e 01 : 12 Δ t .
Since Δ t = 26892 s , one has
Γ e = 4.5 0.5 26892 s × 10 11 el / m 3 1.49 × 10 7 el · Hz / m 3 .
The ratio Γ e N e can be calculated by assuming
N e = N e 08 : 40 + N e 01 : 12 2 = 2.5 × 10 11 el / m 3 ,
so that:
Γ e N e = 1.49 × 10 7 el · Hz / m 3 2.5 × 10 11 el / m 3 6 × 10 5 Hz .
The velocity V e gives another negligible contribution, when considered at the mid latitude of quiet ionosphere as indicated, so that the approximation (147) holds.

Appendix D The B →-parallel vs B →-perpendicular decomposition of velocities

Let us perform step-by-step all the calculations invoked to obtain the equation (200).
The parallel projector applied to W gives:
P i j / / B v j K ε j h k v h B k = P i j / / B v j K P i j / / B ε j h k v h B k = = B i B j B 2 v j K B i B j B 2 ε j h k v h B k = B i B j B 2 v j K B i B 2 v h ε j h k B j B k = = B i B j B 2 v j v i / / ,
so that one can write
V i / / = W j F B j B 2 B i
for every value of the parameters. The orthogonal projection gives:
P i j B v j K ε j h k v h B k = δ i j B i B j B 2 v j K ε j h k v h B k = = δ i j v j δ i j K ε j h k v h B k B i B j B 2 v j + B i B j B 2 K ε j h k v h B k = = v i B i B j B 2 v j K ε i h k v h B k + K B i B 2 v h ε j h k B j B k = = v i B i B j B 2 v j K ε i h k v h B k = δ i j B i B j B 2 K B k ε k i j v j
so that we have:
P i j B K B k ε k i j v j = P i j B W j F .
This is put in vector form:
v K v × B = W F .
Due to the very definition of cross product, only the part orthogonal to B will give a contribution to v × B
v × B = v × B ,
so that (A12) can be re-written as
P i j B v j K B k ε k i j P j h B v h = P i j B W j F ,
and with some manipulation of the indices one has
P i h B v h K B k ε k i j P j h B v h = P i m B W m F .
If one inserts some Kronecker delta
δ i j P j h B v h K B k ε k i j P j h B v h = P i m B W m F
we can write:
δ i j K B k ε k i j P j h B v h = P i m B W m F .
In vector form:
A B , M , ν · v = W F / A i j B , M , ν = δ i j K B k ε k i j .
The inverse of the matrix A B , M , ν must be found in order to obtain the velocity v : when some A 1 B , M , ν such that
A 1 B , M , ν A B , M , ν = 1
is defined, we can write
v = A 1 B , M , ν · W F .
Putting the (local) Cartesian coordinate system with the z axis along the magnetic field
B = B k ^ = B k ^
the tensor A B , M , ν can be written as the following matrix
A B , M , ν = 1 0 0 0 1 0 0 0 1 0 0 0 0 0 e B M ν 0 e B M ν 0 = 1 0 0 0 1 e B M ν 0 e B M ν 1 ,
whose inverse reads:
A 1 B , M , ν = 1 0 0 0 M 2 ν 2 M 2 ν 2 + e 2 B 2 ν e B M M 2 ν 2 + e 2 B 2 0 ν e B M M 2 ν 2 + e 2 B 2 M 2 ν 2 M 2 ν 2 + e 2 B 2 .
From definition (193) one has
A 1 K = 1 0 0 0 1 1 + K 2 K 1 + K 2 0 K 1 + K 2 1 1 + K 2
In order to put everything in a column-like formalism
v = v x v y v z
we will rewrite the tensors A and A 1 as:
A K = 1 K 0 K 1 0 0 0 1 , A 1 K = 1 1 + K 2 K 1 + K 2 0 K 1 + K 2 1 1 + K 2 0 0 0 1 .
For the component of v orthogonal to the magnetic field the solution (A16) has been obtained as in (200).

Appendix E The Electric Current Calculations

In this Appendix the expressions of the electric currents are obtained in the local NeWReF. All the job done in SubSection 4.1 is exploited here, considering the expressions of velocities in terms of the W I vectors, and the expressions of the latter ones in the NeWReF coordinates:
J q = N e V i V e =
= N e W i / / + 1 1 + K i 2 W i + K i 1 + K i 2 W i × B ^ +
W e / / 1 1 + K e 2 W e K e 1 + K e 2 W e × B ^ =
= N e W i / / W e / / + N e 1 1 + K i 2 W i 1 1 + K e 2 W e +
+ N e K i 1 + K i 2 W i K e 1 + K e 2 W e × B ^ = . . .
The two vectors W i and W e are given as
W i = D i H i g ^ e f f + b i E D i N N , W e = D e H e g ^ e f f + b e E D e N N .
This allows for the calculation of the various parts of J q , starting with the component along B :
N e W i / / W e / / =
= N e D i H i g ^ e f f + b i E D i N N / / D e H e g ^ e f f + b e E D e N N / / =
= N e 1 ν i 1 ν e g e f f / / + b i b e E / / D i N D e N N / / =
= N e 1 ν i 1 ν e g z + b i b e E z D i N D e N z N k ^ =
= N e 1 ν i 1 ν e g z D i N D e N z N k ^ + N e b i b e E z k ^ .
The function J z reads:
J z = N e 1 ν i 1 ν e g z e D i D e z N + N e b i b e E z .
Let us evaluate the first part orthogonal to the magnetic field:
N e 1 1 + K i 2 W i 1 1 + K e 2 W e =
= N e 1 + K i 2 D i H i g ^ e f f + b i E D i N N N e 1 + K e 2 D e H e g ^ e f f + b e E D e N N =
= N e 1 + K i 2 1 ν i g e f f + b i E D i N N N e 1 + K e 2 1 ν e g e f f + b e E D e N N =
= N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e g x i ^ + g y j ^ + N e b i 1 + K i 2 N e b e 1 + K e 2 E x i ^ + E y j ^ +
e D i 1 + K i 2 e D e 1 + K e 2 x N i ^ + y N j ^ ,
and archive it:
N e 1 1 + K i 2 W i 1 1 + K e 2 W e = = N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e g x i ^ + g y j ^ + + N e b i 1 + K i 2 N e b e 1 + K e 2 E x i ^ + E y j ^ + e D i 1 + K i 2 e D e 1 + K e 2 x N i ^ + y N j ^ .
The last addendum is evaluated as follows
N e K i 1 + K i 2 W i K e 1 + K e 2 W e × B ^ =
= N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g x i ^ + g y j ^ × k ^ +
+ N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E x i ^ + E y j ^ × k ^ +
e K i D i 1 + K i 2 e K e D e 1 + K e 2 x N i ^ + y N j ^ × k ^ =
= N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g x j ^ + g y i ^ +
+ N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E x j ^ + E y i ^ +
e K i D i 1 + K i 2 e K e D e 1 + K e 2 x N j ^ + y N i ^ = . . .
so that the following holds:
N e K i 1 + K i 2 W i K e 1 + K e 2 W e × B ^ = N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g y + + N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E y e K i D i 1 + K i 2 e K e D e 1 + K e 2 y N i ^ + N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g x + + N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E x e K i D i 1 + K i 2 e K e D e 1 + K e 2 x N j ^ .
The full expression for J x is:
J x = N e b i 1 + K i 2 N e b e 1 + K e 2 E x + N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E y + + N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e g x + N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g y + e D i 1 + K i 2 e D e 1 + K e 2 x N e K i D i 1 + K i 2 e K e D e 1 + K e 2 y N
while J y reads:
J y = N e b i 1 + K i 2 N e b e 1 + K e 2 E y N e K i b i 1 + K i 2 N e K e b e 1 + K e 2 E x + + N e 1 + K i 2 1 ν i N e 1 + K e 2 1 ν e g y N e K i 1 + K i 2 1 ν i N e K e 1 + K e 2 1 ν e g x + e D i 1 + K i 2 e D e 1 + K e 2 y N + e K i D i 1 + K i 2 e K e D e 1 + K e 2 x N .

Appendix F Solving the PDE for ϕ in the Anisotropic Conductor

In this Appendix the anisotropic Poisson PDEs introduced in SubSection 4.5 are discussed and solved.
Since the hessian tensor i j is symmetric, only the symmetric part of the matrix σ i j is needed
σ i j S = 1 2 σ i j + σ j i
then equation (241) reads:
σ i j S i j ϕ = 0 .
The symmetric matrix σ i j S can be diagonalized in suitable coordinates y i , and allows for real eigenvalues s i , so that
i s i y i y i ϕ = 0 .
When σ i j S is written as a diagonal tensor, another suitable transformation can be done
d x i = 1 s i d y i i s i y i y i = 2 .
The scalar potential ϕ will satisfy the proper Poisson equation in these new coordinates
2 ϕ = 0 .
For the conducibility tensor (219) we have
σ = σ P σ H 0 σ H σ P 0 0 0 σ 0 , σ S = σ P 0 0 0 σ P 0 0 0 σ 0 d y 1 d y 2 d y 3 = d x d y d z ,
so that the original coordinates, with z along B , diagonalize σ S , and the final coordinates are
d x d y d z = d x σ P d y σ P d z σ 0 .
With this reasoning the coordinates mapping equation (241) into (A31) can be found, once the tensor (219) is given.
When the ratios σ H , P σ 0 do become small, equation (A31) turns into
x 2 + y 2 ϕ = 0 , z 2 ϕ = 0 .
One possible assumption is that the potential only depends on the coordinates orthogonal to the magnetic line, that is equipotential as repeatedly announced before.
When inhomogeneous regions have to be treated, then the space dependence in σ is to be included in the PDE for ϕ , as in (239), and one may write
i σ i j j ϕ + σ i j i j ϕ = 0 :
if the expression (219) for the conducibility tensor is considered, one has
i σ i j j ϕ = x σ P + y σ H x ϕ + y σ P x σ H y ϕ + z σ 0 z ϕ , σ i j i j ϕ = σ P x 2 ϕ + σ P y 2 ϕ + σ 0 z 2 ϕ ,
so that
σ P x 2 ϕ + σ P y 2 ϕ + σ 0 z 2 ϕ + x σ P + y σ H x ϕ + y σ P x σ H y ϕ + z σ 0 z ϕ = 0 .
As the magnetic line determines some local cylindrical symmetry, it is then sensible to assume that the conducibility only depends on the coordinate along the line. Then
σ = σ z , σ P x 2 ϕ + y 2 ϕ + σ 0 z 2 ϕ + z σ 0 z ϕ = 0 .
If σ P 0 we can write
x 2 ϕ + y 2 ϕ + σ 0 σ P z 2 ϕ + 1 σ P z σ 0 z ϕ = 0
and also
1 σ P σ 0 z 2 ϕ + z σ 0 z ϕ = 1 σ P z σ 0 z ϕ .
The following equation is obtained
x 2 ϕ + y 2 ϕ + 1 σ P z σ 0 z ϕ = 0 ,
and when the transformation
d x d y d z = d x d y σ 0 σ P d z x y z = x y σ 0 σ P z x y σ 0 σ P z = x y z ,
is used in it, we have:
x 2 ϕ + y 2 ϕ + 1 σ P σ P σ 0 z σ 0 σ P σ 0 z ϕ = x 2 ϕ + y 2 ϕ +
+ 1 σ P σ P σ 0 z σ 0 σ P σ 0 z ϕ + σ 0 z σ P σ 0 z ϕ + σ 0 σ P σ 0 z 2 ϕ =
= x 2 ϕ + y 2 ϕ + 1 σ P σ P σ 0 z σ 0 σ P σ 0 z ϕ + 1 σ P σ P σ 0 σ 0 z σ P σ 0 z ϕ +
+ 1 σ P σ P σ 0 σ 0 σ P σ 0 z 2 ϕ = x 2 ϕ + y 2 ϕ + z 2 ϕ + 1 σ P σ P σ 0 z σ 0 z ϕ
+ σ 0 σ P σ P σ 0 z σ P σ 0 z ϕ = x 2 ϕ + y 2 ϕ + z 2 ϕ + z σ 0 σ 0 z ϕ + σ 0 σ P z σ P σ 0 z ϕ =
= x 2 ϕ + y 2 ϕ + z 2 ϕ + z ϕ z ln σ 0 + z ϕ z ln σ P σ 0 = x 2 ϕ + y 2 ϕ + z 2 ϕ +
+ z ϕ z ln σ 0 + z ln σ P σ 0 = x 2 ϕ + y 2 ϕ + z 2 ϕ + z ϕ z ln σ 0 + ln σ P σ 0 =
= x 2 ϕ + y 2 ϕ + z 2 ϕ + z ϕ z ln σ 0 σ P .
Then:
x 2 ϕ + y 2 ϕ + z 2 ϕ + z ϕ z ln σ 0 σ P = 0 .
Equation (A38) may be solved by using some function
ϕ x , y , z = f x , y g z ,
so that:
g z x 2 + y 2 f x , y + f x , y z 2 g z + z g z z ln σ 0 σ P = 0 .
By dividing everything by ϕ = f g one obtains
1 f x , y x 2 + y 2 f x , y + 1 g z z 2 g z + z g z z ln σ 0 σ P = 0
and the usual variable separation technique can be employed here, reminescent of the decomposition R 3 R × R 2 due to the presence of a magnetic field. The two addenda
1 f x , y x 2 + y 2 f x , y
and
1 g z z 2 g z + z g z z ln σ 0 σ P
must give zero when they are summed, for every possible value of x , y and z , so that the following position must be done:
1 f x , y x 2 + y 2 f x , y = C , 1 g z z 2 g z + z g z z ln σ 0 σ P = C , x , y , z R .
The second equation in (A41) may be solved straightforwardly: it is a linear equation with non constant coefficients
g z + z ln σ 0 z σ P z g z + C g z = 0 .
The expression of the square root σ 0 σ P is well fitted by the following mathematical function:
σ 0 z σ P z = c exp c 0 z .
One can finally write an expression deduced from the observation of what σ 0 z and σ P z behave like: σ 0 is basically uniform, while σ P behaves like an exponential, so their geometric average σ 0 σ P does behave like an exponential in turn
z ln σ 0 z σ P z = z ln c + c 0 z = c 0 .
Equation (A42) becomes a constant coefficient linear equation
g z + c 0 g z + C g z = 0 .
The homomorphic polynomial is
α 2 + c 0 α + C = 0 ,
with the two roots
α 1 = 1 2 c 0 + 1 2 c 0 2 4 C , α 2 = 1 2 c 0 1 2 c 0 2 4 C .
Equation (A44) is then solved by:
g z = A 1 exp 1 2 c 0 + 1 2 c 0 2 4 C z + A 2 exp 1 2 c 0 1 2 c 0 2 4 C z .
The function g can be written as a function of the old physical coordinate z:
g z = A 1 exp 1 2 c 0 + 1 2 c 0 2 4 C σ P z σ 0 z z + + A 2 exp 1 2 c 0 1 2 c 0 2 4 C σ P z σ 0 z z .
Let us come to the first of equations (A41)
x 2 + y 2 f x , y = C f x , y .
The function f x , y is simply some eigenfunction of the R 2 Laplacian operator x 2 + y 2 , with eigenvalue C given in (A41): in the particular case, this eigenvalue will be determined by theboundary conditions. f can be surely given via Fourier analysis, since the Laplacian
Δ 2 = x 2 + y 2 = x 2 + y 2
has eigenfunctions as
f k x , y = f 0 exp i k x x + i k y y .
Equation (A47) does then read
Δ 2 f = C f k 2 = k x 2 + k y 2 , k 2 f = C f ,
so that
C = k 2 ,
and
g k z = A 1 exp 1 2 c 0 + 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z + + A 2 exp 1 2 c 0 1 2 c 0 2 + 4 k 2 σ P z σ 0 z z .
The electrostatic potential that solves equation (A36) in the hypothesis (A43) can be written as the Fourier series, in equation (242) at page 52.

References

  1. Materassi, M. , Forte, B., Coster, A., Skone, S. (editors), The Dynamical Ionosphere - A Systems Approach to Ionospheric Irregularity, 2020 Elsevier Inc., ISBN: 978-0-12-814782-5.
  2. M. Materassi, G. Consolini, E. Tassi, Sub-Fluid Models in Dissipative Magneto-Hydrodynamics, in “Topics in Magnetohydrodynamics”, Editor L. Zheng, March 2012. [CrossRef]
  3. K. Davies, Ionospheric radio, Peter Peregrinus Ltd. 1990.
  4. M-B. Kallenrode, Space physics, Springer, 2000.
  5. Vincenzo Carbone, Mirko Piersanti, Massimo Materassi, Roberto Battiston, Fabio Lepreti, and Pietro Ubertini, A mathematical model of Lithosphere-Atmosphere coupling for seismic events, Scientific Reports, (2021) 11:8682. [CrossRef]
  6. L.D. Landau, E.M. Lifshiz, Mcanique des fluides, ditions Mir, Moscou 1971.
  7. A.R. Choudhuri, The Physics of Fluids and Plasmas, Cambridge University Press, 1998.
  8. van der Vaart, A.W. , Asymptotic Statistics, 9780521496032, Cambridge series on statistical and probabilistic mathematics, https://books.google.it/books?id=fiX9ngEACAAJ, 1998, Cambridge University Press.
  9. T. L., Hill, Thermodynamics of Small Systems, The Journal of Chemical Physics, Vol. 36, No. 12, June 1962.
  10. Gregory L. Eyink, Katepalli R. Sreenivasan, Onsager and the theory of hydrodynamic turbulence, Rev. Mod. Phys., Vol. 78, No. 1, January 2006.
  11. N. Taccetti, Course of Electromagnetism and Optics, academic year 1991/1992, University of Florence.
  12. M. Materassi, Variabili canoniche collettive e relative per il campo di Klein-Gordon classico, master degree thesis, University of Firenze (Italy), 1996.
  13. R, Bǎlescu, Statistical Dynamics - Matter Out of Equilibrium, April 1997, doi.org/10.1142/p036, World Scientific Publishing.
  14. R. Tozzi, A. De Sanctis, Quaderni di Geofisica n.8, 2000, Evoluzione di sistemi dinamici non lineari, il caso del campo geomagnetico and references therein.
  15. K. Schlegel, J.P. St.-. K. Schlegel, J.P. St.-Maurice, Anomalous heating of the polar E region by unstable plasma waves. 1. Observations, Journal of Geophysics Research, 1981, 86, 1447.
  16. A. Fasano, S. Marmi, Meccanica analitica, Bollati Boringhieri, 1992.
  17. R. Leitinger et al., The upper atmosphere, Springer 1996.
  18. M.C. Kelley, The Earth’s Ionosphere, Academic Press Inc., 1989.
  19. A.S. Jursa, Handbook of geophysics and the space environment, Air Force Geophysics Laboratory - Air Force Systems Command - USA Air Force, 1985.
  20. E.I. Yakubovich, D.A. Zenkovich, Matrix fluid dynamics, arXiv: physics/0110004 v1 and references therein.
  21. I. Antoniou, G.P. Pronko, On the hamiltonian description of fluid mechanics, arXiv: hep-th/0106119 v2 and references therein.
  22. P. Lynch, Hamiltonian methods for geophysical fluid dynamics: an introduction, IMA, University of Minnesota, Preprint 1838 (February 25, 2002). Webpage: https://maths.ucd.ie/~plynch/Publications/HMGFD.pdf.
  23. M. Materassi, Metriplectic Algebra for Dissipative Fluids in Lagrangian Formulation, Entropy 2015, 17(3), 1329-1346. [CrossRef]
  24. J.A. Dutton, Dynamics of the atmospheric motion, Dover Publications Inc., 1995.
  25. M.W. Zemansky, Heat and thermodynamics, McGraw-Hill Book Company, Inc., 1957.
  26. S. Weinberg, The quantum theory of fields, vol.1, Cambridge University Press 1995.
  27. M. Materassi, C.N. Mitchell, J.P.S. Spencer, Ionospheric imaging of the northern crest of the Equatorial Anomaly, Journal of Atmospheric and Solar-Terrestrial Physics, Volume 65, Issues 16–18, November–December 2003, Pages 1393-1400. [CrossRef]
  28. H. Rishbeth, O.K. Garriott, Introduction to Ionospheric Physics, Academic Press, New York and London, 1969.
  29. C.A. Gonzales, M.C. Kelly, D.L. Carpenter, T.R. Miller, R.H. Wand, Simultaneous measurements of ionospheric and magnetospheric electric fields in the outer plasmasphere, Geophys. Res. Lett. 7, 517 (1980).
  30. M. Materassi, M. Piersanti, G. D’Angelo, E. Papini, G. Consolini, A. Bazzano, J. G. Rodi, P. Ubertini, Gamma-ray burst induced fast dynamics of the electric field in the top-side ionosphere, in preparation.
  31. Materassi, M. and Consolini, G.: Turning the resistive MHD into a stochastic field theory, Nonlin. Processes Geophys., 15, 701–709. [CrossRef]
  32. M. Materassi, G. Consolini, Magnetic Reconnection Rate in Space Plasmas: A Fractal Approach, Phys. Rev. Lett. 99, 175002 – October, 2007. [CrossRef]
  33. M. Materassi, Study of the Northern Crest of the Equatorial Anomaly via a Multi Instrumental Data Analysis System, internal research report of I.R.O.E. ”Nello Carrara” CNR, RR/ATM/02.01.
  34. M. Materassi, L. Ciraolo, P. Spalla, C.N. Mitchell, P.S.J. Spencer, The effect of the Northern Crest of the Equatorial Anomaly on the propagation delay at GPS frequencies, in the proceedings of the GNSS ’2001 Congress, Sevilla 10 May 2001.
Figure 1. A cartoon of the smooth evolution of the classical continuum C from the initial configuration C t 0 to a later one C t . If ψ K is any local proxy describing the system (e.g., its mass density), the continuous hypothesis is that the function ψ : D t C × R K is mathematically smooth, both with respect to the local coordinate r D t C and to time t R . In the picture, the border of D t 0 C and D t C are highlighted through a bold black curve, namely indicating D t 0 C and D t C .
Figure 1. A cartoon of the smooth evolution of the classical continuum C from the initial configuration C t 0 to a later one C t . If ψ K is any local proxy describing the system (e.g., its mass density), the continuous hypothesis is that the function ψ : D t C × R K is mathematically smooth, both with respect to the local coordinate r D t C and to time t R . In the picture, the border of D t 0 C and D t C are highlighted through a bold black curve, namely indicating D t 0 C and D t C .
Preprints 144566 g001
Figure 2. A sketch of the “macroscopic infinitesimal”, i.e. “parcel”, d C of the continuum C , containing a certain number of particles, indicated as red dots inside it, the centre-of-mass of which is the point P. The centre-of-mass velocity V of all the particles within this d C P is attributed to P as V P .
Figure 2. A sketch of the “macroscopic infinitesimal”, i.e. “parcel”, d C of the continuum C , containing a certain number of particles, indicated as red dots inside it, the centre-of-mass of which is the point P. The centre-of-mass velocity V of all the particles within this d C P is attributed to P as V P .
Preprints 144566 g002
Figure 3. The cartoon of how smoothness of macroscopic variables in FD should be intended: d C P and d C P are two “nearby” parcels, as the one portrayed in Figure 2, with centre-of-mass P and P respectively, separated by a “small” displacement δ r . d V and d V are respectively the “small” volumes of d C and d C .
Figure 3. The cartoon of how smoothness of macroscopic variables in FD should be intended: d C P and d C P are two “nearby” parcels, as the one portrayed in Figure 2, with centre-of-mass P and P respectively, separated by a “small” displacement δ r . d V and d V are respectively the “small” volumes of d C and d C .
Preprints 144566 g003
Figure 4. The O N A particles of the I-th species within the fluid parcel: the red dots represent the elementary components of the I-th species, while the red continuous line is the trajectory of d C I as a whole.
Figure 4. The O N A particles of the I-th species within the fluid parcel: the red dots represent the elementary components of the I-th species, while the red continuous line is the trajectory of d C I as a whole.
Preprints 144566 g004
Figure 5. The NeWReF versus the SERF: one defines the NeWReF through the transformations (180) and (183), then one puts the local Cartesian axes with z along B , while the so-oriented NeWReF translates rigidly with respect to the SERF (without rotating). By the way, consider that the axes used in the transformations (184) and (185) are not the ones used in SubSection 4.2, where z ^ is parallel to B and not to U .
Figure 5. The NeWReF versus the SERF: one defines the NeWReF through the transformations (180) and (183), then one puts the local Cartesian axes with z along B , while the so-oriented NeWReF translates rigidly with respect to the SERF (without rotating). By the way, consider that the axes used in the transformations (184) and (185) are not the ones used in SubSection 4.2, where z ^ is parallel to B and not to U .
Preprints 144566 g005
Figure 6. Cartoon of a geomagnetic line of the field B with the electrostatic ionospheric field E : as the conductivity along B is always much higher than that perpendicular to it in the ionosphere, at the electrostatic equilibrium vectors E will lie perpendicularly to the lines of B .
Figure 6. Cartoon of a geomagnetic line of the field B with the electrostatic ionospheric field E : as the conductivity along B is always much higher than that perpendicular to it in the ionosphere, at the electrostatic equilibrium vectors E will lie perpendicularly to the lines of B .
Preprints 144566 g006
Figure 7. Sketching the prevalent component of the geomagnetic field, which is dipolar, one sees that a point P of height h taken along the field line at high latitudes, dragged along the same field line down to low latitudes, will result in a corresponding position P of bigger height h > h .
Figure 7. Sketching the prevalent component of the geomagnetic field, which is dipolar, one sees that a point P of height h taken along the field line at high latitudes, dragged along the same field line down to low latitudes, will result in a corresponding position P of bigger height h > h .
Preprints 144566 g007
Figure 8. Comparison between the “conservation” of wide electrostatic structures transported along a geomagnetic line from the Polar to the Equatorial Region, and the “dispersion” of small electrostatic structures, see the text.
Figure 8. Comparison between the “conservation” of wide electrostatic structures transported along a geomagnetic line from the Polar to the Equatorial Region, and the “dispersion” of small electrostatic structures, see the text.
Preprints 144566 g008
Figure 9. The cartoon of the system described in SubSection 4.5.2, in which some electric current is injected along B at high latitudes.
Figure 9. The cartoon of the system described in SubSection 4.5.2, in which some electric current is injected along B at high latitudes.
Preprints 144566 g009
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated