© The Authors 2024

1 Introduction

Acoustic-gravity waves (AGWs) are omnipresent features of the stratified atmosphere and have been extensively studied on various planets and moons, such as Earth (e.g., Yeh & Liu 1974; Yiǧit et al. 2008, 2009, 2012, 2014, 2016, 2021a; Yiǧit & Medvedev 2012, 2015, 2016, 2010, 2017; Medvedev & Yiǧit 2019; Miyoshi & Yiǧit 2019; Nayak, & Yiǧit 2019; Medvedev et al. 2023), Venus (e.g., Piccialli et al. 2014; Silva et al. 2024), Jupiter (e.g., Young et al. 2005; Lian & Yelle 2019), Saturn (e.g., Müller-Wodarg et al. 2019; Brown et al. 2022), Titan (e.g., Cui et al. 2013, 2014; Wang et al. 2020; Huang et al. 2022b, 2023), and exoplanets (e.g., Watkins & Cho 2010; Rogers et al. 2012; Peralta et al. 2014). AGWs have also been observed on Mars via remote and in situ measurements over the past several decades (e.g., Fritts et al. 2006; Creasey et al. 2006; Yiǧit et al. 2015a,b, 2018, 2021b; Medvedev et al. 2016; England et al. 2017; Terada et al. 2017; Siddle et al. 2019; Nakagawa et al. 2020; Heavens et al. 2020; Starichenko et al. 2021, 2024; Li et al. 2021; Saunders et al. 2021; Ji et al. 2022; Leelavathi & Rao 2024; Pankine et al. 2024). AGWs are formed when air parcels in the stratified atmosphere are displaced from their equilibrium positions by small perturbations, with buoyancy and gravity as the restoring forces, resulting in oscillation of the atmosphere at a characteristic frequency called the buoyancy frequency or the Brunt-Väisälä frequency (Vaisala 1925; Brunt 1927). When the wave frequency is lower than the buoyancy frequency, the potential energy associated with the buoyancy force is dominated. In this scenario, the wave tends to become an internal gravity wave. When the wave frequency is much greater than the buoyancy frequency, the gravitational effect is negligible. In this case, the wave tends to become an ordinary sound wave (e.g., Yeh & Liu 1972). Generally, the period of AGWs ranges from a few minutes to tens of minutes.

Acoustic-gravity waves can be triggered at lower heights. Topographical winds, instability of weather systems, and convection near the Martian surface are potential wave sources (e.g., Parish et al. 2009; Heavens et al. 2020; Yiǧit et al. 2021b; Huang et al. 2022a; Wu et al. 2022; Shaposhnikov et al. 2022). The AGWs generated in the lower atmosphere can propagate to the upper atmosphere. Wave amplitudes can exponentially grow with height due to decreasing neutral density. The vertically propagating AGWs can be dissipated mainly by molecular diffusion and thermal conduction in the upper regions, causing momentum and energy carried by waves to be released into the background atmosphere (Hines 1960; Vadas & Fritts 2005; Yiǧit & Medvedev 2019). Consequently, AGWs can result in substantial thermal variations of the upper atmosphere and eventually contribute to water loss and atmospheric escape on Mars (Medvedev & Yiǧit 2012; Walterscheid et al. 2013; Medvedev et al. 2015; Yiǧit et al. 2016, 2021; Yiǧit 2021, 2023; Roeten et al. 2022; Kuroda et al. 2015, 2020).

Acoustic-gravity waves can also impact the Martian ionosphere. AGWs originating from lower regions can also propagate into ionospheric heights and interact with plasma through neutral–ion and neutral–electron collisions, directly or indirectly creating various structures in the ionosphere. Recently, the plasma variations driven by AGWs in the thermosphere–ionospheric region on Mars were directly measured by the Neutral Gas and Ion Mass Spectrometer (NGIMS) on board the Mars Atmosphere and Volatile EvolutioN (MAVEN) mission (Leelavathi et al. 2023). Highly variable ionospheric structures at Mars have also been observed via remote and in situ measurements. Gurnett et al. (2010) and Harada et al. (2018) observed ionospheric irregularities in the Martian ionosphere by analyzing the sounding data from Mars EXpress (MEX). MAVEN also observed ionospheric irregularities on Mars (Fowler et al. 2017, 2019, 2020; Tian et al. 2022; Wan et al. 2024a,b). The results of these latter authors suggest that AGWs arising from the lower atmosphere are seeding mechanisms of the ionospheric irregularities due to their comparable wavelengths and amplitudes, and broad temporal, spatial, and geographic distributions. Collinson et al. (2019) first observed extraterrestrial traveling ionospheric disturbances (TIDs) in the Martian ionosphere by MAVEN and argued that wavelike oscillations in plasma density and periodic fluctuations in magnetic field associated with the TIDs are mostly driven by AGWs. Sporadic E-like layers or rifts (Collinson et al. 2020), plasma depletions or plasma bubble-like structures (Basuvaraj et al. 2022), and magnetic depletions (Fowler et al. 2022) in the ionosphere of Mars are also observed by MAVEN. The possible formation mechanisms of these events include electromagnetic gradient drift instability (e.g., Ecklund et al. 1981; Woodman et al. 1991), shear-driven Kelvin-Helmholtz instability (e.g., Bernhardt 2002; Penz et al. 2004; Ruhunusiri et al. 2016; Poh et al. 2021; Wang et al. 2022), and gravity-driven Rayleigh-Taylor instability (e.g., Ott 1978). AGWs are believed to trigger these instabilities in the ionosphere (e.g., Kelley et al. 1981; Huang et al. 1993, 1994; Huang & Kelley 1996a,b,c; Singh et al. 1997; Taori et al. 2011; Didebulidze et al. 2021; Zawdie et al. 2022). Therefore, AGWs play an important role in the neutral atmosphere–ionospheric coupling processes, which is an important issue to be explored thoroughly on Mars.

Up to the present, the wave-driven variable structures in Earth’s ionosphere have been extensively studied, such as the Earth’s ionospheric irregularities (e.g., Hines 1960; Hooke 1968), TIDs (e.g., Hooke 1970b,a,a, 1971; Hunsucker 1982; Hocke & Schlegel 1996; Balthazor & Moffett 1997; MacDougall et al. 2009; Azeem et al. 2017; Miyoshi et al. 2018), sporadic-E layers (e.g., Kirkwood & Collis 1989; Woodman et al. 1991), and plasma bubbles (e.g., Singh et al. 1997; Taori et al. 2011; Takahashi et al. 2009; Abdu et al. 2009; Paulino et al. 2011; Ajith et al. 2020). The plasma oscillations driven by AGWs on the ionosphere of other celestial bodies have also been investigated, such as those on Jupiter (Matcheva et al. 2001; Barrow & Matcheva 2011), Saturn (Matcheva & Barrow 2012; Barrow & Matcheva 2013) and Titan (Wu et al. 2022). For Mars, although the effects of AGWs on the neutral atmosphere have been simulated in earlier works (e.g., Parish et al. 2009; Medvedev et al. 2011, 2013; Walterscheid et al. 2013; Imamura et al. 2016; Kuroda et al. 2016, 2019; Srivastava et al. 2022; Roeten et al. 2022; Liu et al. 2023), the coupling processes between AGWs and the ionosphere are poorly explored. Martian ionospheric irregularities have been studied by theoretical (Keskinen 2018) and numerical methods (Jiang et al. 2021, 2022), but the impacts of AGWs on ionospheric irregularities are not taken into consideration in the previous works. Wang et al. (2023) proposed a linearized wave model adopting a Wentzel-Kramers-Brillouin (WKB) approximate method, which has been used in the Jovian and Saturnian ionosphere (e.g., Matcheva et al. 2001; Barrow & Matcheva 2013), only to study the wave–electron interaction in the dayside ionosphere of Mars. However, the wave-induced ion variations in the presence of photochemical reactions were not considered.

The aim of the present study is to examine the coupling processes between AGWs and the Martian ionosphere by constructing a comprehensive wave model. This model can be used to quantitatively investigate the wave-induced variations in the ion/electron density, velocity, and temperature combined with photochemical reactions. In addition, the perturbations in the electromagnetic field driven by AGWs can also be calculated via this model. Unlike Earth, Mars presently lacks a strong global dipole magnetic field (Acuna et al. 1998). However, there is a crustal magnetic field located in the southern hemisphere. Beyond the crustal magnetic field region, the magnetic fields are dominated by the draped interplanetary magnetic field, whose magnitude is tens of nT. As a result, the beta of plasma, the ratio of plasma pressure to the magnetic pressure, is close to 1 in the Martian ionosphere, while this value is much less than 1 at Earth. The apparent wavelike oscillations in the ambient magnetic field associated with the ionospheric disturbances are observed on Mars (Collinson et al. 2019; Fowler et al. 2019). This unique feature clearly demonstrates that the magnetic field in the extraterrestrial ionosphere could be more easily influenced by plasma behavior if AGWs are propagated into the ionospheric heights. We therefore also investigated the wave-driven thermal variations (wave heating/cooling effect) in the Martian ionosphere–thermospheric region.

This paper presents a detailed study of the ionospheric response to AGWs. We focus on the following questions. First of all, we are interested in how AGWs impact the motions of electrons and ions. Second, we would like to identify the differences between the wave-driven plasma density variations in the absence of photochemical reactions and those in the presence of photochemistry (equivalently, in the dayside or night-side of the Martian ionosphere). In addition, we examine the wave-driven variations in the electromagnetic field and thermal structure (wave heating and cooling effect) in the ionosphere-thermospheric region.

The layout of this paper is as follows. The model details are presented in Sect. 2. We adopt a comprehensive coupled model to describe the wave propagation in the ionosphere-thermospheric region and to calculate the wave perturbations associated with the dissipative processes (including molecular viscosity and thermal conductivity). In addition, we study the coupling processes between AGWs and plasma. The model results, including the wave-induced plasma density, velocity, pressure, and temperature variations, the wave-driven electromagnetic field, and the wave-associated thermal variations, are displayed in Sect. 3. We then present a comparison between our model results and spacecraft observations in Sect. 4. Finally, we summarize and discuss our findings and draw conclusions in Sect. 5.

2 Model description

We describe the model to investigate the coupling processes between AGWs and the ionosphere in this section. Sect. 2.1 shows the governing equations describing the behaviors of neutral species and plasma. The corresponding assumptions are also given. A coupling linear system is described in detail in Sect. 2.2. The neutral wave solutions associated with boundary conditions are shown in Sect. 2.3.

2.1 Governing equations and assumptions

We investigate the wave behaviors in the ionosphere-thermospheric region. The governing equations include two parts. The first part of a system of equations described the behavior of neutral atmosphere is displayed in Sect. 2.1.1. The second part involves the ionospheric plasma behaviors associated with photochemistry and electromagnetic field and is shown in Sect. 2.1.2.

2.1.1 Neutral atmosphere

To characterize the wave propagation within a realistic atmosphere, we adopt a linear wave model, called a “full-wave” model, as the first part of our coupled model. The full-wave model has been proposed to describe the wave propagation within a dissipative, compressible, and inhomogeneous atmosphere. The term “full-wave” means that this model fully takes into account the effect of wave reflection, Coriolis force, eddy diffusion, ion-drag effect, and the variation in the background temperature with height. This model can describe wave propagation in the lower, middle, and upper atmosphere. Eddy turbulence is a predominant controlling factor for the composition of the lower and middle atmosphere. It is usual to describe the turbulence by using the concept of eddy diffusion (e.g., Lettau 1951; Colegrove et al. 1965, 1966). Without lose generality, the eddy effect is considered in our model. The full-wave model has been applied on both Earth (e.g., Volland 1969b,a) and Jupiter (e.g., Hickey et al. 2000). We modify the published linear full-wave model by adding neutral-plasma collisional effects and photochemical heating to describe the wave dissipation in the ionosphere. Therefore, the equations governing the behavior of the neutral atmosphere in our revised full-wave model include mass continuity equation, momentum equation, energy equation, and ideal gas law, and can be expressed as the following:

$$\frac{\partial {\rho }_{\text{n}}}{\partial t}+\left(V_{\text{n}}\cdot \nabla \right){\rho }_{\text{n}}+{\rho }_{\text{n}}\nabla \cdot V_{\text{n}}=0$$(1a)

$${\rho }_{\text{n}}\frac{\partial V_{\text{n}}}{\partial t}+{\rho }_{\text{n}}\left(V_{\text{n}}\cdot \nabla \right)V_{\text{n}}+\nabla p_{\text{n}}-{\rho }_{\text{n}}g=F_{\text{n}},$$(1b)

$$\frac{\partial s}{\partial t}+\left(V_{\text{n}}\cdot \nabla \right)s=\frac{Q_{\text{n}}}{T_{\text{n}}},$$(1c)

and

$$p_{\text{n}}={\rho }_{\text{n}}R{T}_{\text{n}}\text{, }$$(1d)

respectively, where ρ_n, p_n, and T_n represent the mass density, pressure, and temperature of the neutral atmosphere, respectively; g is the gravitational acceleration; and subscript “n” denotes neutral species.

The background profiles of the neutral atmosphere, including the density, temperature, and pressure, are assumed to only vary with altitude. The vertical structures of the background atmosphere in the absence of AGWs are displayed in Fig. 1. The major neutral compositions of the Martian atmosphere are composed of CO₂, CO, N₂, and O, whose number densities are plotted in Fig. 1a and marked by $n_{{\text{CO}}_2},n_{\text{CO}},n_{{\text{N}}_2}$, and n_O, respectively. The neutral temperature is plotted in Fig. 1c and marked by T_n. These smooth background neutral density and temperature profiles are measured from Deep-dip 2 Campaign of MAVEN spacecraft and adopted from Wu et al. (2020). The neutral pressure, p_n, can be calculated by Eq. (1d) and plotted in Fig. 1d.

V_n in Eqs. (1a)–(1c) denotes the neutral velocity and can be given by V_n = (u, υ, w) in a local geographic coordinate system, O – xyɀ. As shown in Fig. 2, the positive x-axis in this coordinate system is toward the east, the positive y-axis points toward the north, and the positive z-axis points upwards. u, υ, and w represent the zonal, meridional, and vertical wind speeds, respectively.

F_n in Eq. (1b) represents the force per unit volume and is given by:

$$\begin{array}{l}{F}_{\text{n}}=-2{\rho }_{\text{n}}\Omega \times V_{\text{n}}+\nabla \cdot {\sigma }_{\text{m}}\\ +\nabla \cdot \left({\mu }_{\text{eddy}}\nabla V_{\text{n}}\right)-K_{\text{R}}{\rho }_{\text{n}}{V}_{\text{n}}-f_{\text{nc}}.\end{array}$$(2)

The terms on the right-hand side represent: (a) Coriolis force; (b) molecular viscous force; (c) eddy viscous force; (d) Rayleigh friction; and (e) neutral-plasma collisional force. Ω is the planetary angular velocity; µ_eddy represents the eddy momentum diffusivity; K_R is the Rayleigh friction coefficient.

Considering AGWs are a series of small-scale waves, their wavelengths are smaller than the planetary radius. An ‘f-plane’ approximation is applied to this model. Under this consideration, the Coriolis force in the geographic coordinate system can be expressed as:

$$F_{\text{c}}=-2{\rho }_{\text{n}}\Omega \times V_{\text{n}}=(-vf,uf,0),$$(3)

where f = 2Ω sin ψ with ψ being geographic latitude.

σ_m in Eq. (2) is the viscous stress tensor and can be represented by:

$$\begin{array}{l}{\sigma }_{\text{m}}={\mu }_{\text{m}}\nabla V_{\text{n}}+{\mu }_{\text{m}}{\left(\nabla V_{\text{n}}\right)}^T\\ -\frac{2}{3}{\mu }_{\text{m}}\mathrm{diag}\left\{\nabla \cdot V_{\text{n}},\nabla \cdot V_{\text{n}},\nabla \cdot V_{\text{n}}\right\},\end{array}$$(4)

where the superscript “T” indicates the transform of a matrix, and µ_m is the molecular dynamic viscosity.

As a result, the molecular viscous force per unit volume can be expressed as:

$$F_{\text{v}}=\nabla \cdot {\sigma }_{\text{m}}={\mu }_{\text{m}}{\nabla }^2{V}_{\text{n}}+\frac{{\mu }_{\text{m}}}{3}\nabla \left(\nabla \cdot V_{\text{n}}\right)$$(5)

This model includes the collisional effects between neutrals and plasmas in the ionosphere-thermospheric region. ƒ_nc in the Eq. (2) is the collisional force between neutrals and plasmas:

$$f_{\text{nc}}=m_{\text{i}}{n}_{\text{i}}{v}_{\text{in}}\left(V_{\text{n}}-V_{\text{i}}\right)+m_{\text{e}}{n}_{\text{e}}{v}_{\text{en}}\left(V_{\text{n}}-V_{\text{e}}\right),$$(6)

where subscripts ‘n’, ‘i’ and ‘e’ represent neutral, ion, and electron, respectively; m denotes mass and n is number density; ν_in and ν_en are the ion-neutral and electron-neutral collisional frequency, respectively.

Equation (1c) describes energy conservation in the neutral atmosphere. s in this equation represents atmospheric entropy: s = c_p ln ϑ = c_pln T_n − R ln p_n + Rln p_s with c_p being the specific heats at constant pressure and R = c_p − c_v being the specific gas constant; c_v is the specific heat at constant volume; $T_{\text{n}}{\left(p_{\text{s}}/p_{\text{n}}\right)}^{R/c_{\text{p}}}$ represents the potential temperature that an air parcel would attain if it were brought down adiabatically from a region with pressure p_n to a height where the pressure is p_s. In a realistic atmosphere, the wave propagation processes are not adi-abatic, indicating the presence of an external heat source or sink. The heat input per unit mass on the right-hand side of Eq. (1c), Q_n, is given by (e.g., Hickey et al. 1998, 2022; Schubert et al. 2003):

$$\begin{array}{l}{Q}_{\text{n}}=\frac{1}{{\rho }_{\text{n}}}[{\sigma }_{\text{m}}:\nabla V+\nabla \cdot \left({\lambda }_{\text{m}}\nabla T_{\text{n}}\right)+\frac{c_{\text{p}}{T}_{\text{n}}}{\vartheta }\nabla \cdot \left({\rho }_{\text{n}}{\kappa }_{\text{eddy}}\nabla \vartheta \right)\\ -{\rho }_{\text{n}}{c}_{\text{v}}{K}_{\text{N}}{T}_{\text{n}}+Q_{\text{ni}}+Q_{\text{ne}}+Q_{\text{pc}}].\end{array}$$(7)

The terms on the right-hand side represent (a) molecular viscous heating; (b) thermal conduction heating; (c) eddy viscosity heating; (d) Newtonian cooling heating; (e) neutral–ion collisional heating; (f) electron-neutral collisional heating; and (g) photochemical heating of hot oxygen atom. λ_m represents the thermal conductivity; κ_eddy is the eddy thermal diffusivity; K_N denotes the Newtonian cooling coefficient.

The presence of energy exchange between neutrals and plasmas in the ionosphere-thermospheric region suggests that neutral-plasma collisional heating is considered in this model. The collisional heating between neutrals and plasmas can be expressed as (e.g., Huba et al. 2000; Niu et al. 2021):

$$\begin{array}{l}{Q}_{\text{n}\alpha }=\frac{{\rho }_{\text{n}}{\tilde{v}}_{\text{n}\alpha }{k}_{\text{B}}}{m_{\text{n}}}\left(T_{\alpha }-T_{\text{n}}\right)\\ +\frac{1}{2}\frac{m_{\alpha }}{m_{\text{n}}}{\rho }_{\text{n}}{v}_{\text{n}\alpha }\left(V_{\alpha }-V_{\text{n}}\right)\cdot \left(V_{\alpha }-V_{\text{n}}\right),\end{array}$$(8)

where k_B is the Boltzmann constant; subscripts ‘n’ and ‘α’ represent neutral and plasma, respectively; and ${\tilde{\nu }}_{\text{n}\alpha }$ is defined by ${\tilde{\nu }}_{\text{n}\alpha }=3m_{\alpha }{m}_{\text{n}}{v}_{\text{n}\alpha }/{\left(m_{\alpha }+m_{\text{n}}\right)}^2$ with m being the mass and ν_nα being the neutral-plasma collisional frequency.

Photochemical heating is considered in this study. The dissociative recombination of ${\text{O}}_2^{+},{\text{R}}_1:{\text{O}}_2^{+}+{\text{e}}^{-}\to \text{O}+\text{O}$, occurs in the Martian ionosphere. Two energetic oxygen atoms are produced via the reaction between ${\text{O}}_2^{+}$ and ambient electrons. The detailed reaction of R₁ proceeds via 4 different channels and is shown below:

$$\begin{array}{ccc}{\text{O}}_2^{+}+{\text{e}}^{-}\to & \text{O}\left({}^3\text{P}\right)+\text{O}\left({}^3\text{P}\right)+6.99\text{eV}& 26.5\%\hspace{0.17em}\text{(a)},\\ & \text{O}\left({}^1\text{D}\right)+\text{O}\left({}^3\text{P}\right)+5.02\text{eV}& 47.3\%\hspace{0.17em}\text{(b)},\\ & \text{O}\left({}^1\text{D}\right)+\text{O}\left({}^1\text{D}\right)+3.06\text{eV}& 20.4\%\hspace{0.17em}\text{(c)},\\ & \text{O}\left({}^1\text{D}\right)+\text{O}\left({}^1\text{S}\right)+0.84\text{eV}& 5.8\%\hspace{0.17em}\text{(d)},\end{array}$$(9)

where O(³P), O(¹D), and O(¹S) are the ground and excited electronic states of atomic O; the numbers in units of eV are the associated exothermicities and the percentages refer to the branching ratios of different channels (e.g., Lillis et al. 2017).

The photochemical heating of the hot oxygen atom, Q_pc, mentioned on the right-hand side of Eq. (7) can be calculated by $Q_{\text{pc}}=R_{\text{O}}{\displaystyle \sum _j{E}_j}{L}_j$, where subscript ‘j’ denotes the jth channel of reaction mentioned by Eq. (9); E_j and L_j are the associated exothermicity and the branching ratio of the jth channel; R_O is the production rate of the hot atomic O and is given by $R_{\text{O}}=n_{\text{e}}{n}_{{\text{O}}_2^{+}}{\alpha }_1$; here $n_{\text{e}},n_{{\text{O}}_2^{+}}$, and α₁ are the undisturbed number density of electron and ${\text{O}}_2^{+}$, and the rate coefficient of reaction R₁ , respectively.

The governing equations mentioned above describe the wave propagation in the upper atmosphere. The atmosphere is assumed compressible, suggesting that the acoustic waves have not yet been ruled out. The local sound speed can be calculated by $c_{\text{S}}=\sqrt{\gamma R{T}_{\text{n}}}$ and displayed in Fig. 3a, where γ = c_p/c_v is the rate of specific heat. The acoustic cut-off frequency, ω_a, is defined as: ω_a = γ𝑔/2c_s. Correspondingly, the acoustic cut-off period is τ_a = 2π/ω_a and is plotted in Fig. 3b. The air parcels in the atmosphere can be displayed from their equilibrium positions, with buoyancy and gravity as the restoring force. Under this situation, the atmosphere oscillates at the Brunt-Väisälä frequency. The definition of the Brunt-Väisälä frequency, ω_b, can be given by (Nappo 2002):

$${\omega }_{\text{b}}^2={\Gamma }_{\text{ad}}\frac{\partial s}{\partial z}=\frac{g\left({\Gamma }_{\text{ad}}-\Gamma \right)}{T_{\text{n}}},$$(10)

where Γ_AD = 𝑔/c_p is the adiabatic lapse rate and Γ = −∂T_n/∂ɀ is the negative temperature gradient. ${\omega }_{\text{b}}^2>0\text{ if }\Gamma <{\Gamma }_{\text{ad}}$ if Γ < Γ_ad. The atmosphere in this situation is typically stable, causing air parcels to oscillate at the Brunt-Väisälä frequency. The Brunt-Väisälä period, τ_b = 2π/ω_b, is displayed in Fig. 3b. The region where τ > τ_b domained by internal gravity wave branches while acoustic wave branches domain the region where τ < τ_a, as displayed in Fig. 3b. The region where τ_a < τ < τ_b is domained by evanescent waves or surface waves, which can only exist at the atmospheric boundary. The average values of τ_a and τ_b in the Martian atmosphere are ~10 min and ~15 min, respectively.

The presence of molecular viscosity and thermal heat conduction causes wave attenuation in the upper atmosphere. The molecular dynamic viscosity, µ_m, and thermal conductivity, λ_m, can be given by empirical equations:

$${\mu }_{\text{m}}=\frac{{\displaystyle \sum A_k}{n}_k}{{\displaystyle \sum n_k}}{T}_{\text{n}}^{0.69},{\lambda }_{\text{m}}=\frac{{\displaystyle \sum C_k}{n}_k}{{\displaystyle \sum n_k}}{T}_{\text{n}}^{0.69}$$(11)

where A_k and C_k are the numerical coefficient and can be taken from Rees (1989); n_k is the number density of kth neutral constituent. The molecular kinetic viscosity, ν_m = µ_m/ρ_n, increases exponentially with altitude, as shown in Fig. 3c. Similarly, the thermal diffusivity is given by κ_m = λ_m/c_pρ_n and displayed in Fig. 3c. The eddy effect is unimportant in the thermosphere compared with the molecular viscosity and thermal conduction. We ignore the eddy effect in this study since we only focus on the wave propagation in the ionosphere-thermospheric region.

The altitude in our model ranges from 80 km to 300 km. A sponge layer applied at the upper boundary of the model to rule out the effect of spurious wave reflections from the upper boundary (e.g., Hickey et al. 2000; Walterscheid et al. 2013; Srivastava et al. 2022). Newtonian cooling and Rayleigh friction coefficient, mentioned in Eqs. (2) and (7), respectively, are used to implement the upper sponge layer. The Rayleigh friction, associated with zonal wind drag force, was introduced to explain the anomalous temperature distribution in the upper mesosphere and lower thermosphere (Matsuno 1982; Holton 1982). Earlier studies indicated that wave amplitude growth could be influenced by the Rayleigh friction (e.g., Matsuno 1982; Holton 1982; Smith & Lyjak 1985). The Martian atmosphere is primarily dominated by CO₂ constituent. Therefore, radiative damping due to CO₂15 µm band is significant in the Martian atmosphere. The CO₂ infrared radiation can cool the background temperature and should be included in the realistic modeling of AGWs (e.g., Eckermann et al. 2011; Medvedev et al. 2015; Srivastava et al. 2022). The CO2 infrared radiative cooling can be approximated by Newtonian cooling (Lindzen & Goody 1965; Leovy 1964a,b). The Rayleigh friction and Newtonian cooling effects are applied to the sponge layer at the upper boundary in our model, following previous full-wave model given by Hickey et al. (2000) and Walterscheid et al. (2013). The coefficients K_N and K_R are assumed numerically equal in the full-wave model and can be given by (Walterscheid et al. 2013):

$$K_{\text{N}}=K_{\text{R}}=\omega \exp \left(\frac{z-z_{\text{u}}}{H_{\text{u}}}\right),$$(12)

where ω is the wave frequency; ɀ_u refers to the height of the sponge layer and equals 200 km in this study; H_u represents the damping scale height.

2.1.2 Ionosphere

The governing equations describing the behavior of plasmas in the ionosphere also contain the mass continuity equation, momentum equation, energy equation, and ideal gas law, and are given by the following:

$$\frac{\partial n_{\alpha }}{\partial t}+\left(V_{\alpha }\cdot \nabla \right)n_{\alpha }+n_{\alpha }\nabla \cdot V_{\alpha }=P_{\alpha }-L_{\alpha },$$(13a)

$$\begin{array}{l}{m}_{\alpha }\left[\frac{\partial V_{\alpha }}{\partial t}+\left(V_{\alpha }\cdot \nabla \right)V_{\alpha }\right]=-\frac{\nabla p_{\alpha }}{n_{\alpha }}+q_{\alpha }\left(E+V_{\alpha }\times B\right)\\ +m_{\alpha }{v}_{\alpha n}\left(V_{\text{n}}-V_{\alpha }\right)+m_{\alpha }g,\end{array}$$(13b)

$$\frac{\partial p_{\alpha }}{\partial t}+\left(V_{\alpha }\cdot \nabla \right)p_{\alpha }+{\gamma }_{\alpha }{p}_{\alpha }\left(\nabla \cdot V_{\alpha }\right)+\left({\gamma }_{\alpha }-1\right)\nabla \cdot q_{\alpha }=\left({\gamma }_{\alpha }-1\right)Q_{\alpha n},$$(13c)

and

$$p_{\alpha }=n_{\alpha }{k}_{\text{B}}{T}_{\alpha }\text{, }$$(13d)

respectively, where subscript ‘α’ represents ions or electrons.

Equation (13a) describes the plasma mass continuity. n_α is the number density of plasma; the two terms on the right-hand side, P_α and L_α, are the plasma production rate and loss rate per unit volume due to the photochemical reactions in the ionosphere. In Eq. (13b), m_α and p_α are the plasma mass and pressure, respectively; q_α denotes the plasma electric charge; q_e = −e for electron and q_i for ion; E = −∇ϕ and B are the electric field and magnetic field, respectively. ϕ denotes the electric potential. Equation (13c) describes the plasma energy conservation. γ_α = 5/3 for electron and ion; q_α is the plasma thermal flow vector. q_i = −λ_i∇T_i and q_e = −λ_e∇T_e − β_ej with λ_α being the plasma thermal conductivity coefficient and β_e being the thermoelectric coefficient; T_α is the plasma temperature; j denotes the current density (Rees 1989). Q_αn in Eq. (13c) represents the plasma-neutral collisional heating and is given by (e.g., Huba et al. 2000):

$$\begin{array}{l}{Q}_{\alpha n}=n_{\alpha }{\tilde{v}}_{\alpha n}{k}_{\text{B}}\left(T_{\text{n}}-T_{\alpha }\right)\\ +\frac{1}{2}{m}_{\text{n}}{n}_{\alpha }{v}_{\alpha n}\left(V_{\text{n}}-V_{\alpha }\right)\cdot \left(V_{\text{n}}-V_{\alpha }\right)\end{array}$$(14)

where ${\tilde{v}}_{\alpha n}=3m_{\alpha }{m}_{\text{n}}{v}_{\alpha n}/{\left(m_{\alpha }+m_{\text{n}}\right)}^2$ with m being the mass and ν_αn being the plasma-neutral collisional frequency.

The ionosphere is strongly influenced by the chemical composition of the relevant neutral atmosphere. CO₂ is the most abundant species in the Martian neutral atmosphere. In the Martian dayside ionosphere, plasma is predominantly produced by photoionization by absorption of solar extreme ultraviolet (EUV) radiation. The chemical reactions of major compositions in the Martian ionosphere are listed in Table A.1. The major charged compositions in the ionosphere include, ${\text{O}}^{+},{\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$. We assume that the plasma production rate equals the loss rate due to the chemical equilibrium in the lower-altitude ionosphere (below ~200 km). Based on this assumption, the background number density of ions and electrons are calculated and plotted in Fig. 1b. The smooth temperature profiles of ions and electrons are adopted from Wu et al. (2021) and marked by T_i and T_e in Fig. 1c, respectively. We assume that all the ion species have the same temperature: $T_{\text{i}}=T_{{\text{O}}^{+}}=T_{{\text{O}}_2^{+}}=T_{{\text{CO}}_2^{+}}$. The ion and electron pressures are shown in Fig. 1d. We note that the major ions in the ionosphere is ${\text{O}}_2^{\text{+}}$, we only plot the pressure of ${\text{O}}_2^{\text{+}}$ in Fig. 1d.

The acceleration terms on the left-hand side of Eq. (13b) can be ignored since the plasma response timescale is comparable to the period of AGWs and much larger than plasma gyration period and plasma-neutral collision relaxation timescale. In this case, the plasma momentum balances yield the following relation (Rishbeth & Garriott 1969):

$$V_{\alpha }=V_{\text{n}}+{\mu }_{\alpha }\cdot F_{\alpha },$$(15)

where µ_α refers to the plasma mobility tensor; F_a is given by:

$$F_{\alpha }=-\frac{\nabla p_{\alpha }}{n_{\alpha }}+q_{\alpha }\left(E+V_{\text{n}}\times B\right)+m_{\alpha }g$$(16)

We note that the electric charge of ions is qi = +e since the ions considered in this study include O⁺, ${\text{O}}^{+},{\text{O}}_2^{+},{\text{CO}}_2^{+}$.

The mobility tensor for ions and electrons, µ_i and µ_e, in a geomagnetic coordinate system (O−x′y′ɀ′) can be written as (Rishbeth & Garriott 1969)

$${\mu }_{\text{i}}=\left(\begin{array}{ccc}{\mu }_{\text{pi}}& {\mu }_{\text{Hi}}& 0\\ -{\mu }_{\text{Hi}}& {\mu }_{\text{pi}}& 0\\ 0& 0& {\mu }_{0i}\end{array}\right)$$(17a)

and

$${\mu }_{\text{e}}=\left(\begin{array}{ccc}{\mu }_{\text{pe}}& -{\mu }_{\text{He}}& 0\\ {\mu }_{\text{He}}& {\mu }_{\text{pe}}& 0\\ 0& 0& {\mu }_{0\text{e}}\end{array}\right),$$(17b)

respectively.

The geomagnetic coordinate system O–x′y′ɀ′, is obtained by clockwise rotation of the geographic coordinates by an angle of π/2 + I around the y-axis and an angle of π/2 − D around the ɀ-axis. as illustrated in Fig. 2. Here I denotes the geomagnetic dip angle and D is the geomagnetic declination angle, defined by the incidence of the magnetic field with respect to the geographic north. In coordinates O- x′y′ɀ′, the positive ɀ′-axis is a magnetic field-aligned axis, the positive x′-axis is toward Pedersen direction, and y′ -axis towards the Hall direction and completes the right-handed system,

The elements in the mobility tensor mentioned in Eq. (17a) and (17b) are given by:

$${\mu }_{0\alpha }=\frac{1}{m_{\alpha }{v}_{\alpha n}}=\frac{{\eta }_{\alpha }}{eB},$$(18a)

$${\mu }_{\text{p}\alpha }=\frac{{\mu }_{0\alpha }}{1+{\eta }_{\alpha }^2}=\frac{1}{eB}\frac{{\eta }_{\alpha }}{1+{\eta }_{\alpha }^2},$$(18b)

$${\mu }_{\text{H}\alpha }=\frac{{\mu }_{0\alpha }{\eta }_{\alpha }}{1+{\eta }_{\alpha }^2}=\frac{1}{eB}\frac{{\eta }_{\alpha }^2}{1+{\eta }_{\alpha }^2},$$(18c)

where µ_0α is the magnetic field-aligned mobility per unit charge; µ_pα and µ_Hα are the mobility along the Pedersen and Hall directions, respectively; η_α = Ω_α/ν_αn is the ratio of plasma gyro-frequency to plasma-neutral collisional frequency with Ω_α = eB/m_α being the plasma gyro-frequency.

The plasma-neutral collisional frequencies are adopted from Schunk & Nagy (2009). The profiles of electron- and ion-neutral collisional frequency and gyro-frequency are displayed in Fig. 4a. Figure 4b displays the vertical profiles of ion and electron mobilities. We note that we only consider ${\text{O}}_2^{\text{+}}$ at here as it is the dominantion in the ionosphere. In addition, we assume that the magnetic field magnitude is 20 nT since the magnetic fields are dominated by the draped interplanetary magnetic field outside of the crustal magnetic field region.

The definition of the current density is:

$$j=e{n}_{\text{i}}{V}_{\text{i}}-e{n}_{\text{e}}{V}_{\text{e}}$$(19)

Combining with Eqs. (15) and (19), the current density can be rewritten as:

$$\begin{array}{l}j=e\left(n_{\text{i}}-n_{\text{e}}\right)V_{\text{n}}+{\sigma }_{\text{i}}\cdot \left(-\frac{\nabla p_{\text{i}}}{n_{\text{i}}e}+\frac{m_{\text{i}}g}{e}\right)\\ -{\sigma }_{\text{e}}\cdot \left(-\frac{\nabla p_{\text{e}}}{n_{\text{e}}e}+\frac{m_{\text{e}}g}{e}\right)+\sigma \cdot \left(E+V_{\text{n}}\times B\right),\end{array}$$(20)

where σ_i = n_ie²µ_i and σ_e = n_ee²µ_e are the ion and electron conductivity tensors, respectively; σ = σ_i + σ_e is the total conductivity tensor and can be expressed as the following in the magnetic coordinates O-x′y′z′:

$$\sigma =\left(\begin{array}{ccc}{\sigma }_{\text{p}}& -{\sigma }_{\text{H}}& 0\\ {\sigma }_{\text{H}}& {\sigma }_{\text{p}}& 0\\ 0& 0& {\sigma }_0\end{array}\right),$$(21)

where σ₀, σ_p, and σ_H are the field-aligned, Pedersen, and Hall conductivities, and can be expressed as:

$${\sigma }_0={\sigma }_{0i}+{\sigma }_{0\text{e}}=n_{\text{i}}{\text{e}}^2{\mu }_{0i}+n_{\text{e}}{\text{e}}^2{\mu }_{0\text{e}}$$(22a)

$${\sigma }_{\text{p}}={\sigma }_{\text{pi}}+{\sigma }_{\text{pe}}=\frac{{\sigma }_{0i}}{1+{\eta }_{\text{i}}^2}+\frac{{\sigma }_{0\text{e}}}{1+{\eta }_{\text{e}}^2}$$(22b)

and

$${\sigma }_{\text{H}}={\sigma }_{\text{Hi}}+{\sigma }_{\text{He}}=-\frac{{\sigma }_{0i}{\eta }_{\text{i}}}{1+{\eta }_{\text{i}}^2}+\frac{{\sigma }_{0\text{e}}{\eta }_{\text{e}}}{1+{\eta }_{\text{e}}^2}$$(22c)

respectively.

As displayed in Fig. 4c, the conductivity in the ionosphere is dominated by the field-aligned conductivity, σ₀.

The current density and electric field satisfy:

$$\frac{\partial e\left(n_{\text{i}}-n_{\text{e}}\right)}{\partial t}+\nabla \cdot j=0,\nabla \cdot E=\frac{e\left(n_{\text{i}}-n_{\text{e}}\right)}{{ϵ}_0}$$(23)

where ϵ₀ is the vacuum permittivity.

The ionosphere is quasi-neutral in the absence of AGWs. Generally, the divergence of the undisturbed current density and electric field are assumed to be zero.

To examine the plasma motion in the ionosphere, we first investigate the plasma velocity. Combining Eqs. (15)–(18), yields:

$$V_{\text{i}}=V_{\text{n}}^i+{\mu }_{0i}{\tilde{F}}_{\text{i}}^{//}+{\mu }_{\text{pi}}{\tilde{F}}_{\text{i}}^{\perp }+{\mu }_{\text{Hi}}{\tilde{F}}_{\text{i}}^{\perp }\times e_{\text{B}},$$(24a)

$$V_{\text{e}}=V_{\text{n}}^e+{\mu }_{0\text{e}}{\tilde{F}}_{\text{i}}^{//}+{\mu }_{\text{pe}}{\tilde{F}}_{\text{e}}^{\perp }-{\mu }_{\text{He}}{\tilde{F}}_{\text{e}}^{\perp }\times e_{\text{B}},$$(24b)

where superscript “//” and “⊥” refer to the directions of parallel or perpendicular to the magnetic field, respectively; e_B = B/B is the field-aligned unit vector; and $V_{\text{n}}^{\alpha }$ and ${\tilde{F}}_{\alpha }$ are given by:

$$\begin{array}{l}{V}_{\text{n}}^{\alpha }=V_{\text{n}}+q_{\alpha }{\mu }_{\alpha }\cdot \left(V_{\text{n}}\times B\right)\\ =\frac{1}{1+{\eta }_{\alpha }^2}\left[V_{\text{n}}+{\eta }_{\alpha }^2\left(V_{\text{n}}\cdot e_{\text{B}}\right)e_{\text{B}}+{\eta }_{\alpha }{V}_{\text{n}}\times e_{\text{B}}\right],\end{array}$$(25)

and

$${\tilde{F}}_{\alpha }=-\frac{\nabla p_{\alpha }}{n_{\alpha }}+q_{\alpha }E+m_{\alpha }g,$$(26)

respectively, where $V_{\text{n}}^{\alpha }$ describes the plasma drift motion induced by the neutral wind (e.g., MacLeod 1966; Matcheva et al. 2001), and ${\tilde{F}}_{\alpha }$ contains the effects of the plasma pressure gradient, electric field, and gravity.

As shown in Fig. 4a, the electron-neutral collisional frequency is less than its gyro-frequency: ν_en ≪ Ω_e, but ν_in ≫ Ω_i for ions in a region, called ‘dynamo region’, whose height ranges from ~135 km to ~200 km. This result indicates that the electron motion is dominated by the gyro-motion along the magnetic field above ~135 km. The ions are dragged by the neutral species and are separated from the electrons, forming a strong current in the dynamo region. In this scenario, the ratio of ion gyro-frequency to ion-neutral collisional frequency, ηi, is much less than 1, meaning that µ_Hi ≪ µ_pi ≈ µ_0i; = η_i/eB. Similarly, η_e ≫ 1, leading to the electron mobilities are dominated by the field-aligned component, µ_0e. In this scenario, the electron drift velocity driven by the neutral wind is $V_{\text{n}}^e\approx \left(V_{\text{n}}\cdot e_{\text{B}}\right)e_{\text{B}}$ above ~130 km, which has been used to investigate the gravity wave-induced electron density variations by earlier works (e.g., Hooke 1968; Matcheva et al. 2001; Wang et al. 2023). In addition, the ion mobility is less than the electron mobility, as shown in Fig. 4b, suggesting the ion pressure gradient and electric field acting on ion species are smaller than electrons and can be neglected. Therefore, the ion motions are dominated by the ion-neutral collisions below ~200 km and the neutral wind induced-drift velocity of ions in the dynamo region can be assumed to be equal to the perturbed neutral velocity:

$$V_{\text{i}}\approx V_{\text{n}}.$$(27)

The electron mobilities are dominated by the field-aligned component, indicating that the electron velocity mentioned in Eq. (24b) can be assumed that:

$$V_{\text{e}}=V_{\text{n}}^e+{\mu }_{0\text{e}}{\tilde{F}}_{\text{e}}^{//}\approx \left(V_{\text{n}}\cdot e_{\text{B}}\right)e_{\text{B}}-{\mu }_{0\text{e}}\frac{{\nabla }^{//}{p}_{\text{e}}}{n_{\text{e}}}-{\mu }_{0\text{e}}e{E}^{//},$$(28)

where gravity is not taken into consideration since the electron mass is negligible.

This result indicates that the electron velocity in the dynamo region is determined by the projection of the velocities driven by the neutral wind, pressure gradient, and electric field on the magnetic field lines. The ions are mainly dragged by the neutral species via collisions. The current density in this scenario becomes:

$$j={\sigma }_{0\text{e}}\frac{{\nabla }^{//}{p}_{\text{e}}}{n_{\text{e}}e}+\sigma \cdot \left(E+V_{\text{n}}\times B\right),$$(29)

where the quasi-neutral condition in the absence of AGWs is assumed n_i ≈ n_e.

Additionally, under the assumption of plasma momentum balances, the plasma-neutral collisional force, mentioned in Eq. (6), can be given by:

$$f_{\text{nc }}=\nabla p_{\text{i}}+\nabla p_{\text{e}}-\left(n_{\text{i}}-n_{\text{e}}\right)E-j\times B-m_{\text{i}}{n}_{\text{i}}ɡ-m_{\text{e}}{n}_{\text{e}}ɡ.$$(30)

In the ionosphere-thermospheric region, the neutral pressure is much larger than the plasma pressure. The plasma gravity is negligible. the ionosphere is quasi-neutral and the background electric field is ignored in the absence of AGWs. Hence, the collisional force acting on the neutral constituents is dominated by the Lorentz force and can be expressed as:

$$\begin{array}{c}{f}_{\text{nc}}=-j\times B=-\left[\sigma \cdot \left(V_{\text{n}}\times B\right)\right]\times B\\ =-{\sigma }_{\text{p}}{B}^2\left[\left(V_{\text{n}}\cdot e_{\text{B}}\right)e_{\text{B}}-V_{\text{n}}\right]-{\sigma }_{\text{H}}{B}^2{V}_{\text{n}}\times e_{\text{B}}.\end{array}$$(31)

Equations (1a)–(1d) and Eqs. (13a)–(13d) form a system of coupling equations that describe the wave propagation in the ionosphere-thermospheric region. We investigated the response of AGWs to the plasma behaviors in the dynamo region because the densities of the major charged constituents have the peak values located in the dynamo region (see Fig. 1b). We note that, in this situation, the plasma momentum Eq. (13b) have been replaced by Eqs. (27) and (28), which describe the plasma motions in the dynamo region. We will find the linear wave solution of the system of coupling equations in the next section.

2.2 System of linear equation

The coupling system can be linearized using a small-perturbation method, in which the first order perturbed values are retained. Hence, the variable to be solved, Φ, can be expressed by: Φ = Φ₀ + Φ₁, where Φ₀ represents the variable in a steady-state atmosphere. Φ1 is an unknown variable and can be expressed as a plane wave-like solution in the form of Φ₁(x, y, z, t) = δΦ₁(ɀ)exp[i(k_xx + k_yy − ωt)] in the local geographic coordinates O–xyɀ, where k_x, k_y and ω are the horizontal wavenumbers along the x- and y-axis and wave frequency, respectively. Under this assumption, the differential operators ∇ and ∂/∂t are equivalent to ik and −iω, respectively, where ik = (ik_x, ik_y, ∂/∂z) with k_x and k_y being the horizontal wavenumber along x- and y-axis, respectively.

The linear forms of Eqs. (1a)–(1d) are given by:

$$-i\omega \frac{{\rho }_{\text{n}1}}{{\rho }_{\text{n}0}}-\frac{V_{\text{n}1}\cdot e_z}{\text{H}}+ik\cdot V_{\text{n}1}=0,$$(32a)

$$-i\omega V_{\text{n}1}+i\frac{p_{\text{n}1}}{{\rho }_{\text{n}0}}k-\frac{{\rho }_{\text{n}1}}{{\rho }_{\text{n}0}}ɡ=f_{\text{n}1},$$(32b)

$$-i\omega s_1+\frac{\partial s_0}{\partial z}\left(V_{\text{n}1}\cdot e_z\right)=\frac{Q_{\text{n}1}}{T_{\text{n}0}},$$(32c)

and

$$\frac{p_{\text{n}1}}{p_{\text{n}0}}=\frac{{\rho }_{\text{n}1}}{{\rho }_{\text{n}0}}+\frac{T_{\text{n}1}}{T_{\text{n}0}},$$(32d)

respectively, where e_z is the unit vector along the ɀ-axis; H = (−∂ln ρ_n0/∂ɀ)⁻¹ is the neutral density scale height; the subscripts ‘0’ and ‘1’ refer to the values in the steady-state and the wave-perturbed state, respectively. According to the concept of entropy, the relation between wave-driven entropy and temperature is given by:

$$\frac{s_1}{c_{\text{p}}}=\frac{T_{\text{n}1}}{T_{\text{n}0}}-\frac{R}{c_{\text{p}}}\frac{p_{\text{n}1}}{p_{\text{n}0}}.$$(33)

Here, ƒ_n1 = F_n1/ρ_n0 and Q_n1, the wave-driven force per unit mass and heating input per unit mass, are a linear combination of V_n1, Tn1, and p_n1.

The linear system governing the plasma behaviors induced by AGWs in the dynamo region can be expressed as:

$$-i\omega \frac{n_{\alpha 1}}{n_{\alpha 0}}-\frac{V_{\alpha 1}\cdot e_z}{H_{\alpha }}+ik\cdot V_{\alpha 1}=\frac{P_{\alpha 1}-L_{\alpha 1}}{n_{\alpha 0}},$$(34a)

$$V_{\text{i}}=V_{\text{n}1},V_{\text{e}1}=V_{\text{n}1}^{//}-{\mu }_{0\text{e}}\frac{p_{\text{e}1}}{n_{\text{e}0}}i{\xi }_{\text{e}}-{\mu }_{0\text{e}}e{E}_1^{//},$$(34b)

$$\begin{array}{c}-i\omega \frac{p_{\alpha 1}}{p_{\alpha 0}}-\frac{V_{\alpha 1}\cdot e_z}{H_{\alpha }^{\ast }}+{\gamma }_{\alpha }\left(ik\cdot V_{\alpha 1}\right)\\ +\left({\gamma }_{\alpha }-1\right)\frac{ik\cdot q_{\alpha 1}}{p_{\alpha 0}}=\left({\gamma }_{\alpha }-1\right)\frac{Q_{\alpha n1}}{p_{\alpha 0}},\end{array}$$(34c)

and

$$\frac{p_{\alpha 1}}{p_{\alpha 0}}=\frac{{\rho }_{\alpha 1}}{{\rho }_{\alpha 0}}+\frac{T_{\alpha 1}}{T_{\alpha 0}},$$(34d)

where H_α = (−∂ln n_α0/∂ɀ)⁻¹ and $H_{\alpha }^{*}={\left(-\partial \ln p_{\alpha 0}/\partial z\right)}^{-1}$ are the scale heights of plasma density and pressure, respectively; $V_{\text{n}1}^{//}=\left(V_{\text{n}1}\cdot e_{\text{B}}\right)e_{\text{B}};{\xi }_{\text{e}}=\left(k_{\Vert }-i\sin I/H_{\text{e}}^{*}\right)e_{\text{B}}$, where k_|| = k · e_в; $E_1^{//}=\left(E_1\cdot e_{\text{B}}\right)e_{\text{B}}$ are the field-aligned wave velocity, wavenumber vector, and electric field, respectively. q_α1 refers to the wave-induced plasma thermal conduction; Q_αn1 is the wave-driven plasma collisional heating with neutral constituents and is a linear combination of n_α1, p_α1, and p_n1.

In the dynamo region, the wave-induced ion density does not equal the wave-induced electron density, suggesting that the quasi-neutral condition is not satisfied in the presence of AGWs. In this scenario, the divergence of the wave-induced electric field and current density are not equal to zero and satisfied:

$$\begin{array}{c}-i\omega e\left(n_{\text{i}1}-n_{\text{e}1}\right)+ik\cdot j_1=0,\\ \frac{e\left(n_{\text{i}1}-n_{\text{e}1}\right)}{{ϵ}_0}=ik\cdot E_1=(k\cdot k){\varphi }_1,\end{array}$$(35)

where ϵ₀ denotes vacuum permittivity; ϕ₁ is the wave-driven electric potential and satisfies E₁ = −ikϕ₁; and the wave-induced current density in the dynamo region is

$$j_1=i{\sigma }_{0\text{e}}\frac{p_{\text{e}1}}{n_{\text{e}0}e}{\xi }_{\text{e}}+\sigma \cdot \left(-ik{\varphi }_1+V_{\text{n}1}\times B_0\right)$$(36)

Hence, the relation between current density and electric potential can be given by:

$$\begin{array}{c}\left[-i\omega {ϵ}_0(k\cdot k)+(k\cdot \sigma \cdot k)\right]{\varphi }_1\\ -{\sigma }_{0\text{e}}\frac{p_{\text{e}1}}{n_{\text{e}0}e}{k}_{\left|\right|}{\xi }_{\text{e}}+ik\cdot \left[\sigma \cdot \left(V_{\text{n}1}\times B_0\right)\right]=0\end{array}$$(37)

Combining Eqs. (32a)–(34d) and Eq. (37), yields the following coupling linear system between neutral atmosphere and ionosphere:

$$\tilde{A}\frac{{\partial }^{2}X}{\partial z^2}+\tilde{B}\frac{\partial X}{\partial z}+\tilde{C}X=0,$$(38)

where $\tilde{A},\tilde{B}$, and $\tilde{C}$ are 10 × 10 coefficient matrices; X is the vector to be solved and is given by:

$$X={\left(\frac{n_{\text{i}1}}{n_{\text{i}0}},\frac{p_{\text{i}1}}{p_{\text{i}0}},\frac{n_{\text{e}1}}{n_{\text{e}0}},\frac{p_{\text{e}1}}{p_{\text{e}0}},{\varphi }_1,\frac{p_{\text{n}1}}{p_{\text{n}0}},\frac{s_1}{c_{\text{p}}},u_1,v_1,w_1\right)}^T,$$(39)

where superscript “T” denotes the transform of a matrix.

The system of linear Eq. (38) can be expressed as $\widehat{D}X=0$ with D being a differential operator matrix:

$$\widehat{D}=\tilde{A}\frac{{\partial }^2}{\partial z^2}+\tilde{B}\frac{\partial }{\partial z}+\tilde{C}.$$(40)

The 10 × 10 matrix can be partitioned into four blocks:

$$\widehat{D}=\left(\begin{array}{ll}P\hfill & Q\hfill \\ R\hfill & S\hfill \end{array}\right),$$(41)

where the small four matrices P, Q, R, and S are 5 × 5 matrices. Each matrix is a linear combination of two differential operators: ∂²/∂ɀ² and ∂/∂ɀ. The elements of the four matrices are displayed in Appendix B.

X can also be divided into two parts X_ionos and X_neu, where

$$\begin{array}{c}{X}_{\text{ionos }}={\left(\frac{n_{\text{il}}}{n_{\text{i}0}},\frac{p_{\text{i}1}}{p_{\text{i}0}},\frac{n_{\text{e}1}}{n_{\text{e}0}},\frac{p_{\text{e}1}}{p_{\text{e}0}},{\varphi }_1\right)}^T,\\ X_{\text{neu }}={\left(\frac{p_{\text{n}1}}{p_{\text{n}0}},\frac{s_1}{c_{\text{p}}},u_1,v_1,w_1\right)}^T.\end{array}$$(42)

Hence, the linear system (38) can be rewritten as

$$\widehat{D}X=\left(\begin{array}{cc}P& Q\\ R& S\end{array}\right)\left(\begin{array}{c}{X}_{\text{ionos }}\\ X_{\text{neu }}\end{array}\right)=0.$$(43)

Q and R are called the coupling matrices between the neutral atmosphere and ionosphere. As shown in Fig. 5, Q represents the responses of the plasmas to the neutral constituents via collisions and thermal exchange; When waves propagate into the ionosphere, the plasma motions are governed by the wave velocity, contributing to the fluctuations in the plasma density. In addition, the thermal exchange exist between plasma and neutrals via collision, causing the plasma to obtain the energy from waves. Similarly, R refers to the impacts of the plasmas on the neutral species by collisions and thermal exchange. The waves propagated into the ionospheric region can also be impacted by the plasma behaviors and can lose energy due to thermal exchange. The elements of the coupling matrix R consist of the neutral-plasma thermal exchange and photochemical heating terms, whose rates are defined by a characteristic neutral-plasma energetic exchange frequency, ω_nα, and photochemical heating frequency, ω_pc, respectively (see Appendix B). ω_nα and ω_pc are much less than the Brunt-Väisälä frequency, ω_b, Therefore, the coupling matrix R has an insignificant contribution to the AGWs.

Additionally, waves propagating into the upper atmosphere may dissipate due to ion friction (e.g., Medvedev et al. 2017) and Pedersen and Hall conductivities (e.g., Miesen et al. 1989; Miesen 1991). Throughout this dissipation process, the frictional force acting on the neutral species is primarily governed by the Lorentz force, as defined by Eq. (31). The magnitude of Pedersen and Hall conductivities, as well as the orientation of the local magnetic field, determine the magnitude of the frictional force and eventually dampen wave amplitude. Matrix S consists of this dissipation process, with it quantitative description described by Eq. (B.8). The magnetic field of the present Mars is dominated by the horizontal draped interplanetary magnetic field, with a magnitude in the tens of nT. The Pedersen/(Hall) conductivity damping rate, defined by a characteristic frequency, ω_p/(ω_H), is much less than the Brunt-Väisälä frequency, as displayed in Fig. B.1. This implies that the impact of conductivity damping on the waves can be reasonably ignored when the background magnetic fields are dominated by the draped interplanetary magnetic field. However, in the crustal magnetic field region, the magnetic field topology is complex, with the magnetic field strengths ranging from several hundred to several thousand nT (Acuna et al. 1998). In this scenario, waves are more easily dissipated by the conductivity and orientation of the magnetic field. The wave dissipation process in the crustal magnetic field region is complicated and beyond of this work. Further study is required to fully understand it.

This paper mainly focuses on the ionospheric response to AGWs outside of the crustal magnetic field region. The coupling matrix R can be ignored in this study. Under this assumption, the linear system becomes

$$S{X}_{\text{neu }}=0,$$(44a)

$$P{X}_{\text{ionos }}=-Q{X}_{\text{neu }}.$$(44b)

The linear system can be solved by the following steps. First, we calculate the wave-driven neutral velocity, density, temperature, and pressure by solving Eq. (44a). Next, the ionospheric perturbations driven by AGWs, X_ionos, can be solved by Eq. (44b). We provide a detailed description of the neutral wave solution in Sect. 2.3. The model results, including the wave-driven variations in plasma density, velocity, pressure, and temperature, the perturbed electromagnetic field, and the wave-associated thermal variations, are displayed in Sect. 3.

2.3 Neutral wave solution

In this section, we describe the neutral behaviors driven by the AGWs by solving linear Eq. (44a). As mentioned in the above section, the differential operator matrix S is a linear combination of ∂²/∂ɀ² and ∂/∂ɀ. The linear Eq. (44a) is equivalent to the following:

$$S_1\frac{{\partial }^2{X}_{\text{neu }}}{\partial z^2}+S_2\frac{\partial X_{\text{neu }}}{\partial z}+S_3=0,$$(45)

where S₁, S₂, and S₃ are the 5 × 5 coefficient matrices, whose elements are functions of altitude. These three matrices are too tedious and are shown in Appendix B.

The linear second-order differential Eq. (45) for X_neu can be solved iteratively by a numerical method since their coefficients are not constant. A finite-difference method is applied to solve the differential equation. The upper and lower boundary conditions, called Dirichlet boundary condition, should be given. More details of the numerical method are displayed in Appendix D.

2.3.1 Upper boundary conditions

The upper atmosphere can typically be treated as a quasi isothermal structure, as shown in Fig. 1c. In this scenario, the differential Eq. (45) is equivalent to the following (see Appendix C)

$$\frac{{\partial }^2{\widehat{w}}_1}{\partial z^2}+k_z^2{\widehat{w}}_1=0,$$(46)

where ${\widehat{\omega }}_1$ and k_ɀ are given by

$${\widehat{w}}_1=w_1\exp \left(-{{\displaystyle \int }}^{\text{}}\frac{1}{2H}\text{d}z\right),$$(47)

and

$$\begin{array}{c}{k}_z^2=\frac{k_{\text{h}}^2{\omega }_{\text{b}}^2}{\left(\omega +i{\kappa }_{\text{m}}\alpha \right)\left(\omega +i{v}_{\text{m}}\alpha \right)}\\ +\frac{\omega \left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)}{\left(\omega +i{\kappa }_{\text{m}}\alpha \right)c_{\text{s}}^2}\left[\omega +i{v}_{\text{m}}\alpha \left(1+\frac{1}{3}\right)\right]\\ -k_{\text{h}}^2-\frac{1}{4H^2}\end{array}$$(48)

respectively. Here $k_{\text{h}}=\sqrt{k_x^2+k_y^2}$ and k_ɀ refer to the horizontal and vertical wavenumbers; $\alpha =k_{\text{h}}^2-{\left(i{k}_z+1/2H\right)}^2$.

Equation (48), a sixth-order polynomial of k_ɀ, is the dispersion relation of AGWs associated with the dissipative processes including molecular viscosity and thermal diffusivity. It consists of three pairs of wave solutions describing the upward/downward propagation of AGWs, ordinary viscosity waves, and ordinary heat conduction waves (Volland 1969b). The pair of solutions describing AGWs is supposed to satisfy the condition that $k_z^2\approx k_{\text{h}}^2\left({\omega }_{\text{b}}^2/{\omega }^2-1\right)$.

k_ɀ varies slowly with height, suggesting that Eq. (46) must be solved by using the WKB approximate method (e.g., Yeh & Liu 1974; Matcheva & Strobel 1999; Wang et al. 2020, 2023). The WKB wave solution is:

$$\begin{array}{c}{w}_1(x,y,z,t)={\widehat{w}}_1\exp \left({{\displaystyle \int }}^{\text{}}\frac{1}{2H}\text{d}z\right)\\ =\text{Δ}{w}_1(z)\exp [i\psi (x,y,z,t)]\end{array}$$(49)

where ∆w₁ and ψ are the wave amplitude and phase and can be given by:

$$\text{Δ}{w}_1(z)=\delta \omega \left(z_0\right)\sqrt{\frac{k_{zr}\left(z_0\right)}{k_{zr}(z)}}\exp \left[{{{\displaystyle \int }}^{\text{}}}_{z_0}^z\left(\frac{1}{2H}-k_{zi}\right)\text{d}z\right]\text{, }$$(50)

and

$$\psi (x,y,z,t)=k_{x}x+k_{y}y+{{{\displaystyle \int }}^{\text{}}}_{z_0}^z{k}_{zr}\text{d}z-\omega t,$$(51)

respectively. Here, δw(ɀ₀) denotes the corresponding wave amplitude at the reference altitude, ɀ₀; k_ɀr and k_ɀi are the real and imaginary parts of k_ɀ, respectively.

Other wave-induced variations, including u₁, υ₁, ρ_n1, p_n1, T_n1, and s₁ are related to w₁ and can be obtained by the wave polarization relations. More details of the polarization relations are displayed in Appendix C. This WKB wave solution can be treated as the upper boundary condition.

2.3.2 Lower boundary conditions

In the lower atmosphere, the effects of molecular viscosity and thermal conductivity are negligible. The wave dispersion relation, in this scenario, becomes:

$$k_z^2=\frac{{\omega }^2-{\omega }_{\text{a}}^2}{c_{\text{s}}^2}+k_{\text{h}}^2\frac{{\omega }_{\text{b}}^2-{\omega }^2}{{\omega }^2},$$(52)

which usually describes wave propagation in compressible and non-dissipativeatmosphere(Hines 1960). It consists of two pairs of wave solutions. One pair is the acoustic wave branch and the other pair refers to the internal gravity wave branch.

k_ɀ solved by Eq. (52) only has the real part. the imaginary part equals to zero. Hence, the wave solution in the non-dissipative atmosphere is given by

$$w_1(x,y,z,t)=\delta w_1\left(z_0\right)\exp \left[i\psi +{{{\displaystyle \int }}^{\text{}}}_{z_0}^z\left(i{k}_z+\frac{1}{2H}\right)\text{d}z\right],$$(53)

where ψ = k_xx+ k_yy − ωt.

The polarization relations are given by the following:

$$\frac{p_{\text{n}1}}{p_{\text{n}0}}=\frac{\gamma }{i\omega }\frac{{\omega }^2}{{\omega }^2-c_{\text{s}}^2{k}_{\text{h}}^2}\left[i{k}_z+\left(\frac{1}{2}-\frac{1}{\gamma }\right)\frac{1}{H}\right]w_1,$$(54a)

$$\frac{s_1}{c_{\text{p}}}=\frac{i{\omega }_{\text{b}}^2}{g\omega }{w}_1$$(54b)

$$u_1=\frac{i{k}_x{c}_{\text{s}}^2}{c_{\text{s}}^2{k}_{\text{h}}^2-{\omega }^2}\left[i{k}_z+\left(\frac{1}{2}-\frac{1}{\gamma }\right)\frac{1}{H}\right]w_1,$$(54c)

$$v_1=\frac{i{k}_y{c}_{\text{s}}^2}{c_{\text{s}}^2{k}_{\text{h}}^2-{\omega }^2}\left[i{k}_z+\left(\frac{1}{2}-\frac{1}{\gamma }\right)\frac{1}{H}\right]w_1.$$(54d)

This wave solution in the absence of molecular viscosity and thermal conductivity can be treated as the lower boundary condition.

3 Model results

In this section, we examine the impact of AGWs on the Martian ionosphere. The model input paramters, including wave period and horizontal wavelength, are described in Sect. 3.1. The wave-driven plasma motions are displayed in Sect. 3.2. Correspondingly, the wave-driven plasma density variations associated with the photochemistry, electron pressure gradient, and the electric field are given in Sect. 3.3. The wave-driven motions in the ion and electron can contribute to the variations in the current density and electric field, which ultimately leads to the localized magnetic fluctuation. In Sect. 3.4, we investigated the wave-induced variations in the current density and electromagnetic field. Finally, we also examine the thermal varation driven by waves in Sect. 3.5.

3.1 Input parameters

The input parameters of the coupling model include wave period (τ = 2π/ω), wave horizontal wavelength (λ_h = 2π/k_h), and wave amplitude. τ and λ_h at the lower boundary satisfies the non-dissipative wave dispersion relation (52) in the Martian atmosphere. Earlier spacecraft observations suggest that AGWs at Mars have a wide spatial range of approximately ten to 300 km with an amplitude of tens of percent (e.g., Creasey et al. 2006; Withers 2006; Terada et al. 2017; Siddle et al. 2019). Flow-over topography, weather systems, and convection near the Martian surface are the potential generation mechanisms of AGWs (e.g., Parish et al. 2009; Heavens et al. 2020; Yiǧit et al. 2021b). Topographically generated waves can have a variety of frequencies but have near-zero phase velocity. AGWs resulting from the instability of weather systems typically have lower frequencies with diverse spectrum of horizontal phase speeds. Convectively generated waves can possess either lower or higher frequencies and have a variety of wave phase speeds (e.g., Fritts & Alexander 2003; Heavens et al. 2020). Fig. 6 displays λ_h − c_p parameter space of AGWs, where c_p = ω/k_h refers to the wave horizontal phase speed. The relationship between λ_h and c_p can be obtained from Eq. (52). To represent the realistic behavior of waves on Mars, the horizontal wavelengths are typically constrained within a range of 10 to 300 km. The slope of the λ_h − c_p parameter space represents the wave frequency. For acoustic wave branch, the wave frequency is greater than the acoustic cutoff frequency, ω > ω_a; the horizontal phase speeds greater than the local sound speed, c_p > c_s ≈ 170 m s⁻¹, as mentioned in Sect. 2.1.1. Equivalently, the periods of acoustic waves are generally less than τ_a = 2π/ω_a ≈ 8 min. For gravity wave branch, the horizontal phase speeds are less than c_sω_b/ω_a ≈ 127 m s⁻¹; the wave periods are greater than 2π/ω_b ≈ 10 min. The two solid curves in Fig. 6 are the bounding curves such that the acoustic wave branch and the gravity wave branch are bounded by the two curves along which $k_z^l=0$ km. As λ_h → ∞, the bounding curve approach asymptotically to ω = ω_a for the acoustic wave branch. For the gravity wave branch, the bounding curve approaches to ω = ω_b, as λ_h → 0.

It is worthly note that the vertical components of the wave phase propagation and energy propagation have the same signs for the acoustic wave branch and the opposite signs for the gravity wave branch (e.g., Yeh & Liu 1972). We selected the wave modes whose energy propagated from lower atmosphere into upper levels in this study. In this scenario, the vertical wavenum-ber of the upward propagating waves has positive and negative value for the acoustic and gravity wave branches, respectively.

We choose three wave sets for each acoustic and gravity wave branch to explore the impact of AGWs on the ionosphere. The vertical wavenumbers of the selected acoustic set and gravity set at lower boundary are equal to 0.11, 0.32, and 0.99 km⁻¹ and −0.07, −0.19, and −0.5 km⁻¹, respectively. The selected wave sets corresponding to different vertical wavenumbers are plotted by the dash-dotted curves in Fig. 6. The frequency and phase speed of the selected waves vary widely, ranging from lower to higher frequencies and approximately 10 to 1000 m s⁻¹, respectively. Consequently, these selected wave sets have the capability to represent the realistic behaviors of waves generated near the Martian surface.

In an ideal non-dissipative atmosphere, wave energy remains conserved during propagation. As a result, wave amplitude grow exponentially with increasing height due to the decreasing background density (e.g., Yiğit 2023). However, waves tend to be dissipated by the molecular viscosity and thermal conduction in a realistic atmosphere. As displayed in Fig. 3c, the molecular kinetic viscosity, v_m, and thermal diffusivity, κ_m, increase exponentially with height, resulting in wave attenuation and energy release in the upper atmosphere. Therefore, the wave amplitude begins to decay at a certain height, which is called “damping altitude”. The amplitude profiles of the selected acoustic sets with $k_z^l=0.11$ km⁻¹ and gravity sets with $k_z^l=-0.07$ km⁻¹ are displayed in Fig. 7. Waves have higher damping altitudes when they have faster phase speeds (equivalently, longer wavelengths or lower periods). To explore the effect of AGWs on the ionosphere-thermospheric region, we should choose a wave mode with a longer horizontal wavelength or faster phase speed to ensure that the wave damping altitude is located in the dynamo region (~ 135–200 km, see Fig. 4). We selected one wave mode each in the acoustic set and gravity set to examine the ionospheric response to AGWs. The period, horizontal wavelength, phase speed, and damping altitude of the selected acoustic wave mode are ~4.2 min, 80 km, 321 m s⁻¹, and 160 km, respectively. Compared with the selected acoustic mode, the selected gravity mode has a higher period and a longer wavelength, which are ~39.4 min and 210 km, respectively. Correspondingly, its horizontal phase speed and damping altitude are 89 m s⁻¹ and 153 km. To simplify the description of each wave mode, the selected acoustic mode and gravity mode are named by “AW” and “GW", and marked by “★” and “×” in Fig. 6, respectively.

The neutral density (ρ_n1), pressure (p_n1), temperature (T_n1) driven by the AW and GW are calculated by the neutral wave model mentioned in Sect. 2.3 and their ratio amplitudes are displayed in Figs. 8a–c, and marked by |ρ_n1/ρ_n0|, |p_n1/p_n0|, and |T_n1/T_n0|, respectively. We note that the maximum amplitude of wave-associated neutral density ratio is fixed at 10%, which coincides with the spacecraft measurements (Terada et al. 2017; Siddle et al. 2019). We assume that wave propagates in the x-ɀ plane. Correspondingly, the wave azimuth angle, Φ, which defined by the incidence of the wavenumber vector with respect to geographic north (see Fig. 2), is assumed to be 90°; The x and y component of wavenumber vector are equal to k_h and zero (k_x = k_h sin Φ = k_h, k_y = k_h cos φ = 0), respectively; The y component of wave-driven neutral wind velocity, υ₁, equals to zero. The amplitudes of the perturbed neutral wind speed along x and ɀ-axis, u₁ and w₁ are plotted in Figs. 8d and e, respectively. The vertical wavelengths, λ_z = 2π/k_ɀ, of AW and GW are displayed in Fig. 8f.

We explore the electron and ion densities response to AW and GW in the following section. But the plasma motion behaviors driven by waves in the dynamo region should be first investigated.

3.2 Wave-driven plasma velocity

As displayed in Sect. 2.1.2, in the dynamo region, whose height ranges from ~ 135 km to 200 km, ions are primarily influenced by neutral particles through drag or neutral-ionic collision, whereas electrons are mainly subject to gyro-motion along the magnetic field. The wave-associated motion of ions and electrons in the dynamo region are depicted in Fig. 9. Ions, which are dragged by the neutral constituents, oscillate along the wave phase lines. The electron velocity is mainly dominated by the field-aligned component. As waves propagate into the dynamo region, electrons accumulate in certain regions and disperse outward into adjacent regions, forming plasma layer structures. In this scenario, the wave-driven plasma velocity can be described by Eq. (34b). It has been suggested that the wave-induced ion velocity in the dynamo region is assumed to be equal to the perturbed neutral velocity: V_i1 = V_n1, whereas the wave-driven electron velocity can be expressed as a combination of three terms: $V_{\text{el}}=V_{\text{n}1}^{\parallel }+V_{\nabla p_{\text{e}1}}^{\parallel }+V_{\varphi 1}^{\parallel }$. The first term $V_{\text{n}1}^{\parallel }=\left(V_{\text{n}1}\cdot e_{\text{B}}\right)e_{\text{B}}$, represents the projection of the perturbed neutral velocity onto the magnetic field lines. The second term denotes the field-aligned electron velocity driven by the perturbed electron pressure gradient, $V_{\nabla p_{\text{e}1}}^{\parallel }=-i{\mu }_{0\text{e}}{p}_{\text{e}1}{\xi }_{\text{e}}/n_{\text{e}0}=i{\mu }_{0\text{e}}e{\varphi }_{\text{am}1}{\xi }_{\text{e}}$, and the third term represents the velocity driven by the perturbed electric potential, $V_{{\varphi }_1}^{\parallel }=i{\mu }_{0\text{e}}e{\varphi }_1{k}^{\parallel }$. We note that the perturbed electron pressure, p_e1, perturbed ambipolar electric potential, ϕ_am1 = −p_e1/n_e0e, and perturbed total electric potential, ϕ₁, are calculated by the coupling linear system mentioned in Sect. 2.2. The dip angle of the magnetic field is assumed to be equal to zero degree since the magnetic field dominated by the draped interplanetary magnetic field, whose orientation is quasi-parallel to the Martian surface.

Earlier works on other planets (Earth, Jupiter, and Saturn) only considered the contribution of the perturbed neutral velocities to plasma motion in the wave-ionospheric coupling process (e.g., MacLeod 1966; Hooke 1968; Matcheva et al. 2001; Matcheva & Barrow 2012). The effects of electron pressure and electric field on the electron motion are first included in the Martian ionosphere. As displayed in Figs. 10a and c, $V_{\nabla p_{e1}}^{\Vert }\text{and }{V}_{\varphi 1}^{\Vert }$ induced by AW and GW modes can reach several hundred m s⁻¹, indicating that the electron pressure and electric field are significant and should not be neglected when waves interact with the Martian ionosphere. $V_{\nabla p_{e1}}^{\Vert }\text{and }{V}_{\varphi 1}^{\Vert }$ are reversible velocities that should have perfect cancelation in the ionosphere without neutral-drag effects. However, $V_{\nabla p_{e1}}^{\Vert }\text{and }{V}_{\varphi 1}^{\Vert }$ are not completely canceled out at all altitudes. The reason is that $V_{\nabla p_{e1}}^{\Vert }$ is only determined by the ambipolar electric field while $V_{\varphi 1}^{\Vert }$ consists of the effects of the ambipolar and neutral wind-driven polarization electric field (see Eq. (37)). Therefore, the wave-driven electron velocity (V_e1) is not exactly equal to the field-aligned wave velocity $V_{\text{n}1}^{\Vert }$ for AW and GW modes, as displayed in Figs. 10b and d, respectively. This result indicates that the effect of wave-driven polarization electric field on the electron motion is significant and should be taken into consideration when the wave propagating into the Martian ionosphere.

3.3 Wave-driven plasma density variations

We examine the wave-induced plasma density variations in this section. The perturbed plasma density, n_α1, can be obtained by solving the linearized plasma continuity Eq. (34a). To investigate the effects of photochemistry, electron pressure, and electric field on the wave-induced plasma density variations, we considered three scenarios, which are summarized in Table 1. We note that in all three cases, the perturbed ion velocity is assumed to be equal to Vn1 . The major ion constituents in the Martian ionosphere, including O⁺, ${\text{O}}_2^{+}{\text{, and CO}}_2^{+}$, are investigated in the wave-ionospheric coupling process.

Case 1. The effects of photochemistry, electron pressure gradient, and electric field are not considered. In this scenario, the wave-associated electron velocity is simply equal to the perturbed neutral velocity projected on the magnetic field, $V_{\text{e1}}=V_{\text{n}1}^{\Vert }$. The analytical soultion of Eq. (34a) for perturbed electron density and perturbed ion density (including O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$) are given by:

$$\frac{n_{\text{e}1}}{n_{\text{e}0}}=\frac{V_{\text{n}1}^{\Vert }}{i\omega }\left[i\left(k\cdot e_{\text{B}}\right)-\frac{e_{\text{B}}\cdot e_z}{H_{\text{e}}}\right],$$(55a)

and

$$\frac{n_{\text{i}1}}{n_{\text{i}0}}=\frac{1}{i\omega }\left[ik\cdot V_{\text{i}1}-\frac{V_{\text{i}1}\cdot e_z}{H_{\text{i}}}\right]\text{,}$$(55b)

respectively, where H_α = (−∂ln n_α0/∂ɀ)⁻¹ refers to the plasma density scale height. Eq. (55a) describes the small-amplitude perturbed electron density without photochemistry, electron pressure gradient, and electric field effects and has been used in earlier work in the Martian ionosphere (Wang et al. 2023) and even in the other planetary ionosphere (e.g., Hooke 1968; Matcheva et al. 2001; Barrow & Matcheva 2011; Matcheva & Barrow 2012).

Case 2. In this situation, the effect of photochemistry on the wave-ionospheric coupling process is taken into consideration but the contributions of the electron pressure gradient and electric field are ignored. The electron loss channels in the photochemical process include recombination reactions, which are described in Table A.1. The main dissociative recombination reaction in the Martian lower ionosphere, $R_1:{\text{O}}_2^{+}+{\text{e}}^{-}\to \text{O}+\text{O}$, is only considered in this case. Under this assumption, the electron loss rate in the absence of AGWs is $L_{\text{e}0}\approx {\alpha }_1{n}_{\text{e}0}^2$, where α₁ refers to the dissociative recombination coefficient of R₁. Therefore, the wave-associated electron loss rate is given by L_e1 = 2α₁n_e0n_e1. The electron production channels include photodissociation reactions, suggesting that the electron production rate is determined by the density of neutral constituents. We simply assume that the ratio perturbed electron production rate equals to the ratio perturbed neutral density, P_e1/P_e0 = ρ_n1/ρ_n0. In this scenario, the corresponding perturbed electron density is given by the following formula:

$$\frac{n_{\text{e}1}}{n_{\text{e}0}}=\frac{1}{i\omega -2{\alpha }_1{n}_{\text{e}0}}\left[\left(ik\cdot e_{\text{B}}-\frac{e_{\text{B}}\cdot e_z}{H_{\text{e}}}\right)V_{\text{n}1}^{\Vert }-\frac{{\rho }_{\text{n}1}}{{\rho }_{\text{n}0}}\frac{P_{\text{e}0}}{n_{\text{e}0}}\right].$$(56)

We note that the wave-driven electron velocity in this case is also assumed to be equal to the field-aligned neutral velocity induced by waves.

Equation (55b), describing the wave-perturbed ion density, should be revised in this scenario, i.e., the photochemical effects need to added. Therefore, the corresponding solution for perturbed density of O⁺, ${\text{O}}_2^{+}\text{\hspace{0.17em}and\hspace{0.17em}}{\text{CO}}_2^{+}$, in this scenario is given by:

$$\frac{n_{\text{i}1}}{n_{\text{i}0}}=\frac{1}{i\omega }\left[ik\cdot V_{\text{i}1}-\frac{V_{\text{i}1}\cdot e_z}{H_{\text{i}}}-\frac{P_{\text{i}1}-L_{\text{i}1}}{n_{\text{i}0}}\right],$$(57)

where P_i1 and L_i1 refer to the perturbed ion production rate and loss rate, respectively. The loss rates of ${\text{O}}_2^{+}\text{\hspace{0.17em}and\hspace{0.17em}}{\text{CO}}_2^{+}$are determined by the electron density. The production rate and loss rate of these three major ion compositions in the absence of AGWs are summarized in Table A.2. The perturbed densities of O⁺, ${\text{O}}_2^{+}\text{\hspace{0.17em}and\hspace{0.17em}}{\text{CO}}_2^{+}$ associated with photochemistry can be obtained by solving Eq. (A.2).

Case 3. The effects of photochemical reactions, electron pressure gradient, and electric field are all taken into consideration. The wave-driven electron velocity, in this comprehensive scenario, is given by $V_{\text{e}1}=V_{\text{n}1}^{\parallel }+V_{\nabla p_{\text{el}}}^{\parallel }+V_{\varphi 1}^{\parallel }$. The perturbed electron density is described by the solution of the coupling system mentioned in Sect. 2.2, while the perturbed ion density in the presence of photochemistry is also given by the solution of the linear system Eq. (A.3).

Figure 11 displays the wave-associated electron density ratio amplitudes driven by AW and GW mode. The perturbed density for case 1, 2, and 3 are marked by #1, #2, and #3, respectively. When the photochemistry effect is neglected (case 1), the perturbed electron density ratio attains a maximum of ~5% and ~37% induced by AW and GW, respectively. We noted that the maximum neutral density ratio amplitude (|ρ_n1 /ρ_n0|) is fixed at 10% for both AW and GW mode. Compared with AW mode, the period of GW is more larger. Therefore, there are more cumulative effects on the electron density when the larger-period waves propagate into the ionosphere, causing a higher perturbed density ratio. However, the wave-driven electron density becomes weaker in the presence of photochemical reaction (case 2 and 3), and its ratio is less than 4% for AW and GW mode, indicating that the electron density can be reduced by the chemical reaction. Moreover, the impact of photochemistry on the electron density perturbation is more pronounced for GW mode compared to AW mode. The reason is that longer-periodic waves accumulate more photochemical effects. Additionally, when waves propagate into the ionosphere, electrons react with the ambient ions in the presence of photochemistry, resulting in perturbations in the electron loss or production rate. These complex processes ultimately lead to a decrease in electron density fluctuations. These results suggest that the photochemistry effect plays an important role when the electron responses to AGWs in the wave-ionospheric coupling process and cannot be neglected. In addition, the effects of electron pressure gradient and electric field on the perturbed electron velocity are nearly canceled out, contributing to the perturbed electron density for case 2 is close to those for case 3.

The perturbed ion density driven by AW and GW mode are displayed in Figs. 12 and 13, respectively. Similarly, the wave-associated density variation is smaller and its ratio is less than 10% when the photochemical effect is taken into consideration. The ion density variation also decreases due to the fluctuations of the ion loss and production late in the presence of photochemistry. In the absence of photochemistry, the perturbed ion density ratio is close to neutral density variation for the AW mode while reaches tens of percent for the GW mode. Because the ions have more cumulative effects when they interact with large periodic waves. In addition, O⁺ has the largest perturbed density ratio while ${\text{O}}_2^{+}$ is smallest for case 1. The reason is that the wave-driven variation in ion density is proportional to the plasma density scale height (H_i), as shown by Eq. (55b). The ion densities in the absence of AGWs increase with altitude below ~140 km (see Fig. 1b), suggesting that the ion density scale height has negative value. The average density scale heights below ~140 km of O⁺, ${\text{O}}_2^{+}\text{ and\hspace{0.17em}}{\text{CO}}_2^{+}$ are −3.2 km, −13 km, and −3.8 km, respectively. Therefore, O⁺ has the largest density variation compared with ${\text{O}}_2^{+}\text{ and\hspace{0.17em}}{\text{CO}}_2^{+}$ in the absence of photochemistry.

The spatial distributions of the total perturbed ion and electron density (n_α0 + n_α1) for case 3 driven by AW and GW modes are presented in Fig. 14. We note that we only displayed the wave-associated density of the major ion constituent, ${\text{O}}_2^{+}$. The spatial density distributions of O⁺ and ${\text{CO}}_2^{+}$ have the similar structures and have not presented. The electronic and ionic contour lines driven by waves are consistent with the ‘Traveling Ionospheric Disturbances’ (TID) events observed on Earth. The ‘oval eye’ distributions located at levels near ~130 −150 km and ripples located in other regions are the unique structure of wave-driven TIDs (e.g., Hooke 1968). Using MAVEN measurements, Collinson et al. (2019) firstly observed TIDs on Martian ionosphere and they suggested that AGWs are the potential mechanisms. The numerical result displayed in Fig. 14 indicates that AGWs on Mars can possibly contribute to the formation of TIDs at Mars. The physical mechanisms of AGW-driven TIDs are complex and beyond the scope of this paper.

3.4 Wave-associated electric field, current density, and magnetic field

In this section, we investigate the electric field, current density, and magnetic field variations driven by AGWs. As mentioned in Sect. 2.1.2, the ions and electrons move independently in the dynamo region, generating electric currents and a localized polarization electric field. When AGWs propagate into the ionosphere, the plasma motion perturbed by waves can cause variations in the current and electric field. Since the magnetic field in the present Mars is weaker, the wave-driven electric currents can contribute to localized magnetic fluctuations, modifying the ambient magnetic field in the Martian ionosphere.

First, we study the variations in the electric potential and electric field driven by AGWs. Figure 15 displays the amplitudes of the ambipolar, the polarization and the total electric potential driven by AW and GW modes. The perturbed ambipolar electric potential, ϕ_am1, is caused by the electron pressure gradient and can be given by ϕ_am1 = −p_e1/n_e0e. The neutral wind-driven polarization electric potential, ϕ_n1, can be obtained by solving Eq. (37) when the effect of electron pressure gradient is ignored. In this scenario, ϕ_n1 is given by the following formula:

$${\varphi }_{\text{n}1}=\frac{ik\cdot \left[\sigma \cdot \left(V_{\text{n}1}\times B_0\right)\right]}{i\omega {ϵ}_0{k}^2+k\cdot (\sigma \cdot k)},$$(58)

which has been used to investigate the gravity wave-driven Rayleigh-Taylor instability in the F-region of Earth’s ionosphere in earlier works (e.g., Huang et al. 1994; Huang & Kelley 1996a). Here $\mathit{k}\cdot (\mathit{\sigma }\cdot \mathit{k})=k_{\Vert }^2{\sigma }_0+k_{\perp }^2{\sigma }_{\text{p}}$, where k_|| and k_⊥ refer to the parallel and vertical components of the wavenumber along the magnetic field, respectively; σ₀ and σ_p are the field-aligned and Pedersen conductivity.

The total electric potential, ϕ₁, the solution of Eq. (37), is the total effect of the ambipolar and neutral wind-driven polarization electric potential and given by:

$${\varphi }_1=\varpi {\varphi }_{\text{aml}}+{\varphi }_{\text{nl}},$$(59)

where ϖ = σ_0ek_||ξ_e/[iωє₀k² + k⋅(σ⋅k)]. This result suggests that ${\mathit{V}}_{\nabla p_{\text{el}}}^{\Vert }=i{\mu }_{0\text{e}}e{\varphi }_{\text{am}1}{\mathit{\xi }}_{\text{e}}\text{ and }{\mathit{V}}_{{\varphi }_1}^{\Vert }=i{\mu }_{0\text{e}}e{\varphi }_1{\mathit{k}}^{\Vert }$ mentioned in Sect. 3.2 have not perfect cancelation when the neutral-wind drag effects are taken into consideration. In the absence of AGWs, ϖ ≈ σ_0e/σ₀ ≈ 1, indicating that ϕ₁ = ϕ_am1 The electric potential in the ionosphere is dominated by the ambipolar electric potential and ${\mathit{V}}_{\nabla p_{\text{e}1}}^{\Vert }\text{and }{\mathit{V}}_{{\varphi }_1}^{\Vert }$ are completely canceled out in this scenario.

As illustrated in Figs. 15a and b, the wave perturbed ambipolar potential is comparable to the neutral wind-driven polarization electric potential, suggesting that the electron pressure gradient effect is significant and cannot be ignored when waves coupled with the Martian ionosphere. In addition, the maximum amplitude driven by the GW mode is larger than that caused by the AW mode for both ϕ_am1 and ϕ_n1. The reason is that the GW mode has more cumulative effects on the ionosphere compared with the AW mode due to its larger period, leading to higher perturbed electric potential. The perturbed total electric potential, shown in Fig. 15c, has two peak amplitudes since ϕ_am1 and ϕ_n1 have different peak altitudes for each wave mode. The maximum amplitude of ϕ₁ driven by the GW mode can reach up to approximately 4 m V at 180 km altitude and ~2.8 mV at 160 km altitude for the AW mode.

The perturbed ambipolar electric field, neutral wind-driven polarization electric field, and total electric field are displayed in Figs. 16a, b, and c, respectively. We note that the ambipolar electric field only has the field-aligned component and can be calculated by $E_{\text{am}1}^{\parallel }=-i{k}^{\parallel }{\varphi }_{\text{am}1}$. As shown in Fig. 16a, the maximum amplitude of $E_{\text{am}1}^{\parallel }$ driven by the AW and GW modes is approximately 0.22 mV km⁻¹ at 160 km altitude and 0.12 mV km⁻¹ at 180 km altitude, respectively. The reason is that the field-aligned wavenumber of the AW mode is larger than that of the GW mode, leading to a higher perturbed ambipolar electric field. The neutral wind-driven polarization electric field, E_n1 = −ikϕ_n1, is comparable to the perturbed ambipolar electric field. Their peak amplitudes caused by the AW and GW modes are approximately 0.14 mV km⁻¹ and 0.22 mV km⁻¹, respectively. Similar to the total perturbed electric potential, the total electric field, E₁ = −ikϕ₁, also has two peak amplitudes for each wave mode. However, the peak value of E₁ driven by the AW mode can reach to approximately 0.27 mV km⁻¹, which is comparable to the peak amplitude caused by the GW mode. This result indicates that the effect of the perturbed ambipolar electric field is significant and cannot be ignored.

The wave-driven variations in the current density can be given by Eq. (36) and illustrated in Fig. 17. $j_1^{\parallel }$ is the field-aligned current density and their Pedersen and Hall components are denoted by $j_1^{\text{p}}\text{ and }{j}_1^{\text{H}}$, respectively. The peak amplitudes of current density caused by AW and GW mode are approximately 0.37 and 0.36 A km⁻² at the same altitude (~ 150 km), respectively. The three components of j₁ driven by AW mode are comparable. However, as displayed in Fig. 17b, the current density driven by the GW mode is dominated by the field-aligned and Hall components. The reason is that the field-aligned and Hall conductivities are larger than the Pedersen conductivity in the dynamo region of the Martian ionosphere, as shown in Fig. 4c. In addition, the GW mode has a greater cumulative effect due to its larger period, contributing to the formation of the current density along the field-aligned and Hall directions.

The perturbed current density could modify the ambient magnetic field. The wave-induced magnetic field, B₁, satisfies the linear Maxwell-Ampere’s law and Gauss’s law for magnetism, which can be given by:

$$ik\times B_1={\mu }_0{j}_1$$(60a)

and

$$ik\cdot B_1=0,$$(60b)

respectively, where µ₀ refers to the permeability of vacuum. We note that the differential operator ∇ is replaced by ik. Combining the two equations mentioned above, we have:

$$ik\times \left(k\times B_1\right)=i\left(k\cdot B_1\right)-i(k\cdot k)B_1={\mu }_{0}k\times j_1.$$(61)

Hence, the perturbed magnetic field can be solved and given by the following formula:

$$B_1=\frac{i{\mu }_0}{k^2}\left(k\times j_1\right)$$(62)

Figures 18a and b illustrate the perturbed magnetic field driven by the AW and GW modes, respectively. The field-aligned component is marked by $B_1^{\parallel }$. The Hall and Pedersen components are denoted by $B_1^{\text{H}}\text{ and }{B}_1^{\text{p}}$, respectively. The peak ratio amplitudes of the perturbed magnetic field driven by AW and GW modes are ~23% and 31%, respectively, and occur at 150 km altitude for both wave modes. In addition, the peak value of the Hall component is close to the total perturbed magnetic field, whereas the peak amplitudes of the field-aligned and Pedersen components are near 10%. This result suggests that the Hall magnetic field is more sensitive to AGWs. Moreover, even if the peak ratio amplitude of the neutral density is only 10%, the total magnetic field magnitude caused by waves can still reach approximately 30%, suggesting that its high sensitivity to AGWs in the Martian ionosphere.

3.5 Wave-associated thermal variation

In this section, the wave-induced thermal variation in the ionosphere-thermospheric region is investigated. We firstly examine the perturbed pressure and temperature of neutrals and plasmas, which are the direct outputs of the linear coupling system. The perturbed pressure and temperature ratio (marked by p_α1/p_α0 and T_α1/T_α0) driven by AW and GW mode are displayed in Figs. 19 and 20, respectively. The peak amplitude of the perturbed neutral pressure and temperature ratio induced by the AW mode are approximately 13% and 5%, respectively, as illustrated in Fig. 19. We note that the peak neutral density ratio amplitude is fixed at 10%. The maximum amplitude of the AW mode-driven pressure and temperature of electrons are ~I4% and ~ 10%, respectively. These results reveal that the electron pressure and temperature caused by waves are stronger than those of neutrals. The peak amplitudes of the pressure ratio of ${\text{O}}_2^{+}$ and ${\text{CO}}_2^{+}$ driven by AW mode are also greater than that of neutrals and can reach about 17% and 15%, respectively. The pressure variation of O⁺ driven by the AW mode is less than that of neutrals and reaches a maximum of about 8%. However, as shown in Fig. 19b, the temperature variations of O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$ driven by AW mode have the same values, whose peak amplitude is about 12.5% and greater than that of neutrals and electrons.

The neutral pressure and temperature driven by the GW mode, displayed in Fig. 20, reach a maximum of about 10% and 8.5%, respectively. The electrons are also more sensible than neutrals for the GW mode. The electron pressure and temperature peak amplitudes have the same value, which is approximately 18% and greater than that of neutrals. The GW mode-driven pressure ratio peak amplitudes are ~ 12% for both ${\text{O}}_2^{+}$ and ${\text{CO}}_2^{+}$ and greater than that of neutrals. The GW mode-perturbed pressure variation for O⁺ is weaker and reaches a maximum of about 7%. In addition, the temperature variations for both three ion constituents driven by the GW mode also have the same peak amplitude, which is about 9% and greater than that of neutrals.

The results mentioned above indicate that the variations in the ion and electron temperature are greater than neutrals for both AW and GW modes. The temperatures of ions and electrons are more sensible to waves. The pressure variations of electrons and ions driven by AW and GW modes are also stronger than those of neutrals except for O⁺.

Next, we investigate the wave-driven energy variations in the upper atmosphere. When waves propagate into the upper atmosphere from lower levels, the energy and momentum fluxes carried by the vertically propagating waves could be deposited into the background atmosphere. The reason is that the waves can be dissipated by the molecular viscosity, thermal diffusivity, and ion-drag effects in the upper atmosphere. The total neutral wave heating/cooling rate associated with wave dissipation can be calculated by using the following formula (e.g., Hickey et al. 2000; Wang et al. 2020):

$$H_{\text{wave }}=H_{\text{vis }}+H_{\text{sen }}+H_{\text{press }}+H_{\text{Euler }}$$(63)

where H_vis = − < σ_m1 : ∇_n1 > is the viscous heating rate; H_sen = −d{c_pρ_n0 < w₁T_n1 >} /dɀ refers to the sensible heat flux divergence; H_pгess =< V_nı · ∇p_n1 > represents the work done per unit time by the wave-induced pressure gradient; H_Euler = − < w_n1ρ_n1 > 𝑔 is the work done per unit time by the second-order wave-induced Eulerian drift. The angled bracket ‘<>’ denotes a operator that can calculate a time-average of the product of physical quantities. The calculation expression is given by < AB >= Re(AB*)/2, where A and B are arbitrary complex functions; subscript ‘*’ refers to the complex conjugate; ‘Re’ denotes the real part of the complex value. The total neutral wave energy heat flux can be given by $F_{\text{tot }}=-{{{\displaystyle \int }}^{\text{}}}_z^{\infty }\left(H_{\text{vis }}+H_{\text{sen }}+H_{\text{press }}+H_{\text{Euler }}\right)\text{d}z$.

In addition, the energy exchange between neutral particles and ion constituents or electrons exists in the ionosphere-thermospheric region, suggesting that neutrals can also be heated or cooled in the neutral-plasma collisional processes. The neutral collisional heating rate in the presence of AGWs can be calculated by $H_{\text{n}\alpha }={\overline{Q}}_{\text{n}\alpha 1}$, where Q_nα refers to the neutral-plasma collisional heating and has given by Eq. (8); subscript ‘1’ represents the values in the wave perturbed state; superscript ‘-’ denotes a time-average. The photochemical heating of the hot oxygen atom in the presence of AGWs is also investigated and given by $H_{\text{pc}}={\overline{R}}_{\text{O}}{\displaystyle \sum _j{E}_j}{L}_j$, where subscript ‘j’ refers to the j-th channel of ${\text{O}}_2^{+}$ dissociative recombination (marked by R₁ in Table A.1); E_j and L_j represent the associated exothermicity and the branching ratio of the jth channel; ${\overline{R}}_{\text{O}}$ is the wave-induced hot oxygen atom production rate and is given by ${\overline{R}}_{\text{O}}=<n_{\text{e}1}{n}_{{\text{O}}_2^{+}}^{\prime }>{\alpha }_1$ with α₁ being the rate coefficient of R₁.

Ions and electrons can also be heated or cooled during plasma-neutral collisional processes. The wave-driven plasma collisional heating rate can be given by $H_{\alpha n}={\overline{Q}}_{\alpha n1}$. We quantified the ion heating/cooling rate and electron heating/cooling rate resulting from collisions with neutrals, denoted by H_in1 and H_en1, respectively. In addition, the perturbed current density and electric field generate Joule heating, whose rate can be calculated by H_joule =< j₁ ⋅ E₁ >.

Figuree 21 displays the neutral, ion and electron heating/cooling rate profiles driven by AW and GW mode. As shown in Figs. 21a and c, the total neutral wave heating/cooling rates, H_wave, induced by AW and GW mode can reach a maximum of about 4.5×10³ eV cm⁻³ s⁻¹ near 140 km height and 3.5×10³ eV cm⁻³ s⁻¹ near 135 km height, respectively. The wave heating/cooling rates (degrees per unit time), Q = H_wave/ρ_n0c_p (e.g., Hickey et al. 2000; Schubert et al. 2003; Hickey et al. 2011), for AW and GW wave modes are also calculated and displayed in panels a and c of Fig. 21, respectively. The wave heating rate (degrees per unit time), Q, induced by AW and GW mode can reach a maximum of about 3.29×10² K sol⁻¹ near 164km height and about 97 K sol⁻¹ near 146 km height, respectively (here ‘sol’ represents one Martian day). These results are of a similar order with those calculated by Parish et al. (2009) and Walterscheid et al. (2013) with a 1-D linear full-wave model and by Medvedev & Yiğit (2012) with a 3-D nonlinear general circulation modeling for the Martian atmosphere.

The wave-driven photochemical heating of the hot oxygen atom, H_pc, is dominated, whose peak rates are approximately 10⁴ eV cm⁻³ s⁻¹ for both the AW and GW modes, as illustrated in Figs. 21a and c. However, the neutral collisional heating rate during neutral-plasma collisional process is weaker. The neutral-ion collisional heating rate, H_ni, has a maximum value of ~2 eV cm⁻³ s⁻¹ and ~2.5 eV cm⁻³ s⁻¹ for AW and GW mode, respectively. While the AW mode-driven and GW mode-induced neutral-electron collisional heating rate, marked by H_ne, reach maximum of about 3.5 ×10⁻² eV cm⁻³ s⁻¹ and 3 × 10⁻² eV cm⁻³ s⁻¹. Below 200 km height, the neutral-ion collisions are more frequent than neutral–electron collisions, meaning that H_ni is greater than H_ne. However, the neutral-plasma collisional heating is weakest compared with the total neutral wave heating and the photochemical heating, indicating the collisional heating is insignificant when waves propagating into the ionosphere-thermospheric region.

Since the ionosphere-thermospheric system experiences energy exchange between neutrals and plasmas, neutrals are heated while ions or electrons are cooled during the collisional processes. As illustrated in Figs. 21b and d, the ion cooling rate, H_in, reach a maximum of about 2 eV cm⁻³ s⁻¹ for both the AW and GW modes, which is comparable to the maximum of neutral-ion collisional heating rate. While the maximum of electron cooling rate is approximately 10 eV cm⁻³ s⁻¹ for both two wave modes, which is greater than the neutral-electron heating rate due to the lighter electron mass. In addition, the Joule heating rate generated by the perturbed current density and electric field, marked by H_joule, is less than the plasma-neutral collisional cooling rate. the maximum H_joule driven by AW and GW mode are 4 × 10⁻² eV cm⁻³ s⁻¹ and 6 × 10⁻² eV cm⁻³ s⁻¹, respectively. Matta et al. (2014) suggests that the solar heating rate for electrons and ions below 200 km height in the Martian ionosphere have a maximum of 3.5×10³ eV cm⁻³ s⁻¹ and 5 eV cm⁻³ s⁻¹, respectively. This result reveals that the wave-driven Joule heating is inconsiderable compared with the solar heating mechanism.

Correspondingly, the neutral, ion and electron energy flux profiles driven by AW and GW mode are displayed in Fig. 22. The energy flux is obtained by downwardly integrating the heating/cooling rate along altitude. The maximum photochemical energy fluxes for both two wave modes are about 3.4 × 10¹⁰ eV cm⁻² s⁻¹ , which is greater than the total wave energy flux and neutral collisional energy flux. The ion-neutral collisional energy flux is comparable to the neutral–ion collisional energy flux but in the opposite direction, and its value is close to approximately 10⁷ eV cm⁻² s⁻¹ for both two wave modes. However, the electron-neutral collisional energy flux is two orders of magnitude greater than that of the neutral-electron. The Joule heat flux can reach a maximum of about 1.4 × 10⁵ eV cm⁻² s⁻¹ for both two wave modes. In addition, we also investigate the wave-driven Poynting flux caused by the perturbed electric field and magnetic field, which is calculated by F_pt =< E₁ × B₁ > /µ₀. The Poynting flux is comparable to the Joule heating flux for both two wave modes and its value is in the order of 10⁵ eV cm⁻² s⁻¹, but is still less than the plasma-neutral collisional energy flux. This result indicates that the wave-driven Joule heating and Poynting flux are insignificant compared with the perturbed plasma-neutral collisional effect.

4 Comparison with spacecraft observations and previous modeling works

To validate our wave model, it is necessary to compare our numerical results to the spacecraft observations and previous modeling works. Earlier observations revealed that AGWs have a wide spatial range of ten to about 300 km with an amplitude of tens of percent (e.g., Creasey et al. 2006; Terada et al. 2017; Siddle et al. 2019). Near the Martian surface, AGWs can be triggered by flow-over topography, weather systems, and convection (e.g., Parish et al. 2009; Heavens et al. 2020; Yiǧit et al. 2021b). These waves typically have a variety of wave frequencies and phase speeds (e.g., Fritts & Alexander 2003; Heavens et al. 2020). To represent the realistic behaviors of waves observed by the spacecraft, we restrict the horizontal wavelength, λ_h, to range from 10 km to 300 km. Three wave sets for each acoustic branch and gravity branch, marked by “AWs” and “GWs", respectively, are selected, whose vertical wavenumbers at the lower boundary $\left(k_z^l\right)$ are 0.11, 0.32, and 0.99 km⁻¹ and −0.07, −0.19, and −0.5 km⁻¹, respectively. The λ_h − c_p relationship of the selected wave sets is illustrated in Fig. 6. The frequency and phase speed of the selected wave modes range from lower to higher frequencies and about 10 to 1000 m s⁻¹, respectively, suggesting that these selected waves have the capability to represent the realistic behaviors of waves generated near the Martian surface. To coincide with the spacecraft measurement, the maximum amplitude of the perturbed neutral density ratio for the selected wave modes is fixed at 10% (e.g., Terada et al. 2017). We also assume that the selected waves propagate in the x-ɀ plane. Correspondingly, the wave azimuth angle is equal to 90°. The present Mars lacks a strong dipole magnetic field. There is a crustal magnetic field region mainly located in the Southern hemisphere (Acuna et al. 1998). Out of the crustal magnetic field region, the magnetic field is dominated by the horizontal draped interplanetary magnetic field. In this study, we only focus on the wave propagation outside the crustal magnetic field region. Therefore, the dip angle of the magnetic field is assumed to be 0°.

The wave amplitude is a function of altitude and has a maximum value in a damping altitude due to wave dissipation caused by viscosity and thermal conduction. Therefore, we could only focus on the damping altitude distribution of the selected wave modes. The damping altitudes (ɀ_damping) of AWs and GWs are plotted in Figs. 23a and c, respectively, and range from ~80 km to 170 km for selected wave modes. These damping altitudes are very good agreement with those calculated by the linear wave numerical simulations of Parish et al. (2009) and Walterscheid et al. (2013) and nonlinear general circulation modeling of Yiǧit et al. (2015a). The vertical wavelengths at the damping altitude, λ_ɀ(ɀ_damping), of AWs and GWs are illustrated in Figs. 23b and d, respectively and range from about 8 km to 160 km. To further compare our model results to previous modeling works, the neutral thermal effect caused by waves is studied in this paper. Figures 24a and d display peak altitude and peak magnitude distributions of wave heating/cooling rates driven by AWs, respectively. The result suggests that the acoustic wave branches heat the background atmosphere at all heights, similar to the scenario of Jupiter (e.g., Schubert et al. 2003) and Earth (e.g., Hickey et al. 2001; Schubert et al. 2005; Yiǧit & Medvedev 2009). The maximum of wave heating rate, Q_h, driven by AWs ranges from about 60 to about 540 K sol⁻¹ and occurs at a region whose height ranges from about 80 to 180 km heights. These results are comparable to those calculated by Walterscheid et al. (2013) with a linear full wave model, whose maximum heating rate ranges from about 145 to 410 K day⁻¹ occurred at 60–250 km heights. Gravity wave branches not only heat the upper atmosphere but also cool the upper atmosphere through the effects of sensible heat flux divergence, which are very different from acoustic wave branches (Schubert et al. 2005, and references therein). The peak altitude and peak magnitude distributions of wave heating/(cooling) rates, Q_h/(Q_c), driven GWs are illustrated in panels b/c and e/f of Fig. 24, respectively. The maximum wave heating rate induced by GWs ranges from ~10 to ~230 K sol⁻¹ occurred at approximately 80–250 km heights. The GWs-driven maximum wave cooling rate is up to about 0 to 220 K sol⁻¹ and occurred at 100-170 km heights. This heating/cooling rate driven by gravity wave branches is in very good agreement with the results calculated by earlier works, including the linear full wave model of Parish et al. (2009) and Walterscheid et al. (2013), and the nonlinear general circulation modeling of Medvedev & Yiǧit (2012). The comparison indicates that the first part of our coupled model, describing the neutral behaviors, is effective and consistent with previous modeling work. Next, we compare the results of the second part of our coupled model, which describes the plasma behaviors, with the spacecraft observations.

Recently, highly variable ionospheric fluctuations at Mars have been observed by spacecraft measurements, such as ionospheric irregularities (e.g., Fowler et al. 2017, 2019) and TIDs (Collinson et al. 2019). They argued that AGWs arising from the lower regions are the potential seeding mechanisms of the these variable structures due to their comparable temporal and spatial distributions. In addition, the ion density variations caused by AGWs in the thermosphere-ionospheric region on Mars are directly observed by MAVEN in situ measurements (Leelavathi et al. 2023). In this section, we investigate the variations in the plasma density and magnetic field by selecting a series of wave modes and comparing them with the spacecraft observations.

We first examine the wave-driven plasma density variations for the selected wave modes. We note that the effects of photochemical reactions, electron pressure gradient, and electric field are considered. Similarly, the wave-driven plasma density also has a maximum amplitude at a peak altitude. Figures 25 and 26 display peak altitude and peak amplitude distributions of plasma density ratio (n_α1/n_α0) driven by AWs and GWs, respectively. The first four panels of Figs. 25 and 26 display the peak altitudes of the perturbed plasma density ratio (including electron, O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$) driven by the AWs and GWs and their values range from approximately 80 km to 180 km and 80 km to 230 km, respectively. Figures 25e-h illustrate the AWs-driven peak density ratio amplitudes of electron, O+, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$, with values ranging from ~1% to 10.5%, 0.1% to 5.5%, 4.5% to 12%, and 0.1% to 6%, respectively. The GWs-driven peak density ratio amplitude ranges from ~5% to 80% for both electron and ${\text{O}}_2^{+}$, while ranges from 0.25% to 5.5% for O⁺ and 0.3% to 7% for ${\text{CO}}_2^{+}$, as illustrated in Figs. 26e–h. The result suggests that the wave-driven density variations in electrons are consistent with that of ${\text{O}}_2^{+}$, while the variations in O⁺ density in accordance with those in ${\text{CO}}_2^{+}$ for both AWs and GWs.

The peak altitude and peak amplitude distributions of perturbed magnetic field ratio (B₁ /B₀) are illustrated in Fig. 27. The peak values of B₁ / B₀ occurs at a region whose height ranging from approximately 105 km to 160 km and 114 km to 158 km for AWs and GWs, as shown in Figs. 27a and c, respectively. The peak amplitude of B₁ / B₀ can reach a maximum value of about 75% and 40% for the longer-wavelength AWs and GWs while its value is much less than 1% for the wave modes with shorter wavelength, as shown in Figs. 27b and d.

Based on in situ measurements by MAVEN, a statistical study carried out by Fowler et al. (2019) reveals that the small-scale ionospheric irregularities may be driven by gravity waves. The statistical results suggest that the vertical length scales of the irregular structures range from approximately 5 km to 20 km. The magnitude variations in the magnetic field associated with these irregularities are about 5–15%, while plasma density variations range from ~10% to 75% and occur between heights of ~130 and 160 km. In addition, the statistical study showed that these irregularities occur outside the region of the crustal magnetic field, which is consistent with the horizontal magnetic field assumption mentioned above. To compare with the spacecraft measurements, we display the plasma density and magnetic field variations driven by the waves with λ_ɀ(ɀ_damping) ranging from ~5 km to ~30 km, whose $k_z^l$ equal to 0.32 km⁻¹ and 0.99 km⁻¹ for acoustic branch and −0.19 km⁻¹ and −0.5 km⁻¹ for gravity branch. The variations in the plasma density and magnetic field driven by these waves are summarized in Table 2. The peak density and magnetic field magnitude ratio amplitude of the major ion constituents, ${\text{O}}_2^{+}$, driven by gravity branch range from 10% to 80% occurred at 80–230 km height and 0.2% to 17% occurred at 114–158 km height, respectively. The numerical results of our model are consistent with the plasma density and magnetic field variations associated with the small-scale ionospheric irregularities observed by Fowler et al. (2019), indicating that AGWs are indeed the seeding mechanism of the ionospheric irregularities.

In addition, the plasma density average variation associated with TIDs observed by Collinson et al. (2019) is about 35.2%, with a mean wavelength of ~93 km. As illustrated by Figs. 26e and g, the peak amplitude of electron and ${\text{O}}_2^{+}$ driven by the wave mode with $k_z^l=-0.5$ km⁻¹ and λ_h = 100 km are about 25% and 30%, respectively, which agree with the plasma density variation associated with TIDs observed by Collinson et al. (2019). Therefore, AGWs also have the potential to drive the formation of TIDs at Mars.

Moreover, a statistical result from Leelavathi et al. (2023) based on MAVEN in situ measurements indicates that the ion density amplitudes can reach about 70% with a average value of 5%. The numerical results suggest that the peak density amplitudes of electron, O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$ driven by AGWs range from several to tens of percent, as displayed in Figs. 25 and 26. The electron and ${\text{O}}_2^{+}$ density amplitude driven by the shorter-wavelength gravity wave modes can reach about 80%. In addition, Leelavathi et al. (2023) reveals that the wave-driven ion density variations are more on the nightside than those on the dayside. The wave-driven plasma density variations in the absence or presence of photochemistry have been investigated by considering three scenarios, which described in Sect. 3.3. The results indicate that the wave-driven plasma density variations is weaker when the effect of photochemistry is considered, suggesting that the plasma density variations driven by AGWs on the dayside is less than those on the nightside. The observational results carried out by Leelavathi et al. (2023) also indicate that the ion density amplitudes are moderately correlated with those in the neutrals. We illustrate the relationship between the neutral peak density amplitude and plasma peak density amplitude driven by AW and GW wave modes in Figs. 28a and b, respectively. The density variations of electron, O⁺, ${\text{O}}_2^{+}$, ${\text{CO}}_2^{+}$ increase with the increasing neutral amplitude for both AW and GW wave modes. The comparison mentioned above indicates that the numerical results are in very good agreement with the spacecraft observations and verified that our model described in this study is effective.

5 Discussion and conclusions

Existing spacecraft measurements suggest that AGWs may be a potential mechanism of the highly variable ionospheric structure at Mars (e.g., Gurnett et al. 2010; Harada et al. 2018; Fowler et al. 2017, 2019, 2020; Collinson et al. 2019, 2020; Wan et al. 2024a,b). Although the wavelike oscillations associated with these variable ionospheric structures have been observed, the coupled physical process between AGWs and ionosphere remain unexplored. Earlier works proposed a linearized wave model adopting a WKB approximate method to investigate the plasma density variations in the Jovian and Saturnian ionosphere (Matcheva et al. 2001; Barrow & Matcheva 2011, 2013). Wang et al. (2023) only examine the electron density variations in the Martian ionosphere using the WKB wave model, while the wave-induced ion density oscillations, effects of photochemistry, thermal variations, and electromagnetic field at the ionosphere of Mars are not taken into consideration in the previous works. In the present study, we construct a comprehensive coupled model to investigate the wave-driven variations in the plasma density, temperature, electromagnetic field, and thermal structures in the presence of photochemistry at Mars.

The governing equations of the comprehensive coupled model include two parts. The first part describes the behavior of neutral atmosphere. To characterize the wave propagation within a realistic atmosphere, we adopted a full-wave model rather than a WKB wave model as the first part of our coupled model. The WKB wave model is valid in an isothermal atmosphere. We examined the response of ionosphere to the waves propagating from lower levels to upper regions, in which the temperature is varied with height. The full-wave model, a time-independent linear model that describes the propagation of waves within nonhydrostatic and inhomogeneous atmosphere, has been extensively used to study the atmosphere of Earth (e.g., Hickey et al. 1998), Jupiter (e.g., Hickey et al. 2000), and Mars (Parish et al. 2009; Walterscheid et al. 2013; Huang et al. 2022a). The term “full wave” suggests that this model fully accounts for the effects of wave reflection, Coriolis force, molecular and eddy diffusion, thermal conductivity, ion-drag, and the varying background temperature with height. We improved the full-wave model in our work by adding the collisional force acting on the neutral particles, the neutral–plasma collisional heating, and photochemical heating.

The second part of our coupled model consists of the governing equations describing the plasma behaviors. Earlier studies of the impact of neutral waves on the ionosphere only assumed that the wave-driven plasma velocity is equal to the neutral wave speed projected on the magnetic field lines and calculated the perturbed plasma density by solving the plasma continuity equation in the absence of the photochemical production rate (e.g., Hooke 1968; Matcheva et al. 2001; Barrow & Matcheva 2011; Wang et al. 2023). The effects of the plasma pressure gradient and electromagnetic field are neglected in previous works but are fully taken into consideration in our coupled model of this study. The reason is that Mars presently has a weak magnetic field, meaning that the plasma pressure is comparable to the magnetic pressure. In addition, the wave-like oscillations in the ambient magnetic field have been observed in the Martian ionosphere (Collinson et al. 2019; Fowler et al. 2019). The presence of energy exchange between neutrals and plasmas motivated us to include the plasma-neutral collisional heating in the second part of the coupled model.

The model output consists of the perturbed velocity, density, temperature, and pressure of neutrals and plasmas and wave-associated electric field, current density, and magnetic field. The wave-associated plasma motion can be revealed by the perturbed plasma velocity. Our results indicate that ion constituents are mainly dominated by the neutral-drag effects, meaning that the ions oscillate along the wave phase lines. The wave-driven electron velocity consists of the field-aligned neutral wind speed and the perturbed velocity induced by the neutral wave-driven polarization electric field along the magnetic field lines. Therefore, electrons are primarily subject to the gyro-motion along the magnetic field, which leads to the formation of plasma compressional and rarefactional regions in the Martian ionosphere. In addition, the wave-driven plasma motion can generate the perturbed electric currents and localized polarization electric field, further contributing to the localized magnetic perturbation. Our results indicate that the wave-driven variations in the magnetic field magnitude can reach tens of percent. We also investigated the wave-associated ambipolar electric field driven by the electron pressure gradient and find it to be comparable to the neutral-driven polarization electric field, suggesting that the effect of the electron pressure gradient is significant and cannot be ignored.

We also examined the wave-driven density variations of electron and major ion constituents (including O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$) associated with the photochemistry. The model results reveal that the plasma density variations when the photochemical reactions are not considered are several times larger than those in the presence of photochemistry and larger than the neutral density variations. In the presence of photochemistry, the wave-driven variations of the electron/ion loss and production rate could result in a decrease in the density perturbations of the charge species. We also find that plasma density variations are greater when the wave period is larger; this is because the ions or electrons have more cumulative effects when they respond to the waves with larger periods.

The wave-associated thermal variations are also investigated. The waves can be dissipated in the upper atmosphere and release the momentum fluxes and energy, which lead to the heating and cooling of the background atmosphere. Our results suggest that the acoustic wave branches heat the background atmosphere at all heights, which is similar to the scenario of the Jovian atmosphere and that of Earth (e.g., Hickey et al. 2001; Schubert et al. 2003, 2005; Yiǧit & Medvedev 2009). The gravity wave branches not only heat but also cool the background atmosphere due to the effects of sensible heat flux divergence (Schubert et al. 2005). In addition, the energy exchange exists in the ionosphere–thermospheric region, meaning that the neutrals can be heated by the collisions with ions or electrons. Our numerical results indicate that the wave-driven neutral–ion col-lisional heating rate is greater than neutral–electron collisional heating rate since the neutral–ion collisions are more frequent than neutral–electron collisions. We also examined the photochemical heating of the hot oxygen atom driven by waves, finding that the wave-associated neutral–plasma collisional heating is less than the neutral-wave heating and the wave-driven photochemical heating. The heat sources of the neutral atmosphere are dominated by the photochemical processes and the neutral–plasma collisional heating is insignificant when waves propagate into the ionosphere–thermospheric region.

In the collisional process, neutral species are heated but ions and electrons are cooled. The model results suggest that the ion–neutral collisional cooling rate is comparable to the neutral–ion collisional heating rate. The electron-neutral collisional cooling rate is greater than the neutral–electron heating rate due to the lighter electron mass. In addition, the electrons are cooler than ions in the plasma-neutral collisional process. The wave-driven Joule heating and Poynting energy flux generated by the perturbed current density and electromagnetic field are also examined. The numerical results reveal that the Poynting flux is comparable to the Joule heating flux but is weaker than the plasma-neutral collisional heating flux.

We compared the numerical results to the spacecraft observations and previous modeling works in order to validate our coupled model. First, we investigated the neutral thermal effect caused by waves to compare to the previous modeling studies. Our numerical results indicate that the maximum heating rate caused by the acoustic wave branches ranges from approximately 60 to about 540 Kelvin per sol, which is consistent with the previous model results calculated by Walterscheid et al. (2013). The maximum heating rate induced by the gravity wave branches ranges from about 10 to 230 Kelvin per Martian day. While the maximum cooling rate driven by the gravity wave branches ranges from about 0 to 220 Kelvin per Martian day. Our numerical results regarding the wave heating and cooling rate are consistent with the results calculated by the linear full wave model of Parish et al. (2009) and Walterscheid et al. (2013), and the nonlinear general circulation modeling of Medvedev & Yiǧit (2012).

Second, we compared the plasma variations calculated with our model to the spacecraft observations. Highly variable ionospheric structures at Mars were recently observed (e.g., Gurnett et al. 2010; Harada et al. 2018; Fowler et al. 2017, 2019, 2020; Tian et al. 2022; Collinson et al. 2019, 2020; Wan et al. 2024a,b). The observational results suggest that AGWs are a potential seeding mechanism of the variable ionospheric structure due to their comparable temporal and spatial distributions. The observations carried out by Fowler et al. (2019) led the authors to argue that the plasma density and magnetic field variations associated with the small-scale ionospheric irregularities are ~ 10–75% and ~5–15%, respectively. The vertical length scale of the irregularities ranges from ~5 km to 20 km. We specifically selected several typical wave modes whose vertical wavelengths at damping altitude ranging from about 5 km to 30 km to examine the variations in the plasma density and magnetic field magnitude. The amplitudes of major ion density and magnetic field driven by the selected waves range from about 10% to 80% and ~0.2% to 17%, respectively. Collinson et al. (2019) studied TIDs using MAVEN measurements and suggested that the plasma density variation associated with TIDs is about 35.2% with a wavelength of ~93 km. Our model results reveal that the peak amplitude of electrons and ${\text{O}}_2^{+}$ induced by a gravity wave mode with a horizontal wavelength of 100 km are approximately 25% and 30%, respectively. Moreover, the plasma density variations associated with AGWs were directly measured by MAVEN (Leelavathi et al. 2023). These authors suggested that the ion density amplitudes can reach as high as 70% with an average value of 5%. Our model result indicates that the density amplitudes of electrons, O+, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$ driven by AGWs range from several to tens of percent all the way up to 80% for the shorter wavelengths. Leelavathi et al. (2023) also suggests that the wave-driven variations in the plasma density are more on the nightside than on the dayside and density amplitudes in ions are about two times greater than those in neutrals. The photochemical reactions occur in the dayside ionosphere of Mars. The model results indicate that the plasma density amplitudes in the presence of the photochemistry are less than those in the absence of the photochemistry, suggesting that the plasma variations are weaker on the dayside. Additionally, we also find that the plasma density amplitudes are several times greater than the neutral density amplitudes when the photochemical effects are not taken into consideration. Observations made by Leelavathi et al. (2023) also suggest that the density amplitudes in ions show moderate correlation with those in neutrals. We examined the peak amplitude relationship between plasma and neutrals and find that the density variations of electrons, O+, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$ increase with the increasing neutral amplitudes. The comparison mentioned above indicates that our model results are in excellent agreement with spacecraft observations and verify that the coupled model is effective. AGWs indeed play a significant role in modifying the ionospheric structure at Mars.

The coupled model in this study only considers the linear wave process in the upper atmosphere. Nonlinear behaviors, such as wave breaking and saturation, frequently occur in the real atmosphere. Waves propagating into the thermosphere can be dissipated not only by molecular viscosity and thermal conduction but also by the nonlinear processes (e.g., Fritts 1984; Franke & Robinson 1999; Fritts & Alexander 2003; Yiǧit et al. 2008). Recent studies suggest that nonlinear wave dissipation plays an important role in the influence of gravity waves on the formation of high-altitude Martian CO₂ ice clouds and thermospheric gravity wave activity observed by MAVEN during solar minimum (Yiǧit et al. 2018, 2021b). These results indicate that the nonlinear wave dynamics may play a crucial role in the wave-ionospheric coupling process. It is necessary to investigate the impact of nonlinear wave behaviors on the AGW-ionospheric coupling in future work.

In this paper, we propose a comprehensive coupled model to investigate the ionospheric response to neutral waves. Earlier works only consider the effect of neutral waves on thermo-spheric variations (Parish et al. 2009; Medvedev et al. 2011, 2013; Medvedev & Yiǧit 2012; Walterscheid et al. 2013; Yiǧit et al. 2015a; Kuroda et al. 2015, 2016, 2019, 2020; Imamura et al. 2016; Huang et al. 2022a; Roeten et al. 2022; Srivastava et al. 2022). Our model can be used to examine the wave-driven velocity, density, temperature, and pressure of both neutrals and plasmas. Additionally, the variations in the wave-associated electromagnetic field and thermal structure can also be obtained by our model. The wave heat source effects are included in our coupled model. Therefore, our model can be applied to study the ionospheric response to the waves generated by many mechanisms, such as dust storms near the Martian surface. The gravity wave activity during a dust storm in the Martian ther-mosphere was recently observed and simulated (e.g., Medvedev et al. 2013; Kuroda et al. 2020; Leelavathi et al. 2020; Yiǧit et al. 2021; Shaposhnikov et al. 2022; Wu et al. 2022), but the impact of the dust storm-enhanced gravity waves on the ionosphere remains undetected and needs to studied in future work. Our model also includes the background magnetic field, suggesting that it can be applied to examine the wave-driven plasma motions in the crustal magnetic field region. MAVEN observations revealed that the correlation between density amplitudes in neutrals and ions is weaker when waves propagate in a strong crustal magnetic field region (Leelavathi et al. 2023). In this study, we only examine the ionosphere response to waves outside the crustal magnetic field region by assuming that the dip angle of the magnetic field is equal to zero. The wave-induced plasma behavior associated with the Martian crustal magnetic field is an interesting topic and should be investigated in future work. Our model can be used to examine whether AGWs are the seeding mechanisms of the highly variable ionospheric structures at Mars. Although various structures have been observed in the lower and middle ionosphere of Mars, such as sporadic E-like layers and rifts and plasma depletions (Collinson et al. 2020; Basuvaraj et al. 2022), the underlying processes are unknown. Additionally, earlier works studied the ionospheric irregularities using theoretical and numerical methods (Keskinen 2018; Jiang et al. 2021, 2022), but the effects of AGWs are ignored in their works. We plan to investigate the underlying coupling processes between the variable ionospheric structures and neutral waves using our coupled model in future work. In addition, our model can be applied to other planets or satellites with atmosphere, such as Earth, Venus, Jupiter, Saturn, and Titan, in order to further investigate how AGWs interact with the ionosphere in these celestial bodies.

Acknowledgements

We acknowledge the MAVEN team for providing data. MAVEN data used for this study are publicly available at the NASA Planetary Data System (https://pds-atmospheres.nmsu.edu/data_and_services/atmospheres_data/MAVEN/maven_main.html). This work is supported by the National Natural Science Foundation of China (NSFC) under grants 42241112 and 42122061, the Science and Technology Development Fund, Macau SAR (File No.0003/2022/AFJ and 0098/2022/A2), and by the Specialized Research Fund for State Key Laboratories.

Appendix A Plasma densities in undisturbed and perturbed state

Plasma is predominantly produced by the photoionization by absorption of solar extreme ultraviolet (EUV) radiation in the Martian dayside ionosphere. The chemical reactions of the major charged species (O⁺, ${\text{O}}_2^{+}$, and ${\text{CO}}_2^{+}$) in the ionosphere are summarized in Table A.1. Electrons are mainly produced by the photoionization reactions L₁ and L₂ and mainly lost by the recombination reactions R₁. Under the quasi-neutrality condition, the electron loss rate can be assumed to be equal to $L_{\text{e}}\approx {\alpha }_1{n}_{\text{e}}^2$ with α₁ being the coefficient of reaction R₁ and n_e being the undisturbed electron number density. Due to the chemical equilibrium in the lower-altitude ionosphere, we can assume that the plasma production rate equals the loss rate. Therefore, the electron number density in the dayside of ionosphere can be given by $n_{\text{e}}=\sqrt{P_{\text{e}}/{\alpha }_1}$, where $P_{\text{e}}=P_{ph}^O+P_{ph}^{C{O}_2}$ refers to the electron production rate. $P_{ph}^O\text{ and }{P}_{ph}^{C{O}_2}$ are the photoionization rates of O and CO₂ and can be calculated by the following:

$$P_{ph}^j(z)={\displaystyle {\int }_0^{{\lambda }_{\max }}{\sigma }_j}(\lambda )n_j(z)F(z,\lambda )d\lambda ,$$(A.1)

where λ denotes the wavelength of solar EUV and X-ray; σ_j is the absorption cross-section of j-th neutral species; F represents the solar energy flux at altitude ɀ and given by F(ɀ, λ) = F_∞(λ)exp(−τ/cos χ). Here F_∞ is the solar energy flux near the topside atmosphere of Mars. χ is the solar zenith angle $(\text{SZA})\cdot \tau (ɀ,\lambda )={\displaystyle \sum {\sigma }_j}(\lambda ){\displaystyle {\int }_{ɀ}^{\infty }{n}_j}(ɀ)dɀ$ is the optic depth.

The production and loss rates of the major ion species are summarized in Table A.2. The number densities of O⁺, ${\text{O}}_2^{+}{\text{, and CO}}_2^{+}$ in the absence of acoustic-gravity waves can be solved by the following matrix equation under the chemical equilibrium assumption:

$$\left(\begin{array}{ccc}{k}_1{n}_{{\text{CO}}_2}& 0& -k_3{n}_{\text{O}}\\ -k_1{n}_{{\text{CO}}_2}& {\alpha }_1{n}_{\text{e}}& -k_2{n}_{\text{O}}\\ 0& 0& {\alpha }_2{n}_{\text{e}}+\left(k_2+k_3\right)n_{\text{O}}\end{array}\right)\left(\begin{array}{c}{n}_{{\text{O}}^{+}}\\ n_{{\text{O}}_2^{+}}\\ n_{{\text{CO}}_2^{+}}\end{array}\right)=\left(\begin{array}{c}{P}_{ph}^O\\ 0\\ P_{ph}^{C{O}_2}\end{array}\right),$$(A.2)

where k₁ , k₂, and k₃ are the reaction rate of O-CO₂ charge exchange reaction C₁ , C₂, and C₃, respectively; α₁ and α₂ are the recombination reaction rate coefficient of R₁ and R₂ (see Table A.1); n_O and ${\text{n}}_{{\text{CO}}_2}$ are the number density of O and CO₂.

In the presence of acoustic-gravity waves, the perturbed density of three ion constituents associated with photochemistry can be obtained from Eq. (57). In this situation, the perturbed number densities of O+, ${\text{O}}_2^{+}{\text{, and CO}}_2^{+}$ can be solved by the following linear system:

$$G{n}_{\text{ions }}^{\prime }=f,$$(A.3)

where ${\mathit{n}}_{\text{ions}}^{\prime }={\left(n_{{\text{O}}^{+}}^{\prime },n_{{\text{O}}_2^{+}}^{\prime },n_{{\text{CO}}_2^{+}}^{\prime }\right)}^T$ is the perturbed density variable of three major ion constituents. G and f are the coefficient matrices and are given by the following:

$$G=\left(\begin{array}{ccc}-i\omega +k_1{n}_{{\text{CO}}_2}& 0& -k_3{n}_{\text{O}}\\ -k_1{n}_{{\text{CO}}_2}& -i\omega +{\alpha }_1{n}_{\text{e}0}& -k_2{n}_{\text{O}}\\ 0& 0& -i\omega +{\alpha }_2{n}_{\text{e}0}+\left(k_2+k_3\right)n_{\text{O}}\end{array}\right)\text{, }$$(A.4a)

and

$$f=\left(\begin{array}{c}-{\aleph }_{{\text{O}}^{+}}{n}_{{\text{O}}^{+}}+{\tilde{P}}_{ph}^O+k_3{n^{\prime }}_{\text{O}}{n}_{{\text{CO}}_2^{+}}-k_1{n}_{{\text{O}}^{+}}{n^{\prime }}_{{\text{CO}}_2}\\ -{\aleph }_{{\text{O}}_2^{+}}{n}_{{\text{O}}_2^{+}}+k_1{n}_{{\text{O}}^{+}}{n^{\prime }}_{{\text{CO}}_2}+k_2{n^{\prime }}_{{\text{CO}}_2^{+}}{n^{\prime }}_{\text{O}}-{\alpha }_1{n}_{{\text{O}}_2^{+}}{n}_{\text{e}1}\\ -{\aleph }_{{\text{CO}}_2^{+}}{n}_{{\text{CO}}_2^{+}}+{\tilde{P}}_{ph}^{C{O}_2}-{\alpha }_2{n}_{\text{e}1}{n}_{{\text{CO}}_2^{+}}-\left(k_2+k_3\right){n^{\prime }}_{\text{O}}{n}_{{\text{CO}}_2^{+}}\end{array}\right),$$(A.4b)

respectively, where subscript ‘′’ refers to the perturbed value; ℵ_i = ik · V_i1 − (V_i1 · e_ɀ)/H_i with V_i1 being the perturbed ion velocity and H_i being the ion density scale height. n_e1 is given by the solution of Eq. (56). ${\tilde{P}}_{ph}^O\text{\hspace{0.17em}and\hspace{0.17em}}{\tilde{P}}_{ph}^{C{O}_2}$ are the perturbed photoionization rate of O and CO₂, respectively. We note that $n_{\text{O}}^{\prime }/n_{\text{O}},n_{{\text{CO}}_2}^{\prime }/n_{{\text{CO}}_2},{\tilde{P}}_{ph}^O/P_{ph}^O,\text{and}{\tilde{P}}_{ph}^{C{O}_2}/P_{ph}^{C{O}_2}$ are assumed to be equal to ρ_n1/ρ_n0.

Appendix B Elements of the system of linear equation

The elements of the small four 5 × 5 matrices P, Q, R, and S in the absence of photochemical production and loss rates are displayed in this section. P is written as follows:

$$P=\left(\begin{array}{cc}{P}_{\text{i}}& 0\\ 0& P_{\text{e}}\end{array}\right),$$(B.1)

where P_i and P_e are 2 × 2 and 3 × 3 matrices and can be written as

$$P_{\text{i}}=\left(\begin{array}{cc}-i\omega & 0\\ p_{21}^i& p_{22}^i\end{array}\right),P_{\text{e}}=\left(\begin{array}{ccc}-i\omega & p_{12}^e& p_{13}^e\\ p_{21}^e& p_{22}^e& p_{23}^e\\ 0& p_{32}^e& p_{33}^e\end{array}\right).$$

The elements are

$$\begin{array}{l}{p}_{21}^i=-\left({\gamma }_{\text{i}}-1\right){\omega }_{\text{in}}-\left({\gamma }_{\text{i}}-1\right)\left[{\lambda }_{\text{i}}{k}^2-\left(k\cdot e_z\right)\partial {\lambda }_{\text{i}}/\partial z\right]T_{\text{i}0}/p_{\text{i}0},p_{22}^i=-i\omega +\left({\gamma }_{\text{i}}-1\right){\tilde{v}}_{\text{in}}+\left({\gamma }_{\text{i}}-1\right)\left[{\lambda }_{\text{i}}{k}^2-\left(k\cdot e_z\right)\partial {\lambda }_{\text{i}}/\partial z\right]T_{\text{i}0}/p_{\text{i}0}\\ p_{12}^e={\mu }_{0\text{e}}{p}_{\text{e}0}{\xi }_{\text{e}}{\zeta }_{\text{e}}/n_{\text{e}0},p_{13}^e={\mu }_{0\text{e}}e{k}_{\Vert }{\zeta }_{\text{e}},p_{21}^e=-\left({\gamma }_{\text{e}}-1\right){\omega }_{\text{en}}-k_{\Vert }\left({\gamma }_{\text{e}}-1\right)\left({\lambda }_{\text{e}}{k}_{\Vert }-i\partial {\lambda }_{\text{e}}/\partial b\right)T_{\text{e}0}/p_{\text{e}0},\\ p_{22}^e=-i\omega +\left({\gamma }_{\text{e}}-1\right){\tilde{v}}_{\text{en}}+{\mu }_{0\text{e}}{\xi }_{\text{e}}\left[\left({\gamma }_{\text{e}}-1\right)k_{\Vert }+{\xi }_{\text{e}}\right]p_{\text{e}0}/n_{\text{e}0}+k_{\Vert }\left({\gamma }_{\text{e}}-1\right)\left({\lambda }_{\text{e}}{k}_{\Vert }-i\partial {\lambda }_{\text{e}}/\partial b\right)T_{\text{e}0}/p_{\text{e}0}-ie{\mu }_{0\text{e}}{\xi }_{\text{e}}\left({\gamma }_{\text{e}}-1\right)\left(i{k}_{\Vert }{\beta }_{\text{e}}+\partial {\beta }_{\text{e}}/\partial b\right),\\ p_{23}^e=-{\mu }_{0\text{e}}e{k}_{\Vert }\left[\left({\gamma }_{\text{e}}-1\right)k_{\Vert }+{\xi }_{\text{e}}\right],p_{32}^e=-k_{\Vert }{\xi }_{\text{e}}{\sigma }_{0\text{e}},p_{33}^e=-i\omega {ϵ}_0{k}^2+{\sigma }_0{k}_{\Vert }^2+{\sigma }_{\text{p}}{k}_{\perp }^2,\end{array}$$

where ${\omega }_{\alpha n}={\tilde{v}}_{\alpha n}{\tilde{T}}_{\text{n}\alpha }$ denotes a characteristic ‘plasma-neutral energetic exchange frequency’; Here ${\tilde{T}}_{\text{n}\alpha }=T_{\text{n}0}/T_{\alpha 0}$ is the temperature ratio between neutral and plasma; k = (k_x, k_y, − i∂/∂ɀ) is the wavenumber vector in the local geographic coordinates O -xyz; k_|| = k · e_B is the field-aligned wavenumber; e_B is a unit vector of the magnetic field and written as e_B = (cos I sin D, cos I cos D, = sin I) in O -xyz coordinate system; Here I is the geomagnetic dip angle and D is the geomagnetic declination angle; k_⊥ = k = k_|| refers to the vertical component of wavenumber vector; ${\xi }_{\text{e}}=k_{\Vert }-i\sin I/H_{\text{e}}^{*}$, ζ_e = k_|| − i sin I/H_e, where H_e and $H_{\text{e}}^{*}$ are the electron density and pressure scale height, respectively; λ_i and λ_e are the thermal conductivity of ions and electrons; β_e denotes the thermoelectric coefficient; ∂/∂b = e_B · ∇ is a field-aligned differential operator.

The coupling matrix Q is given by

$$Q=\left(\begin{array}{ccccc}0& 0& i{k}_x& i{k}_y& \partial /\partial z-1/H_{\text{i}}\\ -\left({\gamma }_{\text{i}}-1\right){\omega }_{\text{in}}& 0& i{k}_x{\gamma }_{\text{i}}& i{k}_y{\gamma }_{\text{i}}& {\gamma }_{\text{i}}\partial /\partial z-1/H_{\text{i}}^{*}\\ 0& 0& i{\zeta }_{\text{e}}\cos I\sin D& i{\zeta }_{\text{e}}\cos I\cos D& -i{\zeta }_{\text{e}}\sin I\\ -\left({\gamma }_{\text{e}}-1\right){\omega }_{\text{en}}& 0& i{\widehat{\xi }}_{\text{e}}\cos I\sin D& i{\widehat{\xi }}_{\text{e}}\cos I\cos D& -i{\widehat{k}}_{\text{e}}\sin I\\ 0& 0& q_{53}& q_{54}& q_{55}\end{array}\right),$$(B.2)

where ${\widehat{\xi }}_{\text{e}}=\left({\gamma }_{\text{e}}-1\right)k_{\Vert }+{\xi }_{\text{e}};H_{\text{i}}\text{ and }{H}_{\text{i}}^{*}$ are the ion density and pressure scale height, respectively; q₅₃, q₅₄, and q₅₅ are

$$\begin{array}{cc}& q_{53}=i{B}_0\left[\left(k_{\text{H}}{\sigma }_{\text{p}}-k_{\text{p}}{\sigma }_{\text{H}}\right)\sin I\sin D-\left(k_{\text{p}}{\sigma }_{\text{p}}+k_{\text{H}}{\sigma }_{\text{H}}\right)\cos D\right],\\ & q_{54}=i{B}_0\left[\left(k_{\text{H}}{\sigma }_{\text{p}}-k_{\text{p}}{\sigma }_{\text{H}}\right)\sin I\cos D+\left(k_{\text{p}}{\sigma }_{\text{p}}+k_{\text{H}}{\sigma }_{\text{H}}\right)\sin D\right],\end{array}$$

and

$$q_{55}=i{B}_0\left(k_{\text{H}}{\sigma }_{\text{p}}-k_{\text{p}}{\sigma }_{\text{H}}\right)\cos I\text{,}$$

respectively. B₀ is the magnetic field magnitude in the absence of acoustic-gravity waves. k_p and k_H are the Pedersen and Hall components of the wavenumber vector, respectively.

The wavenumber vector in geographic and geomagnetic coordinates can be expressed as k = k_xe_x + k_ye_y − i(∂/∂z)e_ɀ and k = k_pe_p + k_He_H + k_||e_B, respectively. As illustrated in Fig. 2, the geomagnetic coordinates O − x′y′ɀ′, is obtained by clockwise rotation of the geographic coordinates by an angle of π/2 + I around y-axis and an angle of π/2 − D around z-axis. Therefore, the relation of the wavenumber components between two coordinates satisfies (k_p, k_H, k_||)^T = M_yM_ɀ(k_x, k_y, −i∂/∂ɀ)^T, where M_y = R_y(−π/2 − I) and M_ɀ = R_ɀ(−π/2 + D) with R_y and R_ɀ being the rotation matrix around y- and ɀ-axis, respectively. The total rotation matrix 𝔐 = M_yM_ɀ is given by the following

$$\begin{array}{l}\mathfrak{M}=M_y{M}_z=\left(\begin{array}{ccc}-\sin I& 0& -\cos I\\ 0& 1& 0\\ \cos I& 0& -\sin I\end{array}\right)\left(\begin{array}{ccc}\sin D& \cos D& 0\\ -\cos D& \sin D& 0\\ 0& 0& 1\end{array}\right)\\ =\left(\begin{array}{ccc}-\sin I\sin D& -\sin I\cos D& -\cos I\\ -\cos D& \sin D& 0\\ \cos I\sin D& \cos I\cos D& -\sin I\end{array}\right).\end{array}$$(B.3)

Similarly, the unit vector and differential operators in two coordinates also satisfy the relation mentioned above and can be expressed as (e_p, e_H, e_B)^T = 𝔐(e_x, e_y, e_ɀ)^T and (∂/∂x′, ∂/∂y′, ∂/∂b)^T = 𝔐(∂/∂x, ∂/∂y, ∂/∂ɀ)^T, respectively.

The coupling matrix R, consists of the neutral-plasma thermal exchange and photochemical heating effect, and can be given by the following

$$R=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ {\omega }_{\text{ni }}-{\omega }_{\text{pc}}& -{\omega }_{\text{ni }}& {\omega }_{\text{ne}}-{\omega }_{\text{pc}}& -{\omega }_{\text{ne}}& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$$(B.4)

where ${\omega }_{\text{n}\alpha }={\beta }_{\text{n}\alpha }{\tilde{T}}_{\alpha n}$ refers to a characteristic neutral-plasma energetic exchange frequency , where ${\beta }_{\text{n}\alpha }=k_{\text{B}}{\tilde{\nu }}_{\text{n}\alpha }/c_{\text{p}}{m}_{\text{n}},{\tilde{T}}_{\alpha n}=T_{\alpha 0}/T_{\text{n}0}$ is the temperature ratio between plasma and neutral; ω_pc = Q_pc/c_pρ_n0T_n0 denotes a ‘characteristic photochemical heating frequency’ with Q_pc being the photochemical heating of the hot oxygen atom. As displayed in Fig. B.1, ω_nα and ω_pc are less than the planetary angular frequency, Ω, and much less than the Brunt-Väisälä frequency, ω_b. Therefore, the effects of plasma behaviors on the acoustic-gravity waves can be reasonably ignored. The coupling matrix R may have a significant contribution to planetary-scale waves rather than the small-scale acoustic-gravity waves.

The differential operator matrix S is a linear combination of ∂²/∂ɀ² and ∂/∂ɀ and can be written as S = S₁∂²/∂ɀ² + S₂∂/∂ɀ + S₃, where S₁, S₂, and S₃ are 5 × 5 coefficient matrices. S₁ and S₂ are given by

$$S_1=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ -R{\kappa }_{\text{m}}/c_{\text{p}}& -{\kappa }_t& 0& 0& 0\\ 0& 0& -v_t& 0& 0\\ 0& 0& 0& -v_t& 0\\ 0& 0& 0& 0& -\left(v_t+v_{\text{m}}/3\right)\end{array}\right)$$(B.5)

and

$$S_2=\left(\begin{array}{ccccc}0& 0& 0& 0& 1\\ b_{21}& b_{22}& 0& 0& 0\\ 0& 0& {ϵ}_2& 0& -i{k}_x{v}_{\text{m}}/3\\ 0& 0& 0& {ϵ}_2& -i{k}_y{v}_{\text{m}}/3\\ p_{\text{n}0}/{\rho }_{\text{n}0}& 0& -i{k}_x{v}_{\text{m}}/3& -i{k}_y{v}_{\text{m}}/3& {\widetilde{ϵ}}_2\end{array}\right),$$(B.6)

respectively. κ_t = κ_m + κ_eddy and ν_t = ν_m + ν_eddy denote the total thermal diffusivity and total viscosity.

ϵ₂ = ϵ₂₀ + ϵ₂₁ and ${\widetilde{ϵ}}_2=4{ϵ}_{20}/3+{ϵ}_{21}$, where ${ϵ}_{20}=-{\rho }_{\text{n}0}^{-1}\partial {\mu }_{\text{m}}/\partial ɀ\text{ and }{ϵ}_{21}=-{\rho }_{\text{n}0}^{-1}\partial {\mu }_{\text{eddy}}/\partial ɀ$. Two elements b₂₁ and b₂₂ are

$$b_{21}=\frac{2R}{c_{\text{p}}}\frac{{\kappa }_{\text{m}}}{H^{\ast }}-\frac{R}{c_{\text{p}}}\frac{\partial {\kappa }_{\text{m}}}{\partial z}-\left({\kappa }_{\text{m}}+\frac{{\kappa }_{\text{eddy}}}{\gamma }\right)\frac{{\omega }_{\text{b}}^2}{g}$$

and

$$b_{22}=\left({\kappa }_{\text{m}}+\frac{{\kappa }_{\text{eddy}}}{\gamma }\right)\frac{1}{H^{\ast }}-\frac{\partial \left({\kappa }_{\text{m}}+{\kappa }_{\text{eddy}}\right)}{\partial z},$$

respectively. Here H* = (−∂lnp_n0/∂ɀ)¹ denotes the neutral pressure scale height.

S₃ can be divided into two parts and is given by ${\mathit{S}}_3={\tilde{\mathit{S}}}_3+\delta {\mathit{S}}_3$. The first part is

$${\tilde{S}}_3=\left(\begin{array}{ccccc}-i\omega /\gamma & i\omega & i{k}_x& i{k}_y& -1/H\\ c_{21}& -i\omega +{\kappa }_t{k}_{\text{h}}^2& 0& 0& {\omega }_{\text{b}}^2/g\\ i{k}_x{p}_{\text{n}0}/{\rho }_{\text{n}0}& 0& \tilde{d}+v_{\text{m}}{k}_x^2/3& v_{\text{m}}{k}_x{k}_y/3-f& i{k}_x{ϵ}_{20}\\ i{k}_y{p}_{\text{n}0}/{\rho }_{\text{n}0}& 0& v_{\text{m}}{k}_x{k}_y/3+f& \tilde{d}+v_{\text{m}}{k}_y^2/3& i{k}_y{ϵ}_{20}\\ -Rg/c_{\text{p}}& -g& -2i{k}_x{ϵ}_{20}/3& -2i{k}_y{ϵ}_{20}/3& \tilde{d}\end{array}\right),$$(B.7)

where $\tilde{d}=-i\omega +K_{\text{R}}+\left(v_{\text{m}}+v_{\text{eddy}}\right)k_{\text{h}}^2$. The element c₂₁ is

$$c_{21}=\frac{K_{\text{N}}}{\gamma }-\frac{{\omega }_{\text{b}}^2}{g}\frac{\partial \left({\kappa }_{\text{m}}+{\kappa }_{\text{eddy}}\right)}{\partial z}-\left({\kappa }_{\text{m}}+{\kappa }_{\text{eddy}}\right)\frac{\partial }{\partial z}\left(\frac{{\omega }_{\text{b}}^2}{g}\right)+\frac{R}{c_{\text{p}}}{\kappa }_{\text{m}}{k}_{\text{h}}^2+\frac{1}{H^{\ast }}\left[\left({\kappa }_{\text{m}}+\frac{{\kappa }_{\text{eddy}}}{\gamma }\right)\frac{{\omega }_{\text{b}}^2}{g}+\frac{R}{c_{\text{p}}}\frac{\partial {\kappa }_{\text{m}}}{\partial z}\right]-\frac{R}{c_{\text{p}}}\frac{{\kappa }_{\text{m}}}{H^{\ast 2}}\left(1+\frac{\partial H^{\ast }}{\partial z}\right)$$

The second part, consists of the neutral-plasma collisional energetic exchange effect and conductivity damping effect, can be given by the following

$$\delta S_3=\left(\begin{array}{cc}\delta S_{31}& 0\\ 0& \delta S_{32}\end{array}\right)$$(B.8)

δS₃₁ , which describes the energetic exchange effect, is a 2 × 2 matrix:

$$\delta S_{31}=\left(\begin{array}{cc}0& 0\\ {\tilde{c}}_{21}& {\tilde{c}}_{22}\end{array}\right)$$(B.9)

where elements ${\tilde{c}}_{21}=-\left({\omega }_{\text{ni}}+{\omega }_{\text{ne}}\right)/\gamma +\left({\beta }_{\text{ni}}+{\beta }_{\text{ne}}\right)\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}{\tilde{c}}_{22}=-\left({\omega }_{\text{ni}}+{\omega }_{\text{ne}}\right)$.

3 × 3 matrix δS₃₂ denotes the conductivity damping effect acting on the neutral waves and can be written as δS₃₂ = ω_pM_p + ω_HM_H, where ${\omega }_{\text{p}}={\sigma }_{\text{p}}{B}_0^2/{\rho }_{\text{n}0}\text{ and }{\omega }_{\text{H}}={\sigma }_{\text{H}}{B}_0^2/{\rho }_{\text{n}0}$ are characteristic frequency caused by Pedersen and Hall conductivity, respectively. M_p and M_H are symmetric and anti-symmetric matrix and given by the following

$$M\text{p}=\left(\begin{array}{ccc}1-{\cos }^{2}I{\sin }^{2}D& -{\cos }^{2}I\sin D\cos D& \sin I\cos I\sin D\\ -{\cos }^{2}I\sin D\cos D& 1-{\cos }^{2}I{\cos }^{2}D& \sin I\cos I\cos D\\ \sin I\cos I\sin D& \sin I\cos I\cos D& {\cos }^{2}I\end{array}\right)$$(B.10)

$$M_{\text{H}}=\left(\begin{array}{ccc}0& \sin I& \cos I\cos D\\ -\sin I& 0& -\cos I\sin D\\ -\cos I\cos D& \cos I\sin D& 0\end{array}\right)$$(B.11)

As shown in Fig. B.1, ω_ni, ω_ne, ω_p, and ω_H are much less than the Brunt-Väisälä frequency, ω_b, suggesting that the impacts of energetic exchange and conductivity damping on the acoustic-gravity waves can be reasonably ignored. We can assume that δS₃ ≈ 0 for small-scale neutral waves.

Appendix C Wentzel-Kramers-Brillouin solution

Acoustic-gravity waves are mainly dissipated by the molecular viscosity and thermal conduction in the upper atmosphere. The dissipation time-scale associated with molecular viscosity and thermal conduction can be defined by ω_υis = ν_m/H² and ω_hc = κ_m/H², respectively. Here H refers to the neutral density scale height. As displayed in Fig. B.1, ω_υis and ω_hc increase exponentially with height. In the lower region, ω_υis, ω_hc ≪ ω_b, suggesting the effect of the molecular viscosity and thermal conduction is insignificant and wave amplitude increases with altitude. ω_υis and ω_hc comparable to ω_b in the upper levels, resulting in wave attenuation and energy release. The wave amplitude begines decay at the damping altitude. The upper atmosphere can typically be treated as iso-thermal structure. The background parameters, including background temperature, molecular dynamic viscosity (µ_m), thermal conductivity (λ_m), slowly vary with height. In addition, the eddy and Coriolis effects are neglected. In this situation, the differential operator matrix S in the upper atmosphere can be simplified and could be partitioned into nine blocks:

$$S=\left(\begin{array}{lll}\widehat{A}\hfill & \widehat{B}\hfill & \widehat{C}\hfill \\ \widehat{D}\hfill & \widehat{E}\hfill & \widehat{F}\hfill \\ \widehat{G}\hfill & \widehat{H}\hfill & \widehat{J}\hfill \end{array}\right)$$(C.1)

The small matrices $\widehat{\mathit{A}}-\widehat{\mathit{J}}$ are given by

$$\begin{array}{l}\widehat{A}=\left(\begin{array}{cc}-i\omega /\gamma & i\omega \\ R{\kappa }_{\text{m}}\alpha /c_{\text{p}}& -i\omega +{\kappa }_{\text{m}}\alpha \end{array}\right),\widehat{B}=\left(\begin{array}{cc}i{k}_x& i{k}_y\\ 0& 0\end{array}\right),\widehat{C}=\left(\begin{array}{c}\partial /\partial z-1/H\\ {\omega }_{\text{b}}^2/g\end{array}\right),\\ \widehat{D}=Hg\left(\begin{array}{cc}i{k}_x& 0\\ i{k}_y& 0\end{array}\right),\widehat{E}=\left(\begin{array}{cc}\tilde{c}+v_{\text{m}}{k}_x^2/3& v_{\text{m}}{k}_x{k}_y/3\\ v_{\text{m}}{k}_x{k}_y/3& \tilde{c}+v_{\text{m}}{k}_y^2/3\end{array}\right),\widehat{F}=-\frac{v_{\text{m}}}{3}(\begin{array}{c}i{k}_x\\ i{k}_y\end{array})\frac{\partial }{\partial z},\\ \widehat{G}=\left(\begin{array}{cc}Hg\frac{\partial }{\partial z}-\frac{Rg}{c_{\text{p}}},& -g\end{array}\right),\widehat{H}=-\frac{v_{\text{m}}}{3}\left(\begin{array}{cc}i{k}_x,& i{k}_y\end{array}\right)\frac{\partial }{\partial z},\widehat{J}=((4\tilde{c}-\tilde{d})/3),\end{array}$$

where $\alpha =-{\partial }^2/\partial z^2+k_{\text{h}}^2,\tilde{c}=-i\omega +v_{\text{m}}\alpha ,\tilde{d}=-i\omega +v_{\text{m}}{k}_{\text{h}}^2$. We note that $H=H^{*}=R{T}_{\text{n}0}/g,{\omega }_{\text{b}}^2=Rg/c_{\text{p}}H,\text{ and }{p}_{\text{n}0}/{\rho }_{\text{n}0}=Hg=c_{\text{s}}^2/\gamma$ in the iso-thermal atmosphere.

Matrix S can be decomposed into two factors, a lower triangular matrix L and an upper triangular matrix U and can be written as S = LU, where

$$L=\left(\begin{array}{ccc}I& 0& 0\\ L_{21}& I& 0\\ L_{31}& L_{32}& I\end{array}\right),U=\left(\begin{array}{ccc}\widehat{A}& \widehat{B}& \widehat{C}\\ 0& U_{22}& U_{23}\\ 0& 0& U_{33}\end{array}\right)$$(C.2)

Here I denotes the unit matrix. Other element matrices are

$$\begin{array}{l}{L}_{21}=\widehat{D}{\widehat{A}}^{-1},L_{31}=\widehat{G}{\widehat{A}}^{-1},U_{22}=\widehat{E}-L_{21}\widehat{B},U_{23}=\widehat{F}-L_{21}\widehat{C},\\ L_{32}=\left(\widehat{H}-L_{31}\right)U_{22}^{-1},U_{33}=\widehat{J}-L_{31}\widehat{C}-L_{32}{U}_{23}\end{array}$$

Since the determinant of matrix L equals 1, the second-order differential Eq. (44a) is equivalent to UX_neu = 0. Here X_neu = x₁, x₂, w₁)^T, where x₁ = (p_n1/p_n0, s₁/c_p)^T , x₂ = (u₁,υ₁)^T . To find a ‘non-trivial’ solution of this linear equation, we assume that the coefficient matrix determinant equals zero, $\left|\mathit{U}\right|=\left|\widehat{\mathit{A}}\right|\left|{\mathit{U}}_{22}\right|\left|{\mathit{U}}_{33}\right|=0$, where

$$\begin{array}{ll}\hfill & \left|\widehat{A}\right|=\frac{\omega }{\gamma }\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right),\hfill \\ \hfill & \left|U_{22}\right|=\tilde{c}\left[\left|\widehat{A}\right|\tilde{c}+k_{\text{h}}^2\left(\frac{v_{\text{m}}}{3}\left|\widehat{A}\right|-iHg\left(\omega +i{\kappa }_{\text{m}}\alpha \right)\right)\right],\hfill \\ \hfill & \left|U_{33}\right|=\frac{\tilde{c}\mathfrak{L}}{\left|U_{22}\right|}\hfill \end{array}$$

When |Â| = 0 and |U₂₂| = 0, the solutions describe extraordinary heat conduction waves and extraordinary viscosity waves (Volland 1969b; Hickey & Cole 1987). In this study, the acoustic-gravity wave solution is only considered and can be obtained by only assuming |U₃₃| = 0. We first solve w₁, which satisfies U₃₃w₁ = 0. This equation is equivalent to 𝔏w₁ = 0. 𝔏 is a sixth-order differential operator:

$$\mathfrak{L}=\left(\omega +i{\kappa }_{\text{m}}\alpha \right)\left(\omega +i{v}_{\text{m}}\alpha \right)\left(\frac{{\partial }^2}{\partial z^2}-\frac{1}{\text{H}}\frac{\partial }{\partial z}-k_{\text{h}}^2\right)-Q_1$$(C.3)

where

$$Q_1=k_{\text{h}}^2{\omega }_{\text{b}}^2+\frac{\omega }{c_{\text{s}}^2}\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)\left(\omega +i{v}_{\text{m}}\alpha \right)\left[\omega +i{v}_{\text{m}}\alpha \left(1+\frac{1}{3}\right)\right]$$

We can reasonably assume that the wave solution in the iso-thermal atmosphere satisfies X_neu(ɀ) ∝ exp (∫iK_zdɀ), where Kɀ = k_ɀ − i/2H, k_ɀ denotes the vertical wavenumber. In this situation, the differential operator ∂/∂_ɀ is equivalent to iK_ɀ. Therefore, the differential operator 𝔏 is a sixth-order polynomial of k_ɀ. The wave dispersion relation can be obtained by assuming 𝔏(k_ɀ) = 0 and is given by Eq. (48). In addition, the differential equation 𝔏w₁ = 0 can be rewritten as

$$\frac{{\partial }^2{w}_1}{\partial z^2}-\frac{1}{\text{H}}\frac{\partial w_1}{\partial z}-\left[k_{\text{h}}^2+Q(z)\right]w_1=0,$$(C.4)

where Q(ɀ) = Q₁/ (ω + iκ_mα) (ω + iν_mα), here $\alpha =k_{\text{h}}^2-{\left(i{k}_z+1/2H\right)}^2$. By assuming ${\widehat{w}}_1=w_1\exp \left(-{\displaystyle \int 1}/2Hdɀ\right)$, Eq. (C.4) becomes Eq. (46). The Wentzel-Kramers-Brillouin (WKB) solution of Eq. (C.4) is given by Eq. (49) in Sect. 2.3.1. We note that the WKB wave solution (49) satisfies

$$\frac{{\partial }^2{\widehat{w}}_1}{\partial z^2}+k_z(z)[1+d(z)]{\widehat{w}}_1=0,$$(C.5)

where d(ɀ) is the residue of the WKB solution and can be given by the following (Hines & Reddy 1974)

$$d(z)=\frac{1}{2k_z^3}\frac{{\partial }^2{k}_z^2}{\partial z^2}-\frac{3}{4k_z^4}{\left(\frac{\partial k_z}{\partial z}\right)}^2$$(C.6)

Eq. (C.5) approaches Eq. (46) if d(ɀ) ≪ 1. Therefore, the WKB approximation method is valid in the upper atmosphere, in which the vertical wavenumber, k_ɀ, can be slowly vary with height.

Other solutions, including u₁, υ₁, p_n0, and s₁, are related to w₁ and can be given by the wave polarization relations. The first two columns of equation UX_neu = 0 yield the following relation

$$\left(\begin{array}{cc}\widehat{A}& \widehat{B}\\ 0& U_{22}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=-\left(\begin{array}{c}\widehat{C}\\ U_{23}\end{array}\right)w_1$$(C.7)

x₁ and x₂ can be obtained by solving the above equation and given by

$$\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=\left(\begin{array}{c}{\widehat{A}}^{-1}\left(\widehat{B}{U}_{22}^{-1}{U}_{23}-\widehat{C}\right)\\ -U_{22}^{-1}{U}_{23}\end{array}\right)w_1$$(C.8)

Therefore, the polarization relation is given by

$$\frac{p_{\text{n}1}}{p_{\text{n}0}}=\frac{i}{\omega {\beta }_{0}H}\left[\gamma {\beta }_1\left(\omega +i{v}_{\text{m}}\alpha \right)+\frac{i{v}_{\text{m}}{k}_{\text{h}}^2}{3}\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)\right]w_1,$$(C.9a)

$$\frac{s_1}{c_{\text{p}}}=-\frac{R}{c_{\text{p}}}\frac{i}{\omega {\beta }_{0}H}\left\{\left(\omega +i{v}_{\text{m}}\alpha \right)\left[\omega -i\gamma {\kappa }_{\text{m}}\alpha \left(i{k}_{z}H-\frac{1}{2}\right)\right]-k_{\text{h}}^2{c}_{\text{s}}^2+\frac{i{v}_{\text{m}}{k}_{\text{h}}^2}{3}\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)\right\}{w}_1,$$(C.9b)

$$x_2=\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)=\frac{i}{\omega {\beta }_{0}H}\left[\frac{i{v}_{\text{m}}}{3}\omega \left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)\left(i{k}_{z}H+\frac{1}{2}\right)+{\beta }_1{c}_{\text{s}}^2\right]\left(\begin{array}{c}{k}_x\\ k_y\end{array}\right)w_1\text{. }$$(C.9c)

The perturbed neutral density and temperature are given by

$$\frac{{\rho }_{\text{n}1}}{{\rho }_{\text{n}0}}=\frac{1}{\gamma }\frac{p_{\text{n}1}}{p_{\text{n}0}}-\frac{s_1}{c_{\text{p}}}=\frac{R}{c_{\text{p}}}\frac{i}{\omega {\beta }_{0}H}\left[\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right){\beta }_2\frac{c_{\text{p}}}{R}-k_{\text{h}}^2{c}_{\text{s}}^2\right]w_1,$$(C.9d)

$$\frac{T_{\text{n}1}}{T_{\text{n}0}}=\frac{R}{c_{\text{p}}}\frac{p_{\text{n}1}}{p_{\text{n}0}}+\frac{s_1}{c_{\text{p}}}=-\frac{R}{c_{\text{p}}}\frac{i}{\omega {\beta }_{0}H}\left[\gamma \omega \left(\omega +i{v}_{\text{m}}\alpha \right)\left(i{k}_{z}H+\frac{1}{2}\right)-k_{\text{h}}^2{c}_{\text{s}}^2\right]w_1$$(C.9e)

where

$$\begin{array}{ll}\hfill & {\beta }_0=\left(\omega +i\gamma {\kappa }_{\text{m}}\alpha \right)\left[\omega +i{v}_{\text{m}}\left(\alpha +\frac{k_{\text{h}}^2}{3}\right)\right]-\frac{k_{\text{h}}^2{c}_{\text{s}}^2}{\omega }\left(\omega +i{\kappa }_{\text{m}}\alpha \right),\hfill \\ \hfill & {\beta }_1=-\left[\left(\omega +i{\kappa }_{\text{m}}\alpha \right)\left(i{k}_{z}H-\frac{1}{2}\right)+\frac{R}{c_{\text{p}}}\omega \right],\hfill \\ \hfill & {\beta }_2=-\left(\omega +i{v}_{\text{m}}\alpha \right)\left(i{k}_{z}H-\frac{1}{2}\right)+\frac{i{v}_{\text{m}}}{3}{k}_{\text{h}}^2\hfill \end{array}$$

In the lower atmosphere, the effect of the molecular viscosity and thermal conduction is insignificant. The polarization relations in the lower region, given by Eqs. (54a) - (54d), can be obtained by ignoring ν_m and κ_m in Eqs. (C.9a) - (C.9c).

Appendix D Numerical method of second-order ordinary differential equation with variable coefficients

Eq. (45) can be discretized by using a finite-difference method in the region between the lower boundary and upper boundary, Altitude can be divided into N+1 points: ɀ₀, ɀ₁, ɀ₂,…, ɀ_N with equal intervals (δɀ = ɀ_i+1 − ɀ_i). Hence, differential operators can be written as in form of finite-difference:

$$\frac{d^{2}X}{d{z}^2}=\frac{X_{\text{n}+1}-2X_{\text{n}}+X_{\text{n}-1}}{{(\delta z)}^2},\frac{dX}{dz}=\frac{X_{\text{n}+1}-X_{\text{n}-1}}{2\delta z}$$(D.1)

where n = 1, 2,…, N − 1.

A second-order ordinary differential equation with variable coefficients is expressed as:

$$A(z)\frac{d^{2}X}{d{z}^2}+B(z)\frac{dX}{dz}+C(z)X=r(z)$$(D.2)

Eq. (D.2) can be discretized by:

$${\tilde{A}}_{\text{n}}{X}_{\text{n}-1}+{\tilde{B}}_{\text{n}}{X}_{\text{n}}+{\tilde{C}}_{\text{n}}{X}_{\text{n}+1}=r\left(z_{\text{n}}\right),n=1,2,\dots ,N-1$$(D.3)

where

$${\tilde{A}}_{\text{n}}=\frac{A\left(z_{\text{n}}\right)}{{(\delta z)}^2}-\frac{B\left(z_{\text{n}}\right)}{2\delta z},{\tilde{B}}_{\text{n}}=C\left(z_{\text{n}}\right)-\frac{2A\left(z_{\text{n}}\right)}{{(\delta z)}^2},{\tilde{C}}_{\text{n}}=\frac{A\left(z_{\text{n}}\right)}{{(\delta z)}^2}+\frac{B\left(z_{\text{n}}\right)}{2\delta z}$$(D.4)

Dirichlet boudnary condition is that the values at the lower (X₀) and upper (X_N) boundaries are known. the values between the lower and upper boundary, X_n(n = 1, 2, 3,…, N − 1), can be solved by the linear equation given by Eq. (D.5). In our model, the height

$$\left[\begin{array}{cccccccccc}{\tilde{B}}_1& {\tilde{C}}_1& 0& \dots & \dots & 0& 0& 0& 0& 0\\ {\tilde{A}}_2& {\tilde{B}}_2& {\tilde{C}}_2& 0& \dots & 0& 0& 0& 0& 0\\ 0& {\tilde{A}}_3& {\tilde{B}}_3& {\tilde{C}}_3& \dots & 0& 0& 0& 0& 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& 0& \dots & 0& {\tilde{A}}_{\text{N}-3}& {\tilde{B}}_{\text{N}-3}& {\tilde{C}}_{\text{N}-3}& 0\\ 0& 0& 0& 0& \dots & 0& 0& {\tilde{A}}_{N-2}& {\tilde{B}}_{N-2}& {\tilde{C}}_{N-2}\\ 0& 0& 0& 0& \dots & 0& 0& 0& {\tilde{A}}_{N-1}& {\tilde{B}}_{\text{N}-1}\end{array}\right]\left[\begin{array}{c}{X}_1\\ X_2\\ X_3\\ \cdots \\ \cdots \\ X_{N-3}\\ X_{N-2}\\ X_{\text{N}-1}\end{array}\right]=\left[\begin{array}{c}r\left(z_1\right)-{\tilde{A}}_1{X}_0\\ r\left(z_2\right)\\ r\left(z_3\right)\\ \cdots \\ r\left(z_{\text{N}-3}\right)\\ r\left(z_{\text{N}-2}\right)\\ r\left(z_{\text{N}-1}\right)-{\tilde{C}}_{\text{N}-1}{X}_{\text{N}}\end{array}\right]$$(D.5)

resolution is δɀ =1 km. The altitude ranges from 80 km to 300 km. X₀ and X_N are the non-dissipative wave solution and WKB wave solution, respectively.

References

All Tables

All Figures