Considering a one-dimensional wave at time moment t , the particle number density at x in the medium is given by n x , t . Under the conservation of particle number, there is d n d t = 0 , (2) where n could possibly be regarded as the input of sodium sources from Na⁺ through chemical molecule reactions (Cox and Plane, 1998), dynamic processes (Xu and Smith, 2003), or meteor injection (Richter and Sechrist, 1979). Without the non-linear effect, the input of the neutral sodium atoms (which could not be concentrated by electric field or wind shear like the ions) would diffuse away quickly. Fick’s law for the newly created particle flux into the background gas leads to the classical diffusion equation as follows (Schunk and Nagy, 2009): n x , t = N 2 π D t 1 / 2 e − x 2 / 4 D t , (3) where N is the newly generated particle number per unit area (e.g., column density); D = 1.3 k T m ν is the diffusion coefficient; k = 1.381 × 10 − 23 J / K is the Boltzmann constant; T = 180 K is the approximate kinetic temperature deduced from lidar observations; ν = 2 n v ¯ σ = 2 n 8 k T π m π d 2 is the averaged collision frequency between the sodium atoms and atmospheric particles; n = 10 13 m − 3 is the number density of the atmosphere according to results from SABER; m p = 1.66 × 10 − 27 k m and m N a S = 23 a m u ; and d = 2 × 186 × 10 − 12 m is the diameter of the sodium atom.

We assign the parameters with appropriate values as N = 3 × 10 13   m − 2 (assuming at the initial time t₀ the input sodium density n = 6000 c m − 3 and the width of Na_S equal to 5 km). Then, after one second, the input of 6000 c m − 3 would be attenuated to be n 0,1 = 1455 c m − 3 . Within 3.5 min, the input will reduce to nearly zero (e.g., less than 100cm ^-3, which is a minor amount compared with the normal sodium layer density).

After the input, n x , t would possibly undergo the evolution in two ways: a convergence or divergence. The full derivative form for n x , t can be expanded as follows: ∂ n ∂ t + v ∂ n ∂ x = 0 , (4) where d x ∂ t = v is the velocity of the particle movement. This equation has a generalized solution: n x , t = f x − v n t . (5)

This solution indicates each part of the wave has a different speed of v n . When d v d n > 0 , v n increases with an increase in density n and the wave packet front becomes steeper and steeper as the wave propagates, leading to a non-linear effect of convergence (Figure 1A).

On the other hand, a wave propagating in the x direction can be expressed as follows: u x , t = A e x p i k x − ω t , (6) where u is the wave function, A is the amplitude, k is the wave number, and ω is the angular frequency. The phase velocity v p and group velocity v g are given by v p = ω k (7) and v g = ∂ ω ∂ k , (8) respectively. If ω k ≠ c o n s t a n t , each wavelet will have a distinct velocity due to its individual wave vector. Hence, the term ω = ω k indicates a dispersion effect of the wave packet (Figure 1B).

Start with Euler’s equations for ideal fluid: ∂ u ∂ t = − 1 ρ ∂ p ′ ∂ x , ∂ w ∂ t = − 1 ρ ∂ p ′ ∂ z , ∂ u ∂ x + ∂ w ∂ z = 0 . (9)

The boundary conditions are w z = 0 = 0 , ∂ p ′ ∂ t − ρ g w z = h = 0 , (10) where u represents the horizontal velocity, w represents the vertical velocity, and p ′ is pressure.

Let L ≡ ∂ ∂ t ∂ 2 ∂ x 2 + ∂ 2 ∂ z 2 , then the formula for Eq. 9 can be compactly written as follows: L w = 0 . (11)

Let w = W z e i k x − ω t , (12) and we substitute Eq. 12 into Eq. 11: d 2 W d z 2 − k 2 W = 0 . (13)

Consequently, W z = A e k z + B e − k z , (14) where A and B are constants.

Substituting Eq. 14 into Eq. 11, we can get w = A e k z + B e − k z e i k x − ω t . (15)

Then, substituting Eq. 15 into Eq. 9, p ′ = i ω ρ k A e k z − B e − k z e i k x − ω t . (16)

According to the lower boundary conditions, B = − A , (17) and according to the upper boundary conditions, ω 2 k − g e k h − ω 2 k + g e − k h A = 0 . (18)

Thus, we can obtain ω = g k t a n h k h . (19)

Especially, in shallow water conditions, ω = g k k h − 1 3 k h 3 ≈ k 2 c 0 2 1 − 1 3 k 2 H 2 , = k c 0 1 − 1 6 k 2 H 2 = k c 0 − 1 6 k 3 c 0 H 2 , (20) where ω is a real number. Then, phase velocity v p and group velocity v g can be obtained, respectively, as follows: v p = ω k = c 0 − β k 2 , v g = ∂ ω ∂ k = c 0 − 3 β k 2 , (21) whereas β = c 0 H 2 / 6 . (22)

Therefore, when β ≠ 0 , v p ≠ v g , which fully indicates that the β ∂ 3 u ∂ x 3 term of the KdV equation characterizes the dispersion effect.

In addition, d v g d k = − 6 β k . (23)

Therefore, the effect of β ∂ 3 u ∂ x 3 causes wave dispersion. With the increase of β , the wavelength becomes shorter and the wave dispersion becomes stronger; such waves are known as dispersion waves. Of course, for waves with long wavelength (when k is small), it is a weakly dispersive wave, characterized by ω containing only the odd degree term of k .

On the other hand, since the layer density response is highly dependent on the density gradients occurring in the layer, the steady-state layer density profile becomes important (Gardner and Shelton, 1985). Large density gradients encourage non-linearities in the layer response (Gardner and Shelton, 1985). Therefore, when the non-linear effect affected by the gradient of the Na density profile is balanced with the dispersion effect mentioned previously, the wave shows neither dispersion nor non-linear characteristics, but propagates in the form of solitary waves.

So, the dispersion term of a surface wave in incompressible shallow fluid is given by ω k = g h k − 1 6 g h h 2 k 3 , (24) where h is the fluid depth and g is the gravitational acceleration (Belashov and Vladimirov, 2005).

In the complex space, we have ∂ ∂ t ↔ − i ω and ∂ ∂ x ↔ i k . Substituting into Eq. 24, we obtain ∂ u ∂ t + g h + u ∂ u ∂ x + 1 6 g h h 2 ∂ 3 u ∂ x 3 = 0 . (25)

We set u ′ = u + g h and β = 1 6 g h h 2 ; then, the equation is given by ∂ u ′ ∂ t + u ′ ∂ u ′ ∂ x + β ∂ 3 u ′ ∂ x 3 = 0 . (26)

This equation (similar to Eq. 1) is one of the simplest forms of the KdV equation, balanced by both the non-linear term u ∂ u ∂ x and dispersion term ∂ 3 u ∂ x 3 . Hereby, the solution satisfying this equation will undergo no convergence or dispersion effect, and the wave shape could be maintained for a long time.

The travelling wave could be represented by the form u x , t = f x − c t , where u x , t represents a disturbance moving in the negative or positive x-direction if c < 0 or c > 0 , respectively (Wazwaz, 2009). If the solution u x , t depends only on the difference between the two coordinates of the partial differential equations, then the solution keeps its exact shape and, therefore, is called a solitary wave. So, a solitary wave is a travelling wave whose transition from the asymptotic state at ξ = − ∞ to the other asymptotic state at ξ = ∞ is localized in ξ, where ξ = x-ct, and c is the wave speed (Wazwaz, 2009).

Eq. 26 or Eq. 1 has a special solution given by u x , t = u 2 + u 1 − u 2 sec ⁡ h 2 u 1 − u 2 12 β x − c t , (27) where u 1 and u 2 are exhibited by Figure 2A, and sech is referred to the hyperbolic secant function (Zabusky and Kruskal, 1965). This is just the bow wave observed by Russell in those early years, e.g., a solitary wave (Zabusky and Kruskal, 1965). Then, u 2 represents the limiting wave amplitude of Eq. 27 at infinity, and u 1 characterizes the peak of the wave. We set a = u 1 − u 2 and d = 12 β / a , where a is the amplitude and d is the width of the wave (shown by the vertical distance between the two red dashes of Figure 2A). Figure 2A could possibly corroborate some descriptions of solitary wave properties as follows: 1) this wave propagates along the x -direction, e.g., with the form of u x − c t ; 2) this wave is distributed in a limited space, e.g., lim x → ± ∞ u → 0 ; and 3) the shape of the wave does not change with time. This specific kind of non-linear wave is, therefore, called a soliton.

Eq. 27 is the soliton solution of Eq. 26 (similar to Eq. 1), and its evolution image is shown in Figure 3A when taking β as 1 and the parameter conditions are chosen as u 2 = − 1 and u 1 = 2 . The soliton’s shape is unaltered, and it advances uniformly in the positive x -direction. In the observed time series, the Na_S waveform of Figure 3B has a stable and prolonged propagation time, which is consistent with the temporal evolution of solitary waves. Additionally, the sodium density was only distributed over a narrow range of heights during the evolution of the Na_S event; i.e., at long distance from the peak height, the sodium density value is zero, indicating that this sort of Na_S event has properties similar to solitary waves. Some other patterns of isolated waves are affected by the non-linear and dispersion terms. In Section 3, the significance of non-linearity and dispersion in the evolution of Na_S will be discussed in relation to the Na_S occurrences recorded at the Andes lidar station.

It is shown that the sodium density variation with height can be approximated by the Gaussian distribution. n s z = C s 2 π σ s exp − z − z s 2 2 σ s 2 , (28) where n s z is the sodium number density at z , C s is the total column density of the sodium layer, z s is the centroid height, and σ s 2 is the RMS width (Xue, 2007).

When the Na_S occurs, the sodium density suddenly increases in a narrow and confined altitude range and the density profile obviously deviates from the Gaussian distribution. To quantitatively explore this anomaly, the Gaussian distribution function of the background sodium density should be determined first.

The sodium density data observed at one night by the lidar compose of a two-dimensional matrix: the elements on the column vectors of the matrix represent the sodium densities at different heights at a given moment, and the elements on the row vectors represent the results at different moments at a given height. Now, we choose the column vector of the matrix for data processing.

The origin observation data matrix of the day is set to be M , and then, M is represented by the column vector as follows: M = m → 1 , m → 2 , ⋯ , m → n − 1 , m → n , (29) where n is the number of observational points on the day. The maximum element in M is marked as d max . If we select an Na_S event with intensity factor>3, the critical value d_c for determining the anomalies is defined as follows: d c = d max 3 . (30)

When the value of an element in the column vector is greater than d_c, it reflects the anomaly of Na density. Otherwise, it is considered that the column vector conforms to the Gaussian distribution of Na density on that day. The column vector reflecting the Na density anomaly is arranged into a matrix in the order of the observed time as follows: A = a → 1 , a → 2 , ⋯ , a → k − 1 , a → k . (31)

Then, we set the matrix of column vectors matching the Gaussian distribution in the order of observational moments to be G = g → 1 , g → 2 , ⋯ , g → n − k − 1 , g → n − k . (32)

To obtain the Gaussian distribution of Na density on that day, we average all column vectors in G with a mark of g ⇀ a v e : g ⇀ a v e = ∑ i = 1 n − k g ⇀ i n − k , (33) where g ⇀ a v e reflects the distribution of Na density at the confined altitude. As mentioned previously, it is considered that g ⇀ a v e is consistent with the Gaussian distribution, i.e., the result of a Gaussian fit with g ⇀ a v e to the altitude h should satisfy Eq. 28. Taking the Na_S observed on 3 November 2016, for example, the Gaussian fitting is made with the column vector g ⇀ a v e and h ⇀ according to Eq. 28. The fitting results are as follows: g ⇀ h = 4330 exp − h ⇀ − 92.86 2 7.798 , (34) where g ⇀ h represents the modeled value that corresponds to the Gaussian distribution. The fitting result is shown as the red curve in Figure 2C.

Based on the observational data at Andes station from 20 August 2014 to 7 July 2019, 27 Na_S events are distinguished among 147 observation days. We continue to take the case on 3 November 2016 as an example. The column vector of d max is denoted as a ⇀ d   1 ≤ d ≤ k . The profile of a ⇀ d with h ⇀ is shown as the blue dash–dot line in Figure 2C. The observed Na density can be regarded as the sum of Gaussian distribution and the anomaly. Thus, vector p ⇀ reflecting the net anomaly of Na density is given by the following equation: p ⇀ = a ⇀ d − g ⇀ h , (35) and the distribution of p ⇀ with h ⇀ is shown as the blue dotted line in Figure 2D. Then, the peak density of p ⇀ equals to 10180.95 c m − 3 , with p ⇀ settled to about zero at infinity.

Then, the solitary wave fitting is performed on the vector p ⇀ . Since p ⇀ reflects the net anomaly at a determined time (the moment of the peak density profile), the parameter t in the fitted formula is a constant. At this time, the traveling wave ξ = x − c t represents a specific phase. In order to better represent the altitude corresponding to the peak Na density, Eq. 27 could be written as follows: u ξ = u 2 + u 1 − u 2 sec ⁡ h 2 u 1 − u 2 12 β ξ − ξ 0 , (36) where ξ 0 is the altitude corresponding to the maximum amplitude of the solitary wave. The fitting expression is finally deduced as follows: u ξ = − 82.17 + 10263.12 * sec ⁡ h 2 10263.12 19494.36 ξ − 93.63 , (37) which means u 2 = − 82.17 , u 1 − u 2 = 10263.12 (or u 1 = 10180.95 ), 12 β = 19494.36 , and

ξ 0 = 93.63 .

The fitting results are shown as the red curve in Figure 2D. Now, we know u 1 represents the peak density of the modeled curve, and u 2 corresponds to the density at infinity, i.e., u ξ ξ ≤ 80 ≈ u ξ ξ ≥ 110 ≈ lim ξ → ∞ u ξ → 0 . (38)

Let ξ 1 = ± d 2 = ± 3 β a = ± 3 β u 1 − u 2 . (39)

Substituting it into Eq. 37, we have u ξ 1 = u 2 + 4 e e + 1 2 u 1 − u 2 = u ± 0.6891 = 8071.41 , (40) which means the density at 0.6891 km from the peak is predicted to be 8071.41 c m − 3 . From Figure 2D, we can find all these simulated values; u 1 , u 2 and u ξ 1 are close to the observed results.

Furthermore, the height value at the theoretical full width at half maximum (FWHM) of a soliton is given as follows: h d / 2 = ξ 0 + d 2 = ξ 0 + 12 β u 1 − u 2 2   k m . (41)

To test the fitting results, the observed altitude of the selected case needs to be close to the calculated h 0.5 . Since the vertical resolution of the lidar is limited to 0.5 km, we can perform linear interpolation for the observed altitude series. The interpolated vector p ⇀ is denoted as p 1 ⇀ , and the altitude series of p 1 ⇀ is marked as h 1 ⇀ . The Na density at half-width is u h d / 2 ; then, it is always possible to find a value closest to u h d / 2 on both sides of the peak of p 1 ⇀ (denoted as p 1 m and p 1 n ). The corresponding altitudes for p 1 m and p 1 n are noted as h 1 m and h 1 n , respectively. Then, the observed width of the soliton could be given as follows: d ′ = h 1 n − h 1 m . (42)

According to the definition of the soliton, d ′ is also written as follows: d ′ = 12 β ′ max p 1 ⇀ − 0 , (43) where max p 1 ⇀ is the peak of p 1 ⇀ ; 0 is the density at infinity; and β ′ is the β value deduced from the observations. It is shown that β ′ = d ′ 2 max p 1 ⇀ 12 . (44)

The deduced values of h d / 2 , d , d ′ , u h d / 2 , p 1 , β , and β ′ , for each Na_S event, are listed in Table 1. The fitting parameters on 3 November 2016 are also exhibited in Figure 2D. Indeed, the results show the fitting parameters are close to the observed values.

Alternatively, the column vector obtained by Eq. 37 is denoted as q ⇀ ; then, q ⇀ is the predicted value at altitude h. The coefficient of determination ( R 2 < 1 ) and the root mean square error (RMSE) are then utilized to evaluate the fitting quality. They are calculated, respectively, as R 2 = 1 − ∑ i = 1 L q i − p i 2 ∑ i = 1 L p i − p ¯ i 2 (45) and R M S E = 1 L ∑ i = 1 L q i − p i 2 , (46) where L is the length of vector p ⇀ or q ⇀ , i.e., the number of observed altitude points.

The denominator of Eq. 45 denotes the residuals obtained by predicting the net anomaly vector p ⇀ , and the numerator represents the residuals predicted using the modeled vector q ⇀ .

The values of the parameters R 2 and RMSE obtained by fitting p ⇀ and q ⇀ for this Na_S event are calculated, respectively, as R 2 = 0.9784 and R M S E = 315.8402 , showing a good fit with R 2 close to 1. Through the evaluation of fitting quality, it is confirmed that Eq. 37 can fit the net anomaly vector p ⇀ within an allowable error (i.e., with R 2 > 0.9 ). The fitting results and quality evaluation parameters are obtained, with 24/27 events having R 2 > 0.9 as shown in Table 1. Therefore, the fitting results indicate u(ξ) are consistent with the net anomaly p ⇀ .

Figure 4 shows a simple data mining operation in Table 1. The net anomalous sodium density data p ⇀ were used for the fit, and the calculation of the fit results (as described in Eqs 35–44 produced two additional parameters β, d, which describe the soliton characteristics. In contrast, β ′ , d ′ are the corresponding parameters derived directly from the observed data. Assuming that the error bars’ lengths represent 10% of the parameter values, the ratios of β ′ , d ′ within the errors are 17/27 and 20/27, respectively, showing that we have had some success fitting the data.

In Figure 2D, it is shown that some small wavelets appear on both wings of the blue dotted line drawn from the net anomaly vector p ⇀ , which could not be explained by the standard solitary equation with the red fitting curve. These wavelets, however, have similar characteristics to the higher-order solitary wave with a waveform shown in Figure 2B.

Kawahara and Takuji proposed to add a higher-order dispersion term to the KdV equation, which considers dissipation, instability, and higher-order dispersion effects in the fluid medium (Kawahara and Takuji, 2007). This generalized KdV equation is written as follows: ∂ u ∂ t + u ∂ u ∂ x + β ∂ 3 u ∂ x 3 + γ ∂ 5 u ∂ x 5 = 0 , (47) which is also called the Kawahara equation (Kawahara and Takuji, 2007).

The numerical simulation results show that the Kawahara equation has two types of solutions: in the case of γ > 0 and β ≤ 0 , a soliton is formed with monotonic asymptotics similar to the soliton shown in Figure 2A; in the case of γ > 0 and β > 0 , the two wings of the soliton will have oscillation characteristics (Mamun and Shukla, 2009; Mamun and Shukla, 2009). Taking β = 1 and γ = 10 − 4 , the expression of Eq. 47 will be ∂ u ∂ t + u ∂ u ∂ x + ∂ 3 u ∂ x 3 + 10 − 4 ∂ 5 u ∂ x 5 = 0 , (48) whose solution at t = 0.0016 is shown in Figure 2B.

Furthermore, we can study the evolution of the waveform over time using the Fourier transform method. Let the initial conditions be β = 0.0023 , γ = 10 − 4 , and u x , 0 = e − 36 x 2 . Since the Kawahara equation contains two variables, x and t, the simulation results are denoted as a two-dimensional matrix W . Then, we arrange W by column vectors as follows: W = w ⇀ 1 , w ⇀ 2 , ⋯ , w ⇀ n − 1 , w ⇀ n . (49)

Since the vertical scale of Na_S observed by the lidar is generally less than 10 km, while the horizontal scale is often reported to be more than 300 km or even over 1,000 km (Fan et al., 2007; Ma et al., 2019), this scenario is inconsistent with the shallow water model. The solitary waves could particularly appear at the interface between upper and lower stratifications in the fluid medium (Bogucki and Garrett, 1993). They are frequently found in the stratified or sheared places in the ocean and atmosphere (Huthnance, 1989; Doviak et al., 1991; Gan and Ingram, 1992). In fact, in the lower atmosphere, in rotating, magnetized dusty plasma and in the dayside tropical mesosphere, many solitary wave events have also been observed (Christie et al., 1977; Christie et al., 1981; Mushtaq, 2006). This might suggest that solitary waves propagate beyond the mesospheric Na layer (the Na lidar range) with a larger vertical scale. Therefore, it is more effective to look for observational results related to stable stratifications or shears. The Na_S event on 9 April 2019 is selected as an example of this scenario (Figure 5A). This Na_S appears before the beginning of the observation at about 95 km altitude, with a downward propagation. The deduced vertical wind confirms the Na_S occurs in downward vertical wind (shown as the negative value region in Figure 5B). The vertical wind velocity field in this example has a low-velocity band with a width of around 3 km, and this low-velocity band (as shown in Figure 5B) decreases by roughly 9 km (95–86 km) in 10 hours (8 April 2019, 23:00–9 April 2019, 09:00), and the downward speed is roughly 0.9 km/h. Furthermore, the spatial and temporal distribution of Na_S on that day corresponds well with the distribution of this vertical wind low-velocity zone. The results of numerical simulations demonstrate that the non-linear term factor affects the propagation of the solitons along the spatial axis. Moreover, the downward vertical phase of the Na_S is consistent with the aforementioned properties of solitons. This implies that vertical wind speed has a non-linear function in this Na_S event, influencing the sodium density transfer. As stated in previous studies (Chen et al., 2021a; Chen et al., 2021b), Na_S is also modulated by the horizontal wind field. Figures 5C–E show the zonal wind, meridional wind, and the temperature profiles, respectively. These three profiles indicate the Na_S has been located near the stratification with strong shear. Furthermore, the stability of the stratification is determined by the Richardson number R i (Nappo, 2002; Liu and Liu, 2011). Under stable stratification, both convection and turbulence are less likely to develop.

In unit time, due to vertical displacement, for a unit mass of air, the convective flow energy to resist the net Archimedes buoyancy is as follows: W 1 = − d d t 1 2 w 2 ¯ = K H N 2 , (50) where w = dz/dt is the vertical speed; K H = w z ¯ > 0 is the convective conductivity coefficient; N 2 = g T ∂ T ∂ z + g C p is the buoyancy frequency; g = 9.5   m s − 2 is the gravitational acceleration in the mesopause; C p = 1004   J K − 1 k g − 1 is the specific heat at constant pressure; and T is the atmospheric thermodynamic temperature.

On the other hand, for unit mass of air, the horizontal kinetic energy consumed per unit time in the presence of vertical wind shear, i.e., the kinetic energy provided to convective motion, is as follows: W 2 = − d d t 1 2 V ⇀ h 2 ¯ = K M ∂ u ∂ z 2 + ∂ v ∂ z 2 , (51) where u and v are the zonal and meridional wind velocities, respectively. K M > 0 is the momentum transport coefficient and is usually approximately equal to K H .

Thus, the stability of the layer depends on the value of W 1 / W 2 . When W 2 is less than W 1 , it means that the disturbance kinetic energy converted by the basic airflow is less than the disturbance kinetic energy consumed by stable stratification. In this situation, even if convection or turbulence occurs, it will be suppressed or weakened so that the atmosphere is in a stable state. The dimensionless ratio, W 1 / W 2 , is defined as follows: R i = W 1 W 2 = g T ∂ T ∂ z + g C p ∂ u ∂ z 2 + ∂ v ∂ z 2 , (52) which could be deduced from the wind and temperature results observed by the lidar (Figures 5C–E). The calculated Ri with a value > 1.5 are shown as black scatters in Figure 5F. The Na_S locates around the area where scatters are concentrated, i.e., a special stable area near 95 km. It is obvious that the Na_S and stable region evolved synchronously, and finally, both become blurred. Therefore, the lidar observations and the deduced results are all consistent with the appearance of a solitary wave. A stable stratification with R i   > 1.5 exists nearly about 95 km, accompanied by the zonal wind, meridional wind, and temperature being also stratified around 95 km. Once a fluctuation is excited in the vicinity of the stratification and just satisfies the balance of non-linear and dispersion effects, the waveform will maintain through the propagation; i.e., a solitary wave appears.

Furthermore, the evolution of the observed net anomaly for the Na_s throughout whole night is also needed for comparison. In Section 3.1, the initial observation of sodium density is denoted as M = m → 1 , m → 2 , ⋯ m → n − 1 , m → n .

The Gaussian distribution model of the current day is noted as a two-dimensional matrix G ⃡ 0 , and we define G ⃡ 0 as isomorphic to M ⃡ . Then, G ⃡ 0 could be written in the form of a column vector as follows: G 0 = g → h , g → h , ⋯ , g → h , g → h . (53)

Obviously, the net anomaly evolution of the Na_S is obtained by observation data matrix M ⃡ subtracting the Gaussian distribution model G ⃡ 0 , which could be denoted as S ⃡ . There are S = M − G 0 = m ⇀ 1 − g ⇀ h , m ⇀ 2 − g ⇀ h , … , m ⇀ n − 1 − g ⇀ h , m ⇀ n − g ⇀ h = s ⇀ 1 , s ⇀ 2 , … , s ⇀ n − 1 , s ⇀ n , (54) where s ⇀ i ( i = 1 , 2 , … , n − 1 , n ) is the column vector of S ⃡ .

In order to compare the numerical simulation results W ⃡ with the Na_S evolution model S ⃡ more intuitively, the Fourier transform method is used uniformly to disperse a certain period of time on time variable t into n moments in numerical simulation. Thus, both W ⃡ and S ⃡ are two-dimensional matrices with the same column number n. To better show them, a dynamic video of the variation of their column vectors with time is uploaded as Supplementary Movie S1. The five images on the left in Figure 6 (i.e., Figures 6A, C, E, G, I) are single-frame captures from W ⃡ of the video, showing the evolutions of the five-order solitary wave over time. The five images on the right (i.e., Figures 6B, D, F, H, J) are intercepted from S ⃡ , indicating variations of the net anomaly of the Na_S. Figures 6A, B show at this moment, the simulated wave shapes are similar to the observed peak density profile. Figures 6C, D show the huge peaks attenuate gradually. Figures 6E, Fshow the peaks decay to about zero value. Figures 6G, H show the peaks change phase and resume. Figures 6I, J show the peaks recover to a sharp form similar to the initial condition, except with different phases. By comparison, it is verified that the numerical simulation results W ⃡ are in good agreement with the evolution process S of the Na_S. So, the five-order solitary wave theory is a potential candidate in explaining this Na_S event.

However, it is worth noting that the numerical simulation of the higher-order KdV equation is probably only suitable for explaining the events similar to the selected case. These events are typically characterized by occurrence heights below 95 km, longer durations, and descending patterns similar to tidal fluctuations. In contrast, the other events with shorter durations and cloud-like shapes are less consistent with the higher-order simulation results. This discrepancy also implies that Na_S with different characteristics may have different fine structures.

The fluctuations can affect various physical parameters in the upper mesosphere. During the Na_S event on 9 April 2019, the temperature in the mesosphere displayed a distinct wave signature similar to the solitary wave observed in sodium density.

Similar to the approach used to analyze sodium density data, in the current study, the background temperature distribution (illustrated in Figure 7A) was obtained by averaging the temperature data of April 2019, excluding the temperature data recorded during the Na_S event. The temperature anomaly distribution (illustrated in Figure 7B) was obtained by differencing the temperature observations on 9 April 2019 when the Nas event occurred with the background temperature.

Figure 7A displays the background temperature distribution in the 93–97 km region consistently maintaining 190–200 K from −0.8 UT to 2.8 UT (where 9 April 2019 00:00 UT is marked as 0 UT), with no significant anomaly. The distribution of temperature anomalies from 93 to 97 km above the Andes lidar station on 9 April 2019 is shown in Figure 7B. Compared with the background temperature profile, there is a stripe or band of warming region with an amplitude of approximately 15 K within 93–97 km from −0.8 UT to 2.8 UT, overlapping the region of Na_S. Figure 8 presents the temperature plots (black lines) vs. the background temperature profiles (red lines) for 93–97 km, with (a) to (f) corresponding to the time sequence from −0.1 UT, 0.1 UT, 0.3 UT, 0.5 UT, and 0.7 UT, to 0.9 UT. Figures 8A–F show that the temperature anomalies in this region persist at approximately 10–15 K throughout the Na_S event (−0.8 UT∼2.8 UT), indicating a stable propagation pattern/signal similar to that expected for solitary wave.

Moreover, in contrast to Figure 2C, the temperature profile does not undergo any significant changes as the sodium density varies during the occurrence of the Na_S event. This discrepancy could potentially be explained by the following factors: 1) The occurrence of Na_S involves an input source (the sodium ions from E_S could potentially offer a sufficient neutral sodium atom source through recombination with free electrons). This implies that passive temperature variations may not be as sensitive to fluctuations as Na_S. 2) The temperature in the mesopause is also affected by atmospheric dynamical processes such as solar radiation, photochemical reactions, gravity waves, and tidal waves. Therefore, accurately determining the background temperature becomes challenging, which can result in less-pronounced fluctuation features in the temperature profile.