Study of a Gas Disturbance Mode Content Based on the Measurement of Atmospheric Parameters at the Heights of the Mesosphere and Lower Thermosphere

: The main result of this work is the estimation of the entropy mode accompanying a wave disturbance, observed at the atmosphere heights range of 90–120 km. The study is the direct continuation and development of recent results on diagnosis of the acoustic wave with the separation on direction of propagation. The estimation of the entropy mode contribution relies upon the measurements of the three dynamic variables (the temperature, density, and vertical velocity perturbations) of the neutral atmosphere measured by the method of the resonant scattering of radio waves on the artiﬁcial periodic irregularities of the ionospheric plasma. The measurement of the atmosphere dynamic parameters was carried out on the SURA heating facility. The mathematical foundation of the mode separation algorithm is based on the dynamic projection operators technique. The operators are constructed via the eigenvectors of the coordinate evolution operator of the transformed system of balance equations of the hydro-thermodynamics.


Introduction
The idea of a fluid perturbation decomposition goes up to the celebrated paper [1] related to a turbulent flow and studies of the nonlinear viscous and heat conducting compressible gas were formulated in Reference [2], where the wave and non-wave components notions and a perturbation were formulated. Interaction of the components was also introduced.
A general description of a fluid perturbation in respect to equilibrium is naturally divided by projections to linear evolution operator subspaces to be identified as the perturbations modes. Its mathematical formulation is given in Reference [3] for the classic (1, 2, 3 kinds) boundary conditions; see the recent book by Reference [4], as well. Its spectral content is determined by a kind of Fourier transform [5]. There are wave nonzero frequency ω and non-wave (ω = 0) components [6]. The significance of the non-wave (entropy) component of a gas disturbance relates directly to the atmosphere warming phenomenon, for sudden stratospheric events [7] and for warming compendium look [8].
Appearance of the entropy mode is connected either with a heating phenomenon when energy transfer from waves of a large amplitude is implied [4] or when such mode is transported by wind or large scale (e.g., planetary) waves [9,10].
Analogous problems in plasma physics is more complicated due to the basic field components abundance with applications to the sun atmosphere physics [11,12]. A direct theory of the heating phenomenon at plasma with mechanical and thermal losses and electrical resistivity is presented in Reference [13]. The magnetoacoustic heating of plasma caused by periodic magneto-sound perturbations with discontinuities is studied in a quasi-isentropic magnetic gas [14].
In algorithmic approach to mode definition [15,16], the method is formulated via projecting operators, that are defined for each evolution operator subspace. The projecting operators and their energy weight calculation are realized by means of an auxiliary norm introduction. The whole algorithm stages of a disturbance mode content diagnosis are outlined below for the case of exponential atmosphere.
1. There are two possible classes of evolution operators that may be a base for the projecting procedure construction. One of them is evolution in time-say, t−evolution,-for its applications [17]. The second is the z-evolution operator, both acting over the vector state. In our 1D case, it is a three-component one ψ T = (U, p, ρ) with the gas velocity (U), pressure (p), and density (ρ). The z-evolution equation reads namely, such a form that is derived in Section 3 defines the evolution operator L. The notations of the section differ because the dimensionless coordinates and entropy perturbation (instead of temperature) are introduced for further convenience. 2. Let consider now the eigenvalue problem for the t-Fourier transformL of the evolution operator L having the matrix formLψ where ψ →ψ is the Fourier transform of the state vector. 3. For each eigenspaceψ i , i = 1, 2, 3, it is convenient to introduce the matrix projecting operators P I by means of the 3 × 3 eigen matrix Φ αβ = φ 1 α φ 2 β , α, β = 1, 2, 3 so that the matrix elements [4] P i αγ = Φ αβ Φ −1 βγ , are the functions of the frequency ω. 4. After the inverse Fourier transformation, one can apply the projecting technique to the observation data; otherwise, we transform the data by the discrete Fourier transformation; at the next step, we use the corresponding discretized projecting matrices P i γδ to extract the data for the three modes (directed waves and entropy mode). 5. The energies of the directed and the entropy mode are estimated via discretized form of the general expression for energy of a gas perturbation. Its expression follows basic equations of the conservation law.
A role of the zero frequency mode is specific because its propagation is determined by the simultaneous action of the nonlinearity and dissipation. Hence, the mode (entropy mode) diagnostics by local measurements is a challenge from both theoretical point of view (asymmetry of derivatives entrance in basic equations) and from observations, due to its slow evolution [4].
While the entropy mode account enters the main mechanism for some physical processes, such as non-linear heating, in other contexts, it can be safely dropped out. In addition, because of the difficulty of its measurement and estimation, the majority of models do just that. Generally, it is expected that it should be relatively insignificant in this case of the heights range under investigation. However, we still want to obtain the value of it and to check those expectations. The proposed method could also be applied to other problems where its contribution is more pronounced. Even in the situation we examine, the mode growth with altitude may, potentially, lead to the essential change of the content of the atmosphere perturbation.
To solve this problem, we rearrange equations of the basic system so as to be a z-evolution one (Section 1 and Reference [10]). This will allow us to find a projection operators for all three modes by the mentioned algorithm. In addition, we will derive formula for an energy of the mode in the ω-space. Such results open the way to the local space diagnostics of acoustics and entropy modes.
For a source of observation data, we take measurements of the velocity, temperature and the density of the neutral atmosphere in the height range 90-120 km by API (artificial periodic irregularities) technique. The technique is based on the resonant scattering of radio waves by induced periodic irregularities of the ionospheric plasma emerging in the field of a standing wave arising from the interference of the incident and reflected waves from the ionosphere [18,19]. The details of the physical base of the API is described in Section 2. The mathematical statement of the problem is given in the Section 3. Its reformulation as a z-evolution problem is derived in Section 4. An expression for energy density are given at Section 6. The errors of measurements is given at Section 7 and the error of the whole estimation algorithm is described at the Section 8.
The main idea of this work is to prove that the data on atmosphere parameters extracted from the SURA facility allow to estimate contributions of three modes (two waves and entropy). Other sources do not give such an opportunity.
Note that, in Reference [15], it was presented as an alternative approach based on discrete projection operators. In addition, results in this paper were not applied to experimental data. The problem which was described there had only two wave modes without taking the entropy mode into account. The approach used to calculate energy in Reference [15] was also more complex, and we were able to find the means to greatly simplify this task.
Reference [16] also describes a two mode model and does not contain ways to calculate an error estimation. This paper, while it is a logical continuation of the work done in Reference [16], describes a more precise three mode model, which allows for analyzing entropy mode, usually a very hard task to do with real experimental data. It also includes error estimations, both from experimental errors and ones from numerical method errors.
A more general basis of mathematics used in projection operators methods can be seen in book by Reference [17], as applied to several different fields of physics.

Experimental Background
The estimation of the entropy mode contribution relies upon the measurements of three dynamic variables (the temperature, density, and vertical velocity perturbations) of the neutral atmosphere obtained by the method of the resonant scattering of radio waves on the artificial periodic irregularities (APIs) of the ionospheric plasma [18][19][20]. The measurement of the atmosphere dynamic parameters has been carried out on the SURA heating facility (56.11 o N; 46.1 o E), which is situated near Nizhny Novgorod, Russian Federation. The facility creates artificial disturbances in the ionosphere influenced radio wave propagation. Powerful radio waves radiated into the ionosphere disturb the temperature and electron density. This causes many different "heating" effects in the ionospheric plasma, which are described in detail in reviews in Reference [21,22]. The creation of artificial periodic irregularities in the field of a powerful standing radio wave is among them. The API technique is described in detail in Reference [18]. The scattering of probe radio waves by these irregularities has resonant properties, that is, the signal received from the API has significant amplitude when the frequencies and polarizations of the powerful and probe radio waves are the same. After the end of the artificial impact on the ionosphere, APIs disappear during the relaxation time. A large level of scattered signal with a signal-to-noise ratio of the order of 10-100 allows one to determine atmospheric parameters with high accuracy. The altitude-time resolution is 1 km in altitude and 15 s in time and makes it possible to study short-term and long-term processes in the ionosphere and in neutral atmosphere.
The method for determining the temperature and density of the neutral atmosphere based on the analysis of the height dependence of the relaxation time of the API scattered signal is described in detail in Reference [18,19]. At altitudes of 90-120 km, the API relaxation occurs under the influence of the ambipolar diffusion [18,22]. The amplitude and phase of the scattered signal at the stage of relaxation of inhomogeneities are measured. The relaxation time is determined by the expression: where K B is the Boltzmann constant, K = 4π/λ is the wavenumber of the standing wave, λ = λ 0 /n is the wavelength in the propagation environment, n is the refractive index, D is the ambipolar diffusion coefficient, M i is the molecular mass of the ion, T e0 and T i0 are the unperturbed electron and ion temperatures, and ν im is the frequency of collisions between ions and neutral molecules. It is considered that, in the mid-latitude at mesosphere and the lower thermosphere heights, T e0 = T i0 = T, up to the height 120-130 km, where T is the temperature of the neutral component. The expression for τ * underlies the determination of many parameters of the lower ionosphere and the neutral atmosphere, including its temperature and density [18,19]. The atmospheric temperature and density are determined from the API relaxation time, which is inversely proportional to the ambipolar diffusion coefficient. Many results of studying the altitude-temporal variations of temperature and density in different geo-and helio-physical conditions are presented in Reference [23][24][25][26][27][28].
To learn more about the API technique and, in detail, the methodology for determining the atmospheric temperature and density, we recommend reading Reference [18,19]. The velocity of the vertical motion of the plasma is determined from the dependence of the phase on time. This velocity in the mesosphere and lower thermosphere is equal to the velocity of the neutral component, since the plasma at these altitudes is a passive admixture and moves with the neutrals. The measured phase φ of the scattered signal can be written as dφ/dt = 2πF d = 4πV/λ, where F d is the Doppler velocity, and V is the velocity of the vertical movement of the neutral medium. Then, the vertical velocity can be determined by the formula: where f is the frequency of the powerful and probing waves. Positive velocity values correspond to downward movement. To eliminate random errors associated with a change in the phase of the scattered signal due to the influence of natural noise on the scattered signal, the linear part of the approximation of the φ(t) dependence is used for the velocity calculation.

Mathematical Statement of the Problem Formulation
Following Reference [6], we start from the linearized conservation equations describing one-dimensional flow along the vertical axis z in terms of the deviations of pressure p , density ρ , and velocity U , from corresponding equilibrium stationary values p, ρ, U. In the exponentially stratified atmosphere, they take the form: where H is the height scale of stratification, p 0 and ρ 0 are the values of pressure and density at the zero z level, and g is gravity acceleration. This system can be written in the form: ∂p ∂t In the the one-dimensional exponentially stratified model presented in this paper, the external force is described as the constant gravity acceleration g, which is directed opposite to axis z, though it can be generalized to other mass forces, including non-inertial ones. The flow of an ideal gas is considered, whose internal energy , in terms of pressure and density, takes the form: where p = p + p , ρ = ρ + ρ , i.e., full value is the sum of unperturbed value plus perturbation, and γ = C p /C v denotes specific heat ratio. The relation between the equilibrium pressure and density follows from the zero order stationary equality: These linearized 1D equations for perturbations of pressure p , density ρ , and gas velocity U now can be simplified by introducing the new variable φ = p − γ p ρ ρ , and then we can switch to these variables, conventionally used for the exponentially stratified atmosphere model: where z is vertical coordinate, and H is the height of the stationary atmosphere. That transforms the system ((4)-(6)) to the form: where ρ 0 denotes the density at the model lower boundary z = 0. System ((10)- (12)) describes a three-mode problem with two wave modes and one stationary entropy mode. We can rewrite system ((10)-(12)) using the variables which are measured by API technique using relations: ρ 0 = p ρ in exponentially stratified atmosphere and Mendeleev-Clapeyron's law (V * is the gas volume) and, since p = R µ Tρ in Mendeleev-Clapeyron's law (µ, there is an average gas molar mass) for nonpertubed atmosphere and after dropping out non-linear terms, we get It is convenient to rewrite the system ((10)- (12)) in terms of dimensionless functions and variables. To do this, we shall use the uniform atmosphere height H and the velocity of sound c = γgH as dimension parameters, which gives time scale H/c = H γg so that the new dimensionless variables are z = Hξ, t = τ · H/c = τ H γg . Functions are redefined as U = cV = V γgH, P = p 0p , since the Φ, too, has the pressure as dimension (because φ = p − γ p ρ ρ and Φ = φ · e z/2H ), and Φ = p 0φ . After isolation of terms with ∂ ∂τ and taking into account that p 0 = gHρ 0 , we will get: ∂p ∂τ

Atmosphere Gas Disturbance Initiation by a Boundary Regime
Let us transform our task into the boundary mode z-propagation form. It is easy to isolate spatial differentials forp and V: and but the initial equation for ∂φ ∂τ , (17) does not contain ∂ ∂ξφ . We can obtain it from the third equation in the form: ∂φ ∂ξ The integral should not be a problem since it has a rather simple form after Fourier transformation. Finally, the system (18)- (19) and (21) takes the form where L is the z-propagation matrix operator, that gains the form where I τ κ = κdτ. To obtain the mode projection operators, we need to transform (22) into Fourier space. At the domain, the propagation matrix takes the form:

Projection Operators
Propagation matrixL (24) has the eigenvalues 0, 1 which are consistent with expected two directed waves, with up-and downward directions, and one stationary entropy mode. We can calculate projection operators using the standard method P k ij = E ik E −1 kj , where E is a matrix constructed from eigenvectors. They take the form: where K = ), X ± = ±2γω 2 + k ± 1. Variables K, Z ± , and X ± were introduced only to shorten notation of projection operators.

The Mode Energy Density at Omega-Domain
Following Reference [6], we get the formula for the energy density: in the three-mode case. This formula is written using variables from (z, t)-space, so we need to rewrite it, using Fourier-images of these variables from (ω, τ)-space. Firstly, we must transform it into dimensionless form. To do it, we will use the same formulas that we used for main system: z = Hξ, t = τ · H/c = τ H γg ,P = p 0p = ρ 0 gHp, U = cV = V γgH and Φ = p 0φ = ρ 0 gHφ. We get the expression: A transition to the ω-space variables by the conventional Fourier transformations: and We plug it in the expression for energy density (30) and get: Hence, the energy ξ−profile is obtained by the integration over dimensionless time: It depends on ξ instead of τ since we are solving boundary mode ξ-propagation problem. However, by definition, hence,

Neutral Atmosphere Parameters Obtained Using API Technique
As in Reference [16], we used the results of determining atmospheric temperature and density on 26 September 2017 from 12:00 to 16:20 local time (Moscow time) at the height at 100 km, adding to them the results of simultaneous measurements of the vertical velocity. Experimental data take the form of measurements of temperature, density, and vertical gas velocity. Our projection operators are formulated in the terms of V,p, andφ. Vertical gas velocity can be used after calculating difference from average, transforming into dimensionless variables using formula U = V γgH. Then, the discrete Fourier transformation can be applied.
For the wave entropy variableφ, we haveφ =p − γ ρ 0 ρ, where ρ is perturbations of density, and it can be calculated as deviation from average, too.
So, the only question is how to go over from temperature top, which is the perturbation of pressure. To do that, we calculate pressure in each point using Mendeleev-Clapeyron law: and then we calculate the difference with the averaged value over all points to get perturbations values. Experimental data are presented in Table 1. Now, it is possible to apply out projection operators (26,27,28) and then to calculate, separately, energies for each mode, using formula (37), which need to be rewritten in the discrete form by replacement of integrals with sums, i.e., using midpoint rule. Since points are spaced, and we are interested in relative values only, simple sums can be used. It gives us relative values. For entropy mode, we obtain and, for the wave modes, we have where E = E 1 + E 2 + E 3 is the total energy of all modes.
In contrast to Reference [16], we used the data of determination of three atmospheric quantities; we use a set that allows us to find the contribution of the entropy mode. An important fact is that it turned out to be small. In addition, we have demonstrated the procedure for measuring this contribution. Recall that the origin of the entropy mode is heating by waves, which is proportional to the dissipation and the square of the amplitude of this mode. At these altitudes, both are small.

The Vertical Velocity Measurement Error
Our analysis has been carried out with a certain degree of accuracy. The error in the E estimates is the sum of the errors in measuring the parameters of the neutral component (temperature, density, and vertical velocity) by the API technique and the error in the algorithm for calculating the entropy and wave modes.
As can be seen from the formula for the vertical velocity, the absolute error ∆V of a single measurement of the speed is determined by the error of the time derivative of the signal phase. If we find ∆φ/∆t from two dimensions, then the variance of this derivative is expressed by the formula: where the time measurement error is neglected, σ 2 φ 1 and σ 2 φ 2 are the variances of the first and second measurements, and t 1 is the measurement time.
To reduce the error, it is necessary to increase the measurement time interval t 1 , but this is prevented by the exponential decrease in the signal amplitude during the API relaxation period, which leads to a decrease in the signal-to-noise ratio (A/A noise ) and an increase in the variance σ 2 φ 2 . According to Reference [29], for A/A noise > 3, the value of σ 2 φ is σ 2 This value reaches a minimum at t 1 = 0.86τ * and gives the value σ 2 ∆φ/∆t = 3σ φ 1 . A more accurate estimate of the derivative can be obtained by taking into account all measured values of the phase φ. It was shown in Reference [20] that the optimal procedure for finding ∆φ/∆t is a linear approximation of φ(t) by the least squares method with the weight function exp(−t 1 /τ * ), while at t 1 τ * . Ultimately, the error in determining the vertical speed is expressed by the formula: Assuming the relaxation time of the scattered signal equal to τ * = 1 s the time between two measurements ∆t = 20 ms and A/A noise = 20, which corresponds to the API experiments, we find the value ∆V = 0.08 = 0.08 m/s. Thus, the random error of a single measurement does not exceed 0.1 m/s. There is a small systematic error due to the heating of the electron gas by a powerful radio wave. This issue is considered in detail in Reference [18]. It was shown there that, for typical experimental conditions at an altitude of 100 km, the systematic error in measuring the velocity does not exceed ∆V = 0.05 m/s for the extraordinary component of a powerful wave. The total relative velocity error does not exceed 2%.

Errors in Determining Neutral Temperature and Density
In Reference [18,19], a method for determining the temperature and density of the neutral component based on measurements of the height-time dependencies of the amplitude of the signal scattered by the API technique is considered in detail. Let us briefly explain how the errors in determining the atmospheric parameters were calculated. Based on the expression of the barometric dependence of pressure on altitude in an isothermal ρ(H) = ρ 0 exp [−(h − h 0 )/H], H = K B T/mg, where H is the height of a homogeneous atmosphere, K B is the Boltzmann constant, and measuring the relaxation time τ * at two close (to satisfy the isothermal condition in height) heights h 1 and h 2 , taking into account the dependence of the wave vector k on the height h through the refractive index n of a powerful wave in plasma k = k 0 n(h), we obtain for H the expression: The values of the relaxation time τ * are determined by the decay of the amplitude of the scattered signal by a factor of e with a step in height of 1 km.
The neutral temperature T is determined from the relation T = mgh/K B , respectively, the absolute error in determining the temperature is ∆T = mg∆h/K B , and the relative error is δT = δH. To find the value of H, a linear approximation of the function ln (n 2 τ * ) = b 0 h + b 1 over several values of the relaxation time is used. Coefficients b 0 and b 1 are determined by the least squares method by applying the linear regression method to the dependence ln (n 2 τ * ) In this case, the standard (root-mean-square) approximation error is calculated in the usual way [30]. For most measurements by the API method, the relative error δT = δH does not exceed 5%. The error in determining the density ρ is of the same order of magnitude and does not exceed 10% [19,20]. As the long-term measurements of atmospheric parameters by the API method show, only in some cases these errors can reach 20%.

Algorithm Error
The algorithm described in this paper inevitably also introduces an additional error. It is impossible to definitely restore function by the limited set of points without making some propositions about function, which is not always possible to do. We will try to evaluate which error is introduced by this process of restoration. This error appears even for pure waves.
To estimate it, we used a method based on Runge's approach. We recalculated all our energies using only a subset of our data points. The resulting relative error was the estimated as: There, E N denotes the value of energy calculated using N data points. E N/2 , using only odd numbered points. Errors do not noticeably change if we use complimentary subset of points (i.e., even numbered points), which confirm correctness of such a method of error estimation.
Resulting error was around 49% for entropy mode and around 0.9% for wave modes. Therefore, and, for wave modes We have used the results of the energy density estimation at the heights within the nearest vicinity of the experimental conditions (about z = 100 km). The estimation is performed in the frequency domain, which we consider as the upper boundary of the error because of the minimum contribution of white noise while the Fourier transformation is performed.
The relatively large value for error of the entropy mode estimation is the direct result of smallness of its value. In any case, it shows that it is a lot smaller than both wave modes and can be non-issue for a wide range of problems.

Conclusions
The main purpose of this work is the estimation of the entropy mode accompanying a wave disturbance, observed at the atmosphere heights range of 90-120 km. The study is the direct continuation and development of recent results on diagnosis of the acoustic wave with the separation on direction of propagation. The estimation of the entropy mode contribution relies upon the measurements of the three dynamic variables (the temperature, density, and vertical velocity perturbations) of the neutral atmosphere measured by the method of the resonant scattering of radio waves on the artificial periodic irregularities of the ionospheric plasma. The measurement of the atmosphere dynamic parameters was carried out on the SURA heating facility. The mathematical foundation of the mode separation algorithm was based on the dynamic projection operators technique. The operators were constructed via the eigenvectors of the evolution operator of the transformed system of balance equations of the hydro-thermodynamics.
The main result is the algorithmic extraction of the modes of an atmosphere gas disturbance, specifically: the entropy mode from three-component data of atmospheric parameters. Having in mind the vertical movements of the atmosphere gas within a rather small height interval, we applied the 1D exponential model of the atmosphere that results in a relatively simple algorithm of the mode contribution extraction. The entropy mode energy estimation, as calculated in Section 7, gave small values compared to ones of wave modes. We consider it an important result, expecting, however, the mode contribution growth, that may be an important reason for the background temperature and density slow dynamics, which, as mentioned, should be taken into account in thermosphere heating or cooling during strong atmospheric storms. A reason for such growth may be a presence of wave perturbation of high enough amplitude. It is important for atmosphere perturbations modeling, as for prognosis.