Ionization fronts in coupled MHD-gas simulations

Partially ionized plasmas are ubiquitous in both nature and the laboratory, and their behaviour is best described by models which take into account the interactions between the neutral and charged species. We present a new non-linear, 3-dimensional, finite difference Gas-MHD Interactions Code designed to solve simultaneously the time evolution of fluid equations of both species in the conservation form as well as collisional interactions between them via appropriate choices of source term; in particular, we present results from this code in simulating Alfven ionization in a partially ionized plasma. In this fashion, larger changes in the ionization fraction than were addressable in the linear limit are possible. Alfven ionization is shown to impart plasmas with an inherent resistance to rapid recombination, where the recombination itself is significant enough to drive relative motion between the ionised and neutral species at speeds in excess of the critical velocity.


I. INTRODUCTION
Plasmas are ubiquitous in both the laboratory and space environments, and when comparing different plasmas, there is a huge variation in many of the key parameters.Many of these plasmas are not fully ionized with their ionization fraction being dictated by a balance in reaction rates between ionizing reactions, either thermal or otherwise, and recombination reactions.It is clear that these partially ionized plasmas still fit the definition of a plasma by exhibiting the long-range, collective behaviour that we associate with plasmas despite charged particles being, in many cases, the minority species.
The charged particles and the neutral particles in a partially ionized plasma have different governing physics with only the former having any interaction with electric and magnetic fields, both those externally applied and the self fields generated by the motion of the charged particles themselves.However, the neutral and charged components interact via collisions, both elastic and inelastic.As a result, the electric and magnetic fields can indirectly influence the motion of neutral atoms and vice versa.
The partially ionized nature of many astrophysical plasmas has been driving recent efforts to add the additional effects of the neutral species to new models.In the fluid limit, this can be done either by the inclusion of additional terms to represent the neutrals in the model equations of a single plasma fluid or by treating the ion-electrons and the neutrals as two separate fluids with interaction terms.][15] Alfv en ionization (AI) is a mechanism, first proposed in 1942, 16,17 that relies on elastic collisions between a magnetised plasma and an impinging flow of neutral gas to create pockets of charge imbalance which then produce a nonthermal distribution of electrons.If the impinging flow has sufficient kinetic energy, then these electrons can reach the energy required to ionize the gas.For an overview of the mechanism and experimental and observational evidence of AI, see, for example, Lai, 18 Danielsson, 19 and Newell. 20ilson and Diver 21 presented a linear two-fluid finite difference code for partially ionized plasmas that modelled frictional drag between the ion-electron plasma and the neutral gas as well as an Alfv en ionization (AI) term.
Alfv en ionization provides several good reasons to move to a non-linear fluid code: a simple argument is that a kinetic energy dependence exists in AI and kinetic energy is an inherently a non-linear calculation.More practically, nonlinearity is desirable because the most obvious astrophysical applications are in low ionization fraction plasmas where AI can induce large, order-of-magnitude changes in the plasma density such as in the atmospheres of brown dwarfs. 22This places these applications out of reach of a linear simulation.
To model effectively these scenarios, a non-linear version of the gas-plasma interactions code described in Wilson and Diver 21 was written, which itself was based on an existing MHD-gas momentum coupling code from Diver et al. 23 This ground-up approach has the advantage of being tailored specifically for the inclusion of AI, allowing the concepts behind original linear code, given its unique plasmaneutral coupling dynamics, to be built upon.However, there is no reason why appropriate model equations and AI source terms could not be added to a community standard, featurerich MHD codes such as BOUTþþ (Ref.24) or Lare3d. 25he resulting Gas-MHD Interactions Code (GMIC) has its model equations described in Sec.II and the numerical methods in Sec.III before results using GMIC to simulate Alfv en ionization are presented in Sec.IV.

II. MODEL EQUATIONS
GMIC solves the equations of MHD and Euler's equations of fluid dynamics in the conservation form with the addition of a frictional drag term in the momentum equation 23 and an Alfv en ionization term of the form described in Wilson and Diver. 21he model equations for the plasma are 26 q, v, and B are the plasma density, velocity, and the magnetic field, respectively.E is the total energy density which we have defined as The total pressure, P t , is given by, , the sum of the magnetic and thermal pressure where c is the adiabatic constant.We have assumed that the thermal pressure is isotropic.
The model equations for the gas are The ^indicates the corresponding fluid variable for the gas.In both sets of model equations, we have introduced some tensor notation; I is the identity matrix and ab is the dyad of vectors a and b formed by multiplying the column vector a with the row vector b such that ðabÞ ij ¼ a i b j .The system of equations of the gas can reassuringly be recovered by setting B ¼ 0 in the equations for the MHD plasma.
We have introduced source terms to the right hand side of our system of equations that represent the interactions between the two fluids.The first is a frictional drag coupling term, K x;y;z ¼ Cqðv x;y;z À vx;y;z Þ; the momentum coupling is controlled by a coupling constant (C). 23The remaining terms, 6 _ q ai and _ v ai = _ v ai , are the corresponding changes in plasma and gas density and plasma/gas velocity that are dictated by the fluid model for Alfv en ionization.
The formulation of these source terms is described in detail in Wilson and Diver; 21 due to the importance of the AI term to the results that follow, we shall summarise the key points of the approach here.When a neutral fluid flow impinges on a plasma, then neutral atoms will collide with and displace ions.This causes isolated regions of charge imbalance to form stochastically.The maximum potential of the electric field associated with these regions of charge imbalance is the same as the kinetic energy of the incoming neutrals, as at that stage the electric field itself prevents further ions from being displaced.The electrons, unable to immediately follow the displaced positive ions due to presence of a perpendicular (to the direction of the neutral flow) magnetic field, are accelerated under their mutual repulsion up to, and indeed beyond, the electric potential of the charge pocket. 27Some of these electrons will collide inelastically with neutrals and ionize them.
Since ionization can only occur when the electrons are able to reach this ionization energy, the mechanism only takes effect when the incoming neutral flow is at a critical velocity.This critical velocity, assuming only the presence of a single neutral species, is given by where E I and m g are the neutral gas ionization energy and mass, respectively.In laboratory experiments, the observed critical velocity lies very close to the value given by Eq. ( 10). 28,29ach time an ionizing reaction occurs, then the kinetic energy of the free electron is reduced by E I , with the remaining kinetic energy shared between it and the newly created free electron, and an atom that was once part of the neutral species has become a plasma ion.Energy balance is maintained by tapping an internal energy state that is not represented in the model equations.Ultimately, when considering AI, the free electron gained its non-thermal energy from the kinetic energy of the bulk neutral fluid motion and the neutral fluid must therefore have lost this kinetic energy.The mechanism halts when the kinetic energy has been reduced below the point where the pockets of charge imbalance are able to produce ionization energy electrons, or rather when v < v c .A gas moving perpendicular to the magnetic field at velocity v, in the rest frame of the plasma, has only a finite reservoir of energy density, E ai , available for ionization The energy required to ionise a unit mass of gas into plasma is analogous to the energy associated with a more orthodox phase change; just as we have the latent heat for evaporation, we have a latent heat for ionization.In this way, the total number density of atoms that can be ionized by AI for a given relative fluid velocity is just N ¼ E ai /E I and the total mass density ionized is then m gN .Simple substitution with Eq. (10) and recognising that the amount of energy required to ionize a fluid cell of gas with density q is E I q=m g tell us that, assuming 100% efficiency, total ionization is possible when v !ffiffi ffi 2 p v c .While a fluid model is not able to resolve the sub-Larmor radius/period processes involved in AI, our goal is to have a high enough time resolution to be able to resolve a fluid proxy for AI; as a result, we include a scalar Alfv en ionization efficiency factor, f, in order to give a rate of change of the plasma density due to AI to be included in the model equations The value for f can be varied between 0 and 1 and represents to the competing time scales of the fluid and particleparticle interactions.In practical terms, f represents the consequence of the difference between the microscopic scale of ionization and the macroscopic scale of fluid motion.
As the spatial resolution grows, then f is reduced allowing some measure of slippage between the two species before AI acts to reduce the relative velocity.
Although there is no source term in the energy evolution equation, the total energy density does have a contribution from both the thermal energy and the kinetic energy, which both evolve as AI takes place.Energy is conserved but only if the internal energy difference of the freshly created plasma is considered.Note that an atom of gas has a lower internal energy than a corresponding ion, Êinternal ¼ ÀE I ; E internal ¼ 0.
The amount of energy, dE, required to ionise q ai is given by Since the energy required to make this jump in the internal energy (via ionization) comes from the kinetic energy of the flowing neutral, then the perpendicular gas velocity is reduced This term always acts in the direction opposite to the gas velocity in the frame of reference of the plasma.The plasma velocity is also evolving; while the plasma velocity does not directly change due to AI, a newly created ion will have previously had a mean velocity equal to the mean velocity of the gas and not the plasma.The mean velocity of the plasma will therefore move towards that of the gas by an amount proportional to the change in plasma density and, unlike the previous term, this acts on all 3 velocity components.
In this manner, total energy is conserved.

III. NUMERICAL METHOD
From the two arrays of model Eqs.( 1) and ( 7) flux functions are created (see Appendix A), and these flux functions are functions of the entire system of variables (densities, velocities etc.), u.To calculate component-by-component, the divergence of an arbitrary function of a system of variables (r Á gðuÞ) only invokes the x-derivative of the xcomponent of gðuÞ, the y-derivative of the y-component, and the z-derivative of the z-component.For example, the xcomponent flux function will only be operated on by an xdirected finite difference operator, and so on.This makes it a very natural approach to calculate divergence from flux functions.Each flux function has n elements, for the n model equations (8 for the plasma, 5 for the gas).This separation of the model equations into flux functions means that there are now a set of 6 flux functions, one per dimension for both plasma and gas.
These flux functions are solved by a conservative finite difference scheme.The code was initially trialled with two such solvers: the Richtmyer scheme 30 and the MacCormack scheme. 31Both schemes are two-step algorithms equivalent to the Lax-Wendroff scheme for linear problems.First, a predictor step calculates an estimate for the system of variables at the half-timestep before the corrector step, using updated flux functions from the half-timestep, calculates the evolution to the full timestep.
The MacCormack scheme is an extremely popular choice of finite difference scheme.It has been used extensively across a variety of fluid dynamics for example; the simulation of water waves; 32 the breakup of thin liquid sheets; 33 or modelling the sucker-rod pumps of oil wells. 34It is also used extensively in plasma physics; in modelling the MHD pedestal in tokamaks; 35 simulating the solar wind 36 or the Rayleigh-Taylor instabilities in inertial confinement fusion (ICF); 37 and countless other examples.
The extensive use of this scheme, due to its simplicity, allows its strengths and limitations to be well known.][40] The model equations, source functions, the MacCormack scheme itself, and the boundary conditions used are described in Appendix A.

IV. RESULTS
GMIC was used in order to simulate the two-fluid response, including Alfv en ionization in two scenarios.In the first, we set up a high amplitude, cylindrically symmetric overdensity in a gas at a constant temperature resulting in an overpressure which drives outward fluid motion at the local sound speed.This outward fluid motion will result in Alfv en ionization wherever it exceeds the critical velocity in the perpendicular to magnetic field direction.
In the second simulation, we repeat the same initial conditions from the first but with the addition of an equivalent but opposite perturbation in the plasma, an underdensity, such that the total density perturbation is zero.This represents the two-fluid's response to a significant recombination event and the increased relative velocities as the two fluids move opposite directions, which results in significantly higher ionization rates.
The model equations are normalized with respect to the magnetic field which brings out the plasma beta as a free parameter and allows all velocities to be normalized to the Alfv en speed.We set the Alfv en speed to be equal to the basic grid speed such that v A ¼ v 0 ¼ l 0 /t 0 .Plasma and gas densities and the magnetic field are all normalized with respect to their respective equilibriums.The MHD evolution is controlled by a choice for this plasma beta, and the interaction terms are controlled by f, C and the ratio of the two fluid densities.In a standard MHD code, calculating real world values from the output (in code units) can be done by a choice of either B 0 , q 0 , and either l 0 or t 0 .However, since the ratio of the two fluid densities implies a temperature (assuming thermal ionization balance), then choosing a value for B 0 also fixes q 0 and vice versa.Therefore, choosing a value for either B 0 or l 0 and a for either l 0 or t 0 allows the full system of physical variables to be recovered.
In all following results, the equilibrium ratio of fluid densities is kept at 1.0 implying equal parts plasma and gas or, equivalently, an ionization fraction of 0.5.For hydrogen, assuming Saha equilibrium, this would imply a temperature of approximately 9500 K.

A. Ionization from gas explosion
In the first scenario, which we call a cylindrical gas explosion, we place a localized Gaussian shaped density in the centre of the computational domain with a maximum non-dimensional amplitude of 0.4 in the gas.This density overabundance causes a cylindrically symmetric wave to expand outwards at the local sound speed.The plasma is initially stationary, meaning that a relative velocity exists between our two fluids; where this relative velocity exceeds the critical velocity, ionization takes place.The ionization rate is proportional to both f and v rel -v c .
The ionization efficiency factor is f ¼ 0.04 and the critical velocity is 0.03 in units of the Alfv en speed; these values were chosen as typical observable values for illustration purposes.Conventional momentum coupling, where the gas and plasma fluids impart impulse on each other via elastic collisions, is absent from this simulation, i.e., C ¼ 0. This choice is made in order to make the unique dynamics of AI more distinct, but there is no expectation that there could ever be a plasma where the ionized and neutral species interact via AI and not by a frictional drag.The magnetic field is orientated vertically with B ¼ B y ŷ ¼ 1.The equilibrium gas and plasma densities are both 1.0, this equilibrium density is subtracted to obtain a density perturbation (such as shown in Fig. 1), and these densities imply an ionization fraction of 0.5 or a ratio of q=q of 1.0.Both fluids are at the same equilibrium temperature.
The evolution of the perturbation in the gas density is shown in Fig. 1.Since a significant amount of our time dependent gas density will become plasma, the total density provides a good visualisation of the disturbance.In the absence of momentum coupling, all damping of this initial perturbation is a result of AI.
Figure 2 shows a newly created plasma from AI.The higher the plasma density, the more ionization has taken place.Note that this is not ionization rate but rather total integrated ionization as new plasma does not recombine.The localised increase in plasma density due to AI leads to compressional motion in the plasma; the evolution of plasma density seen in Fig. 2 is from a combination of continuing ionization and standard MHD motion.
The location of maximum ionization is no surprise.The nearer to the centre of the domain, the higher the wave amplitude and thus the higher the relative velocity and the higher the degree of the ionization.Since the magnetic field is vertical, there is a symmetry around this axis.Ionization only occurs due to relative velocity in the perpendicular direction; therefore, where the gas motion is mostly in the xdirection, there is a higher rate of ionization.This is seen clearly in Fig. 3 by looking at a plot of relative velocity in the x (cross-field) direction.This is obtained by subtracting the two fluid velocities although in practice the relative velocity is almost equal to the gas velocity for much of the simulation since the plasma takes a while to FIG. 1. Gas density perturbation at select times with C ¼ 0, f ¼ 0.04, v c ¼ 0.03, the magnetic field is orientated vertically.The initial disturbance has a central amplitude of 1.4 compared to the mean density of 1.0 and is shaped as a Gaussian.This maximum amplitude is truncated at t ¼ 1 for clarity and so the full Gaussian shape is not shown.The resultant wave loses amplitude due to dispersion and Alfv en ionization.react due to the lack of elastic coupling and the presence of the magnetic field.
We can see that, by t ¼ 30 although the relative velocity is significantly lower than at t ¼ 10, it remains, for at least some of its span, higher than the critical velocity of 0.03.As such, we should expect ionization to still be occurring at this time.Figure 4 shows the total mass of plasma with the background subtracted.The ionization rate is most significant at early times but continues up until the perturbation leaves the domain.The lack of momentum coupling allows the relative velocity to stay above the critical threshold for longer.The more persistent ionization means that new plasma is spread across a broader area but still peaks strongly near the disturbance.FIG. 2. Perturbation in plasma density due to AI from the driver shown in Fig. 1, i.e., with C ¼ 0, f ¼ 0.04, v c ¼ 0.03, the magnetic field is orientated vertically.The majority of ionization occurs in the centre of the perturbation at early times but continues to some extent until the perturbation exits the domain.FIG. 3. Relative velocity between the gas and plasma in the x-direction (perpendicular direction) taken from the same simulation as Figs. 1 and 2. The colour bar is in units of v A .The total amount of ionization in a given fluid cell is dependent only on the value of v rel .
The total ionization is much larger than can be achieved with our linear code. 21In some places, the plasma density is changed by near 5% of the background plasma density.It is easy to imagine a situation where a region of plasma is being repeatedly perturbed by the consistent wave motion: In such a scenario, the fraction of gas is ionized a little at a time, and if the ionization events are frequent enough, then, due to the different time-scales of the thermal recombination and nonthermal ionization, the freshly ionized plasma is able to persist.In such a scenario, we expect even more significant departures from a thermal ionization equilibrium, greater than the 5% seen here.
In this sense, a partially ionized plasma in the presence of high velocity fluid motion gains some memory of this motion, and the ionization fraction is no longer purely determined by a local ionization and recombination rate equation balance but by a cumulative history of what came before.
Another obvious extension here is that once above the threshold in ionization fraction for AI, which in practice is a very low threshold, the AI rates are not dependent on the plasma density.This means that we would see the same absolute change in the plasma density even with a seed plasma that is many of orders of magnitude smaller in density.In such a case, these same relative velocities from the same perturbation would result in a dramatic change to ionization fraction.
Just as in the linear case, the newly created plasma density causes a fluid plasma motion.We can see that the new plasma is a source of fast and slow magneto-acoustic modes from the separated peaks and troughs evident in later times.Unlike previous linear approaches, we have a disabled momentum coupling so we are seeing the presence of MHD modes rather than MHD-acoustic coupled waves.
Interestingly, the inhibited motion of the plasma across the field means that a large density gradient persists for some time.The simulation was run up to t ¼ 200 or 2000 finite difference timesteps.This corresponds to 5 Alfv en crossing times.By the end of the simulation, the freshly ionized plasma has not dispersed significantly and what redistribution has occurred is mostly in the parallel-to-field direction.

B. Ionization from implosion-explosion
We repeat the previous simulation for a slightly more sophisticated set of initial conditions; all simulation parameters are the same as IV A except for the initial perturbation.Instead of just an over-abundance in the gas, we also include a plasma rarefaction; this rarefaction of the plasma is equal and opposite to the increase in the gas density such that, at t ¼ 0, the total density is 1.0 across the entire domain.This is a combination of the cylindrical gas explosion in Sec.IV A and the inverse, a negative cylindrical implosion in the plasma, which consists of a density rarefaction of the same shape and opposite direction as the previous case.
The idea here is that a total density perturbation of 0 could represent a region where plasma has recombined into gas, simultaneously causing pressure gradients in both the gas and the plasma and driving flows in both fluids.These flows result in relative velocities between the species as plasma flows in and gas flows out from the recombination event.As before, the gas is still moving outwards at the local sound speed, but this time it encounters an inflowing background plasma, which greatly enhances the relative velocity compared to the gas explosion.If the initial recombination is rapid enough, then these flows can result in AI which would re-ionize some fraction of the gas.
In a sense, this may offer plasmas in a particular regime some form of resistance to rapid recombination.When a region of plasma recombines into gas, then the changes to the partial pressures of both fluids will naturally lead to relative velocities and thus some degree of re-ionization.If this ionization is significant when compared to the original perturbation, then we could describe such a plasma as resilient against recombination.
We see the initial perturbations to the plasma density, gas density, and the total density in Fig. 5.
The developing magnetoacoustic perturbation of the plasma (not present initially) makes the ionization signature hard to contrast.Instead, we look at the total amount of ionization done over time.Figure 6 shows the comparison of the total ionization done between the implosion-explosion and the gas explosion simulations.
The total ionization is significantly higher in the implosionexplosion case, causing approximately 3 times more ionization than previously.As the plasma initially moves one way and the gas the other, there are significantly higher relative velocities resulting in the enhanced ionization.
We choose a normalized system of units where the initial integrated plasma density is-100-negative because relative to the equilibrium some of the plasma has recombined, lowering the total plasma density.On this scale, the plasma density, when the first magnetoacoustic wavefront exits the domain at t ¼ 25 (and thus total plasma density becomes a less useful measure), was -91.8.This corresponds to an increase of plasma density of 8.2 of these normalized units over the whole simulation.In other words, this means FIG. 4. Total ionization from gas explosion.Ionization routine is disabled until t ¼ 2 at which point the ionization rate is high.At late times, ionization is still taking place, but the rate has significantly slowed.AI was able to restore 8.2% of the initial deficiency in the plasma density.The amount of total ionization that takes place is increased by lowering the critical velocity threshold or by increasing f; a higher f means more of the total wave energy is extracted before the wave disperses.The total ionization, for both this case and the previous case, can be seen in Fig. 6.This total ionization is low when compared to the initial perturbation but not insignificant.The plasma is inhibited from moving in the x-direction due to the force exerted on it by the magnetic field.Alfv en ionization can only occur due to relative velocity in the perpendicular to field direction and the plasma is unable to contribute as much to the relative velocity in the x-direction as the gas due to the magnetic field.If the perturbation were spherically, rather than cylindrically symmetric, then the total ionization would be greater.In a spherically symmetric setup, two out of three of the dimensions are perpendicular to the field, assuming that equipartition amongst the three dimensions 2/3 of the energy goes into perpendicularly directed motion (compared to 1/2 in the cylindrical case).
In our test case, we were able to re-ionize a reasonable fraction of the initial perturbation, but AI was not close to being able to re-ionize all of the initial recombination globally.If rapid recombination occurs, high amplitude compressional motion may cause AI which in turn acts to slow, but not prevent, the recombination.Note that while the average global level ionization might be modest compared to the initial perturbation, localised changes to the ionization fraction are clearly non-linear.

C. Effect of varying beta
If the plasma beta is increased by reducing the magnetic field, then the motion of the plasma in the direction perpendicular to the equilibrium magnetic field is less restricted and the ionization rate increases due to the greater relative velocities of contra-flowing plasma and gas.This behaviour is shown in Fig. 7.At low beta, less ionization takes place due to the inability of the charged component to flow in the perpendicular-to-field direction; the ionization increases with increasing beta.At low beta, ionization is not reduced to zero since the neutrals are always free to move in the perpendicular direction, regardless of magnetic field strength.At FIG. 5. Initial conditions of the gas and plasma densities, with equilibrium subtracted, for the implosion-explosion setup relative to the equilibrium density of 1 for both.Note that the plasma and gas initial perturbations are equal and opposite giving a total perturbation of 0, the magnetic field is orientated vertically.FIG. 6.Total ionization from the implosion-explosion compared to the gas explosion case.Note that all conditions are as in Figs. 1 and 2 other than the initial perturbation.Ionization routine is disabled until t ¼ 2 at which point the ionization rate is high.The implosion-explosion significantly outperforms the gas explosion, resulting in a factor $3 more ionization at all times.Note that the total elapsed time, t ¼ 25, corresponds to the maximum time for which the perturbation remains inside the computational domain.
high beta, ionization does not increase without limit since while at higher beta the plasma motion is no longer being controlled by the magnetic field, there are diminishing returns to this behaviour: when the plasma beta is already high, each further increase provides less of an increase in the magnitude of the perpendicular-to-field plasma motion and thus to the maximum potential for relative velocities in this direction.The limit here is one of infinite beta, where the plasma behaves like a gas and perpendicular and parallel motions are equivalent.Eventually, as the magnetic field is reduced below a threshold value, there is no Alfv en ionization at all.
For a beta of 2, the total ionization done is $25% of the initial perturbation.So AI even in a reasonably high beta plasma does not manage to completely cancel out the initial recombination.
To see why, suppose that the plasma beta is infinite, the critical velocity was 0 and f ¼ 1, then all energy in the perpendicular direction would be claimed for ionization.Since, due to the symmetry of the problem, 1/2 of the kinetic energy (in 2-dimensions) is in the perpendicular direction, then, in such a scenario, 1/2 of the total perturbation energy would be converted to plasma internal energy.In 3-dimensions, the ratio is slightly improved, assuming the same symmetrical disturbance in 3-d, only 1/3 of the total perturbation is in the parallel to field direction with 2/3 of it in the perpendicular directions.
If we assume that the recombination triggered by the implosion-explosion scenario resulted in radiation that escaped the system, then it is not possible that the total perturbation energy is greater than the internal energy deficit that was caused by the initial recombination.It is therefore impossible for the AI to completely undo the initial recombination because energy has been irredeemably lost as radiation.The maximum potential for ionization is also reduced by the presence of collisional drag between the two fluids: collisional drag becomes increasingly detrimental to AI the closer the typical relative velocities are to the critical velocity.
As beta decreases, then the kinetic energy of the plasma becomes more and more biased to the parallel to field direction and the total possible ionization decreases.If C, the momentum coupling coefficient, is non-zero, then the relative velocity is reduced without ionizing (and if f < 1) and the wave will damp and disperse without complete extraction of the kinetic energy for the purpose of AI.We examine this next.

D. Interplay with momentum coupling
To investigate the importance of the strength of momentum coupling, the implosion-explosion scenario is repeated as before, with b ¼ 1, f ¼ 0.04, v c ¼ 0.03.We vary C from 0.001 to 0.05.The presence of momentum coupling reduces the total Alfv en ionization by reducing the relative velocity between species.
Figure 8 shows the total ionization done in the implosionexplosion case for varying C (the momentum coupling coefficient).We can see that an increase in the momentum coupling has a corresponding decrease in total ionization.This trend continues until the ionization reaches a floor at C ¼ 0.01 and the total ionization done of $2% of the perturbation.The total ionization done never reduces to zero.The limiting behaviour is a consequence of the finite mesh ratio: even in perfectly coupled fluids, we only equilibrate the two species velocity fields every half timestep.In between this equilibration, with the governing equations for the two species being different, there is a chance for more relative velocity to build up between each half timestep; in essence, this motion takes place before either species can communicate to the other that they are moving.If this relative velocity exceeds the critical velocity, then ionization still takes place even in strongly coupled fluids.The magnitude of the relative flows can only be reduced to zero for increasing values of C; if sufficient mesh resolution is employed (beyond the mesh-ratio demonstrated here), then for C ¼ 0.5, the relative velocity will vanish after one timestep.FIG. 7. Total ionization from the implosion-explosion changes as the plasma beta is varied, all other parameters are fixed.Higher beta corresponds to higher total ionization.The total ionization is taken from time t ¼ 25.

FIG. 8. Total ionization from the implosion-explosion for changing C.
Higher momentum coupling corresponds to lower total ionization.The total ionization is taken at t ¼ 30.f is fixed at 0.04.Despite this being a numerical effect, there exists a scale length, at the scale of the Larmor radius/ion-neutral mean free path, at which the two fluids cannot have zero relative motion but, of course, such scale-lengths are unresolvable in the fluid limit.
There is an extremely rapid falloff in total ionization caused by increasing the momentum coupling constant from 0 to 0.01.We attribute this dramatic change to the low value of f (being fixed at 0.04).Since, in each timestep, both f and C reduce the velocity by an amount proportional to their values, then the ratio of C to f dictates what fraction of the relative velocity goes towards frictional drag instead of towards AI.The exact fraction of velocity that is extracted to ionize is not as simple as the ratio of these two numbers because ionization only works to reduce relative velocities when they are above the critical velocity, whereas momentum coupling works to reduce relative velocities at all times.The ratio of relative velocity damped by AI and relative velocity damped by momentum coupling is only ever at most the ratio of f/C and this limit is reached only when v rel ) v c .When dealing with velocities only slightly above the critical velocity, then the vast majority of relative velocity is lost to frictional drag rather than AI even when f is much larger than C and so the total ionization, in this regime, is not strongly dependent on the ratio of f/C.This is seen in Fig. 8; by the time C is 5 times higher than f, the majority of the relative velocity is invested in adjusting the overall flow elastically, rather than ionizing the gas.
Two colour maps of plasma density, one for C ¼ 0 (Fig. 9) and one for C ¼ 0.01 (Fig. 10), show clear differences in the plasma behaviour.Both of these snapshots show t ¼ 30 in their respective simulations.There are some similarities between these figures in that the perturbation on the plasma density remains large for some portion of the plasma at this time.It is hard for the plasma to diffuse quickly into the hole created by the recombination, since the cross field motion is inhibited.This is partly compounded by the slightly slower (20%) speed of sound in the low density regions.
These figures highlight distinct differences in the ionization potential of the two cases.In Fig. 9, the central region of plasma rarefaction has been partly filled in by freshly ionized plasma created by AI resulting in a much higher plasma density in this area compared to the coupled case (Fig. 10).The central region is where the fluid velocities are highest and the majority of ionization takes place.
Figure 10 in contrast shows that the central region remains at a relatively lower plasma density.In this case, the plasma motion is mostly confined to the vertical, parallel direction.While coupling has not reduced the AI rate to zero, due to the neutral component still being able to provide perpendicular velocities greater than v c , it is significantly reduced by the presence of coupling.

V. CONCLUSIONS
We have created a novel non-linear Gas-MHD interaction code (GMIC) and demonstrated its efficacy by simulating new scenarios that have eluded effective modelling until now.GMIC simultaneously solves the fully non-linear fluid equations of an MHD plasma and a coexisting neutral gas by the method of finite difference using the MacCormack method.Appropriate choices of source terms also allow this code to include the collisional interactions between these species in terms of both elastic collisions (momentum coupling) and inelastic collisions (ionization).In particular, the fluid velocities and magnetic field information available in a fluid model allow Alfv en ionization rates to be calculated.While the results presented are 2-dimensional, to capture the FIG. 9. Perturbation in plasma density for the implosion-explosion with b ¼ 1, f ¼ 0.04, and C ¼ 0. Note that the centre region has been filled in significantly more than the coupled case by freshly ionized plasma.The magnetic field is orientated vertically.This snapshot is taken at t ¼ 30.essential parallel and perpendicular behaviour of the system, GMIC has the capacity for full 3-dimensional simulations.
The non-linear nature of this code allows large amplitude and velocity scenarios inaccessible to previous linear approaches.We tackle the expansion of a cylindrically symmetric Gaussian overdensity of gas into a background seed plasma, resulting in large changes in ionization fraction in localised regions.The resultant ionization is strongly asymmetric due to the asymmetry in both the Alfv en ionization rate equations and the plasma response.
Since both elastic momentum coupling and Alfv en ionization depend on, and act to reduce, relative velocities between the ionized and neutral species, they compete.We show, unsurprisingly, that the more strongly the two fluids are coupled together, the less effective Alfv en ionization is.In simulations using discrete timestepping, momentum coupling must always take a finite time to be established, and hence there is always the possibility of ionizing flows.Although this limit is numerical, in a real plasma a similar limit also exists; if we were to have arbitrary spatial resolution of a plasma, then there will be some small scale length below which relative motion is not yet damped by elastic collisions.At this scale length and below, there exists an exploitable AI energy reservoir.
We find the total ionization that occurs due to AI depends strongly on the plasma beta.At low beta, the magnetic field heavily limits cross field plasma motion.This means that almost all the relative velocity in the perpendicular direction comes purely from gas motion.At high values for beta, both fluids are able to move in the cross field enhancing the relative velocities and total ionization seen.
By simulating an extreme recombination event where a significant portion of a high ionization fraction plasma recombines leaving an overdensity in gas and the opposite under density in plasma, we drive motion in the two fluids, via their partial pressures, in opposite directions.This dramatically enhances the ionization rates and is shown to undo a significant fraction (8%) of the initial perturbation.We liken this to such a plasma as having some kind of innate resistance to rapid recombination.
We anticipate that this method of capturing Alfv en ionization rates via a source term in a system of fluid equations will be useful in simulating the dynamics of partially ionized, magnetized plasmas such as the atmospheres of brown dwarfs and gas giants.numerics and, for the purposes of demonstrating Alfv en ionization and gas-MHD momentum coupling in a non-linear regime, we only require simulations involving smoothly varying solutions.

Flux functions
The flux functions for the plasma are given by ; ;

:
(A5) The flux functions for the gas are given by

Source functions
Unlike flux functions which allow evolution of a conserved quantity due to a difference in the incoming and outgoing passage of said quantity through the cell boundary, source terms allow an external term to create/destroy some of this quantity without taking it from somewhere else.For example, ionization can increase the plasma density in a cell without lowering it in another.
The source functions for the plasma are  Note that the source terms here are mostly symmetric, the same amount of momentum gained by one species is lost from the other, and each unit mass of plasma created comes at the expense of a unit mass of gas.The AI term is an exception, the momentum of the gas is reduced, but the energy is gained by the internal energy of the plasma.
There is a final adjustment to the fluid velocities made whenever AI takes place; fresh ions will have been created from neutrals and, as such, will have velocities sampled from the distribution of velocities of the gas, not the plasma; the mean of these distributions is different.As such, in order to conserve momentum across species, newly ionized plasma comes into existence with the fluid velocity of the gas and we carry out a weighted average of the new to old plasma in order to get an updated mean plasma velocity.

Boundary conditions
Numerical simulations often require that computational boundaries are constructed between a region of interest (ROI) and other regions one hopes to neglect.Artificial boundaries, especially those drawn through non-uniform flows, that do not affect the region of interest (and preserve the same order of accuracy as the solver of the governing equations) can be tricky to create.It is also always preferable to minimize the computational requirements of the boundary by keeping them as straightforward as possible and fitting them close to the region of interest.
When the governing equations are linear, concerned only with small amplitude fluctuations around a steady mean, then analytic solutions for the creation of a non-reflecting boundary condition (NRBC) are possible.Theoretical solutions for this have existed for many years (for example, Givoli; 45 Tsynkov, 46 and Hagstrom 47 ), but problems in implementation can still be present.
In contrast, when the fluid cannot be represented as linear perturbations around a static mean, there is very little chance to develop an NRBC.It is possible to proceed with the techniques developed for the linear fluid, but this presents problems, the accuracy of the results will be diminished, and the complexity of implementing these boundary conditions may not justify their use.Further, the discretisation of the equations for computation will result in instabilities in the non-linear case that can hinder efforts. 48t is thus often preferable to use more simplistic, ad hoc methods.Sacrificial regions (SR)-regions surrounding the region of interest (ROI) that are governed by modified versions of the physical equations-can be created in order to minimise the effect of the boundaries on the interior at the cost of additional computational time.The equations in the absorbing layer are written as where r is 0 in the ROI and has a positive value in the SR.Any disturbance that reaches the sacrificial region will be exponentially damped and, if signal reaches the computational boundary, any reflections will be similarly damped on their return to the ROI.Since the boundary between r ¼ 0 and r > 0 is itself reflecting, then we must be careful that sigma is increased sufficiently slowly such that the impedance is as closely matched as possible on either side of the boundary.This minimises the reflection and non-physical distortion of any signal approaching the boundary (be that compressional modes or perturbations in the magnetic field), which acts as if the transition into the boundary layer is the same as entering any other fluid cell.The presence of a sacrificial region provides an effective but computationally inelegant solution. 49e choose an exponentially increasing form for r.In one dimension, if x a and x b are the lower and upper limits on the ROI such that the SR is present on grid cells with x < x a and x > x b , then and where A ¼ 0.2 and B ¼ 0.04 are constants chosen to get the desired effect.The computational cost of adding this layer can be reduced by making the layer as small as possible.In the 2dimensional simulations presented here with a total grid size of 600 Â 600 (consisting of a 400 Â 400 ROI with a SR of 100 cells on all sides), the total runtime was trivial (20-30 min on a 16 core intel Xeon E5).The cost of the SR is increased significantly in 3-dimensional simulations.

FIG. 10 .
FIG.10.Perturbation in plasma density for an implosion-explosion with b ¼ 1, f ¼ 0.01, and C ¼ 0.01.Momentum coupling has reduced the ability for the plasma to be re-ionized by Alfv en ionization, which means a lower minimum in the centre of the domain.The magnetic field is once again orientated vertically.This snapshot is taken at t ¼ 30.