Buoyancy-driven Magnetohydrodynamic Waves in a Partially Ionized Plasma

A magnetohydrodynamic (MHD) fluid description is typically employed to study the magnetized plasma comprising the solar atmosphere. This approach has had many successes in modeling and explaining solar phenomena. Most often, the plasma is assumed to be fully ionized. While this approach is justified in the higher atmosphere, i.e., the solar corona; the temperature in the lower solar atmosphere is such that a large proportion of the fluid may be electrically neutral. This begs the question: to what degree are the results derived from a fully ionized MHD description valid? In this article, we investigate the effect of partial ionization on buoyancy-driven MHD waves (the MHD analog of internal gravity waves) by applying a simplified two-fluid description. We show that previously derived results may be applied, when the fluid is weakly ionized, if the ion–neutral collision frequency is high. We derive dispersion relations for buoyancy-driven MHD waves, which include correction factors and damping rates due to ion–neutral collisions.


Introduction
Magnetohydrodynamic effects are observed in the lower solar atmosphere. Sunspots, pores, bright points, etc., show that magnetic fields do indeed influence the dynamics of the solar plasma in the photosphere. When studying perturbations in magnetic structures, it is often assumed that the particles are fully ionized and so a magnetohydrodynamic (MHD) description of the fluid may be used. To include the effects of the magnetic field on perturbations, the most simple MHD model, a single fluid composed of fully ionized material, is frequently employed. In reality, the solar plasma is not composed of fully ionized hydrogen particles, rather the fluid also consists of neutral hydrogen particles (and some heavier elements that may display varying levels of ionization). To include the more realistic effects of partial ionization, we may have to use a modified description of MHD.
A three-fluid description of a plasma, where distinctions are made between electrons, ions, and neutral hydrogen is discussed, based on first principles, by, e.g., Braginskii (1965). MHD waves in a two-fluid MHD model (where the electrons and ions may be considered as a single fluid, valid for timescales longer than the ion-electron collision time) were studied by, e.g., Zaqarashvili et al. (2011b) and Soler et al. (2013aSoler et al. ( , 2013bSoler et al. ( , 2013c. Note that these alternate MHD models do not include the effects of elemental abundances, though these effects are included in many widely used numerical codes. Alfvén waves in plasma containing hydrogen and helium were studied by Zaqarashvili et al. (2011a). The previous works did not include the effects of buoyancy or gravitational stratification.
Internal gravity waves (IGWs), where buoyancy plays the role of the dominant resorting force, have been observed to propagate through the lower solar atmosphere by, e.g., Komm et al. (1991), Stodilka (2008), and Straus et al. (2008). The theory of linear buoyancy-driven MHD waves (that is, the MHD analog of IGWs) has been developed by, e.g., Hasan & Christensen-Dalsgaard (1992), Barnes et al. (1998), Roberts (2006, and Hague & Erdélyi (2016). These previous works focused on buoyancy-driven MHD waves in a single, fully ionized fluid. A numerical study of IGWs applicable to the solar atmosphere was recently undertaken by Vigeesh et al. (2017).
IGWs modified by magnetic fields in the lower solar atmosphere (that is, slow magnetoacoustic gravity (MAG) waves in the high-β plasma) are likely excited by turbulent motion close to the solar surface. This mechanism was previously suggested by Komm et al. (1991). Numerical simulations by Brun et al. (2013) have shown that g-modes in the radiative interior can be generated by turbulent convection. IGWs generated close to the visible surface of the Sun may propagate higher into the atmosphere. The degree of ionization of the solar plasma changes through the solar atmosphere. The lower atmosphere is such that only a small amount of matter is ionized. Hence, a realistic treatment of atmospheric IGWs requires that we take into account the effects of partial ionization.
In this article, we use a two-fluid MHD model to analytically study the effects of partial ionization on buoyancy-driven MHD waves. Here, we would like to investigate the physical domain where partial ionization and buoyancy both play roles in the wave dynamics. The character of the resulting wave solutions is highly dependent on the timescale of the motion, relative to the ion-neutral collision timescale. We investigate the possible wave solutions analytically and numerically. The inclusion of neutral-ion collisions introduces a damping mechanism for the wave motion.
We begin, in Section 2, by deriving the governing equations for buoyancy-driven MHD waves in a simple two-fluid model. We derive a local dispersion relation for propagating waves using the Wentzel-Kramers-Brillouin (WKB) approximation (see, e.g., Bender & Orszag 1978). The dispersion relation is solved numerically and the solutions are discussed. Asymptotic solutions to the dispersion relation are presented in Section 3. The solutions represent generalizations of previously derived results for buoyancy-driven MHD waves. The solutions also allow us to investigate the conditions under which a fully ionized fluid description may be used to approximate the weakly ionized lower solar atmosphere. Damping rates, due to collisions between ions and neutral particles, are derived in various limits. We conclude in Section 4.

Derivation of the Dispersion Relation
To study the effects of partial ionization, we use a two-fluid model. This approach treats the plasma as a mixture of an ionelectron gas and a neutral hydrogen gas. The ion-electron gas is taken to be quasi-neutral (n i ≈n e ). We are interested in the effect of partial ionization on buoyancy-driven MHD waves, hence we apply a Boussinesq-type approximation to the ionelectron and neutral gases. The approximation applied essentially retains density perturbations only where a buoyancy force is produced. This is sometimes referred to as the anelastic approximation; see, e.g., Lighthill (1978) for a discussion on the approximation (note that the term "Boussinesq approximation" is used in this work). The Boussinesq-type approximation is applicable to high-β plasmas, as discussed by Hague & Erdélyi (2016).
We assume a time-independent, static background stratified under gravity. Gravity is taken to be constant and in the negative z-direction, g=(0, 0, −g), the background magnetic field is assumed constant and parallel to gravity B 0 =(0, 0, B 0 ). The background pressures and densities are, therefore, functions of z-only.
The system of governing equations used in this article are those of Díaz et al. (2012) with the addition of the Boussinesq approximation. For simplicity, we neglect electron inertia, collisions between electrons and neutrals, the effect of gravity on electrons, and nonisotropic components of the pressure tensor as in Zaqarashvili et al. (2011b) and Díaz et al. (2012). Magnetic diffusion is neglected in the induction equation, hence the only terms present are ideal terms (a more complete treatment of the induction equation for partially ionized, stratified plasmas is given by Díaz et al. (2014); the additional effects are not considered here).
The linearized equations for perturbations around this background are where v, p, and ρ represent the fluid velocity, pressure, and density, respectively. The subscript 0 refers to the background quantity, while 1 refers to the perturbed quantity. N is the Brunt-Väisälä frequency, defined by where c s is the sound speed given by c p , and γ is the ratio of specific heats.
The subscripts n and i in the governing equations refer to the neutral and ionized fluids, respectively, and α in is the equilibrium coefficient of friction between neutrals and ions. We assume no y-dependence, so that we may separate Alfvén and slow MAG waves (see Hague & Erdélyi 2016 for a discussion on the equivalence of slow waves and magnetic IGWs). The Alfvén waves, in the current geometry, are represented by the y-component of the velocity field, which we do not consider. First, we Fourier analyze in x and t; that is, we assume that perturbed quantities are proportional to where k x and ω are the horizontal wavenumber and frequency, respectively. Some algebra leads to the coupled governing equations for the vertical velocity components, The above equations are the governing equations, which we will use to study buoyancy-driven MHD waves in a simplified two-fluid model. We study the wave solutions via the local dispersion relation, which may be found using the WKB method (see, e.g., Bender & Orszag 1978). To apply the WKB method, a slowly varying spatial variable z z, Substituting Equation (9) into Equations (7) and (8) and keeping leading-order terms (that is, terms order 0  ( )) leads, after some manipulation to the local dispersion relation, = + The fourth-order polynomial can be solved exactly; though, the expressions themselves do not reveal the nature of the solutions. We note that there are four solutions, hence there are two branches of solutions. Zaqarashvili et al. (2011b) also found four slow solutions. In the high-frequency domain (compared to the collision frequency), the two gases do not couple strongly, so each may support oscillatory wave modes. In the low-frequency regime these modes have zero real part. In this case, the solutions were termed vortex modes; these are solutions to the fluid equations that damp due to collisions. The solutions do not represent oscillatory waves as the ions and neutrals are coupled strongly through collisions. Therefore, the medium behaves as a single fluid, hence there are only two slow solutions. We will investigate these solutions in more detail.
Let us now consider the dispersion relation (10) by plotting the solutions numerically. We assume that the temperature of all three species are equal, that is T i =T e =T n ; hence, c c 2 .
si sn The factor of 2 is due to the fact that, for hydrogen plasma, the mean particle mass of the neutral fluid is twice that of the ionized fluid. This result may also be derived from the ideal gas equation. For simplicity, we consider the gases to be isothermal so the square of the Brunt-Väisälä frequency is given by N g c 1 .
i n s i n , We use these values in the following numerical solutions.
Let us begin by investigating the effect of varying the collision frequency. Figures 1 and 2 show the real and imaginary parts of the frequency when ν in /N varies from 0 to 1.5. We plot the case of 70% neutral gas and 30% ionized gas. These values are chosen primarily for illustrative purposes, where the interesting physical regimes of the solutions can be clearly seen. The values are representative of conditions that may be found in the upper chromosphere, at the base of the transition region (Carlsson & Stein 2002). It should be noted that the assumptions used in this analysis may be used in this region only for a sufficiently weak magnetic field (the Boussinesq approximation may fail for a stronger magnetic field in the chromosphere). The other parameters take the values k z /k x =1, v A k z /N=1. It can be shown that the curves plotted follow the same paths for a large range of ν in /N. The behavior plotted here is, therefore, representative of arbitrarily large ν in . The solid curves, in Figure 1, represent forward and backward propagating slow MAG waves that remain oscillatory for all values of ν in . In the limit of small collision frequency, the frequency of these waves may be determined via asymptotic analysis (Equation (19); this expression is close to the numerical solution for ν in /N<0.2). Figure 1 clearly shows two regimes for the solution. For increasing collision frequency the wave frequency decreases. When ν in /N≈0.55 the real part of the solutions cease to depend on the collision frequency. The frequency is given by Equation (24), valid for "high" collision frequencies. Beyond this value, the plasma acts as a single fluid, a higher collision frequency does not change this fact so the real part of the propagating wave is unaffected by the collision frequency.
There is another solution present that is represented by the dotted line. This solution is a propagating mode in the low collision frequency regime and a purely damped vortex mode when the partially ionized plasma acts as a single fluid. The value of the real part of the frequency when the wave is propagating is given by the second term of Equation (19) for ν in /N<0.4. Beyond this value the real part of the frequency tends to 0. When the partially ionized particles interact weakly, Figure 1. Real part of the frequency for varying collision frequency, normalized to N. The other parameters take the values ξ n =0.7, ξ i =0.3, k z /k x =1, v A k z /N=1. The solid curves represent propagating slow MAG waves in either the ionized fluid (in the two-fluid limit) or the coupled fluid (in the single-fluid limit). The dotted curve represents IGWs in the neutral fluid, which propagate when the two fluids do not couple strongly, that become purely damped vortex modes in the single-fluid limit. Figure 2. Imaginary part of the frequency (damping rate) for varying collision frequency, normalized to N. The other parameters take the values ξ n =0.7, ξ i =0.3, k z /k x =1, v A k z /N=1. The solid and dotted curves represent the damping rate of wave modes described in Figure 1. the dotted curves correspond to IGWs in the neutral fluid, while the solid curves are buoyancy-driven MHD waves in the ionized fluid. When the plasma species interact strongly as a single fluid, the solid curves represent buoyancy-driven MHD waves in the whole fluid. Figure 2 shows the imaginary part of the frequency. We see that all curves take negative values, indicating that the solutions are damped. We, again, see the two regimes where the fluid may be treated as separate fluids or a single fluid. In the low ν in region, the MAG wave (solid curve) is damped more strongly than the IGW (dotted curve). The damping rates for small ν in are given by Equation (20). As expected, the damping rates for the forward and backward waves are equal by symmetry. In the high ν in region, the damping rate of the propagating modes decrease with ν in , tending asymptotically to 0. This is expressed in solution (24), which contains no imaginary part. The physical reason for this is that the wave acts on a greater timescale than the ion-neutral collision timescale and so does not feel the particle interactions in this limit. When the propagating IGWs become imaginary vortex modes, the damping rate bifurcates; there are two vortex solutions that are strongly and weakly damped, respectively.
In the solar atmosphere, the Brunt-Väisälä frequency is typically low compared to the collision frequency. Buoyancy, therefore, plays a key role in the motion of low-frequency waves. We see from the above that single-fluid effects are dominant for low-frequency waves with an IGW character. Two-fluid effects may become important for higher-frequency waves, i.e., waves with a magnetic character. Figures 3 and 4 show the effect of varying the magnetic field strength. Again, the solid curves show the propagating solutions and dotted curves show the vortex/propagating modes. As above, we set the parameters to the values ξ n =0.7, ξ i =0.3, k z /k x =1. The collision frequency is set so that ν in /N=0.5, close to the region where propagating modes become vortex modes. It can be shown that the trends seen in Figures 3 and 4 persist for larger Alfvén frequencies. A key point depicted here is that the behavior of the solutions can be split into magnetically or gravitationally dominated regions. The propagating MAG waves undergo a change in the gradient of the real part of the frequency around v A k z /N=1. The gradient of the imaginary part also undergoes a similar change.
It is interesting to note that the parameters of magnetic field strength and collision frequency produce contrasting effects. That is, for a weak magnetic field, the vortex modes are purely imaginary; increasing the magnetic field causes them to become propagating modes. There is, again, a weakly and strongly damped vortex mode.
Note that for large Alfvén frequency the solutions here do not match up with the analytical results for large vertical wavenumber. This is due to the fixing of the ratio k z /k x =1. It can be shown that plotting the solutions for large Alfvén frequency and large values of k z /k x does indeed agree with the asymptotic results.

Limiting Cases of the Dispersion Relation
Here, we note some limiting forms of the dispersion relation. In the nonmagnetic, collisionless (ν in =0) limit, the dispersion relation becomes k k N k k N 0, 11 x i x n 2 2 2 2 2 2 2 2 which represent IGWs in each of the noninteracting fluids. In the collisionless, homogeneous (N i,n =0) limit the dispersion relation becomes v k 0 , which represents the familiar slow MHD wave. Let us now consider the more interesting cases where the two fluids may interact. First, we consider the plasma to be strongly ionized, ξ n =1. In this limit, the dispersion relation reduces to This represents the buoyancy-driven slow mode with a correction due to the presence of a small amount of neutral fluid. The form of the dispersion equation is similar to the case of the fully ionized plasma, studied previously in, e.g., Hague & Erdélyi (2016). In this limit, the fourth-order equation has reduced to a second-order one, hence some solutions have been removed. These are the modes associated with the neutral gas. The other parameters take the values ξ n =0.7, ξ i =0.3, k z /k x =1, ν in /N=0.5. As in Figure 1, the solid curves represent slow MAG waves in the ionized fluid (two-fluid limit) or the coupled fluid (single-fluid limit). The mode may take on a gravity or magnetic character depending on the magnetic field strength. The dotted curve represents a purely damped vortex mode (single-fluid limit) or an IGW in the neutral fluid (twofluid limit). Perhaps of more interest when considering the lower solar atmosphere, is the limit of a weakly ionized plasma. As discussed in Hague & Erdélyi (2016), we expect the Boussinesq approximation to be more applicable in the lower solar atmosphere, where a significant proportion of the fluid may be neutral. At the base of the photosphere, the level of ionization is placed at around 1% (Carlsson & Stein 2002). The ionization fraction decreases significantly, and rapidly, with height through the photosphere. At the temperature minimum, approximately 500 km above the base of the photosphere, the ionization fraction takes a value of 10 −5 (i.e., 10 −3 % of the fluid is ionized). The ionization fraction steadily rises through the chromosphere, to a value of around 30% ionization. Let us consider the case of ξ i =1, where the dispersion relation takes the form This is the dispersion relation with which we may study photospheric plasma. Again, we have lost a solution associated with the ionized gas. While this equation is less complex than the full dispersion relation, the solutions are not in a simple form. In the interest of simplicity, we will investigate the solutions by making a further approximation, that of small/ large collision frequency relative to the wave frequency. These limits represent the extent to which fluid species interact. For processes that act on timescales shorter than the collision time, the waves do not strongly feel the interaction between the two fluids. On long timescales, the waves are highly influenced by the interaction of the fluids. We begin by assuming that ν in is small (when nondimensionalized appropriately by, e.g., the Alfvén or Brunt-Väisälä frequency). The dispersion relation is singular with regards to ν in , hence the second solution associated with the ionized gas is removed. We use a perturbation series ω=ω 0 +ν in ω 1 /(v A k z )+ .... The leadingorder solution may be found to be k k N , 1 4 x n 0 2 2 2 2 w = ( ) which represents IGWs. The leading-order effect of small collision frequency is that magnetic effects are unimportant. The next term may be found to be In the above analysis, we made the assumption that ξ i =1, which implies that ξ n ≈1. The first term in ω 1 may then be considered to be small, and so contributes to the third term in the perturbation series. The frequency then takes the form We see that for the case of a small collision frequency, the waves may be considered damped IGWs in the mostly neutral plasma. The damping is due to ion-neutral collision and the damping rate is given by the imaginary part of the frequency. The above analysis is valid when the g-mode frequency is larger than the collision frequency. This is not the case for the solar atmosphere, although it may have application in other stars, such as Red Giants, where the the Brunt-Väisälä frequency may become higher. See below for results more applicable to higher-frequency waves in the solar atmosphere.
In the solar atmosphere, the Brunt-Väisälä frequency is relatively small, hence we should also consider the case of a large collision frequency. Zaqarashvili et al. (2011b) estimate the ion-neutral collision frequency in the chromosphere, based on the FAL-3 model (Fontenla et al. 1990), to be 4 Hz. Hence, the case of high collision frequency may be more applicable to the lower solar atmosphere. This limit is of interest for lower frequency, buoyancy-driven MHD waves. That is, slow MAG waves (which may take on a buoyancy or magnetic character) where the wave frequency is smaller than the collision frequency. Let the frequency take the form We see, when the collision frequency is large, the frequency is of the form of a fully ionized plasma and magnetic effects may be important even when the fluid is weakly ionized. There is also a solution ω 0 =0, representing a damped vortex mode (see below). The next term in the series is given by which is imaginary, and so the damping rate is given by v A k z ω 1 /ν in . In the above analysis, we have shown that the frequency, derived for the case of fully ionized fluid, may be used in partially ionized plasmas when either the plasma is strongly ionized or weakly ionized and the ion-neutral collision frequency is high. We have calculated generalizations to these results for either case.
The nature of the solutions may also be investigated by applying the limit of small collision frequency to the full fourth-order dispersion relation, Equation (10). In this limit, the dispersion relation is not singular, hence we have all four solutions. These results are valid for arbitrary proportion of ionized and neutral material. As above, let us use a perturbation series of the form ω=ω 0 +ν in ω 1 /(v A k z )+..., where ν in = v A k z . To lowest order, we find the leading-order solutions to be These solutions represent buoyancy-driven waves in each of the fluids. As the collision frequency is small, the coupling of the fluids is weak. As previously noted, the limit of small collision frequency is applicable to high-frequency waves (relative to the collision frequency). In the solar atmosphere, this limit is not appropriate for low-frequency, IGW-type waves. Slow MAG waves take on a magnetic character in the high-frequency domain. Hence, in the solar atmosphere, twofluid effects are dominant for higher-frequency waves, which are given by the first solution in Equation (19). Regularity in the perturbation series may be lost for the second solution, As the previous analysis suggests, in the large wavenumber limit, there is a vortex mode with zero real frequency. We see that the damping rate is actually twice as fast as previously suggested. This discrepancy is attributed to the fact the perturbation approach is, as discussed previously, not valid in the large vertical wavenumber limit. When the wavenumber is not large, the damping rate is lower. Note also that the solutions lost correspond to the slow waves, the frequency of which is proportional to v A k z ; hence, these represent the singular solutions. A perturbation series, for small collision frequency, is possible to investigate the solutions of Equation (22), resulting in the first two slow solutions in Equation (21). Finally, analytical expressions are available for the frequency in the limit of high collision frequency. The oscillatory solutions are given, to leading order, by In this limit, the vortex modes are purely imaginary. The propagating solutions are purely real in this limit, hence we see that the damping tends to 0 in the limit of high collision frequency. When the collision frequency is high, the plasma acts as a single fluid; the solution here then represents a buoyancy-driven MHD wave in the single fluid composed of ions and neutrals.

Discussion and Conclusions
It has been established observationally that buoyancy-driven IGWs, likely generated through turbulent motion, propagate through the lower solar atmosphere. Magnetic fields play an imported role in the dynamics of solar plasmas. The effect of magnetic fields was studied theoretically by, e.g., Hague & Erdélyi (2016).
Previous works considered, as a first step, that the plasma could be described as a single fluid. This approach is usually applied based on the assumption that the plasma is fully ionized. This is not the case for the lower solar atmosphere. In this article, we have taken the first steps toward providing the theory of buoyancy-driven MHD waves in a partially ionized plasma stratified under gravity. The case of a partially ionized plasma was studied using a two-fluid model, following the approach of Zaqarashvili et al. (2011b) who studied MHD waves in a homogeneous plasma.
The governing equations for buoyancy-driven MHD waves were found. The wave solutions were studied, numerically and analytically, via the local (WKB) dispersion relation. We calculated the corrections to the frequencies and the damping rates in several important physical limits. A key result was that the frequency, derived from the starting point of a fully ionized fluid, is close to that of a weakly ionized fluid when the ionneutral collision frequency is large relative to the wave frequency. The major difference to the fully ionized case is that damping due to ion-neutral collisions is present. We conclude that, magnetic effects may still be very important even to a weakly ionized plasma. In the solar atmosphere, the Brunt-Väisälä frequency is much smaller than the ion-neutral collision frequency, hence this result is applicable to lower frequency IGW-type waves.
When the collision frequency is low, relative to the wave frequency, the two fluids interact weakly so oscillatory solutions are present in each of the gases: buoyancy-driven MHD waves in the ionized fluid and IGWs in the neutral fluid. For large collision frequency, the two gases act as a single fluid. The oscillatory IGWs become purely damped vortex modes and we have MHD waves in the ion/neutral single fluid. This limit is important in the solar atmosphere for higherfrequency waves with a magnetic character. The low collision frequency results may also have application to IGW-type waves in stars such as Red Giants where the Brunt-Väisälä frequency can become larger than typical solar values.
The results derived in this article are likely to be very useful in interpreting and explaining forthcoming observations of oscillations in the lower solar atmosphere. The effect of, for example, a magnetic field on solar atmospheric IGWs has yet to be established observationally. A detailed study is urged to observe the properties of buoyancy-driven MHD waves. Such a study is within the capabilities of current solar telescopes.