Alfvén ionization in an MHD-gas interactions code

A numerical model of partially ionized plasmas is developed in order to capture their evolving ionization fractions as a result of Alfven ionization (AI). The mechanism of, and the parameter regime necessary for, AI is discussed and an expression for the AI rate based on fluid parameters, from a gas-MHD model, is derived. This AI term is added to an existing MHD-gas interactions' code, and the result is a linear, 2D, two-fluid model that includes momentum transfer between charged and neutral species as well as an ionization rate that depends on the velocity fields of both fluids. The dynamics of waves propagating through such a partially ionized plasma are investigated, and it is found that AI has a significant influence on the fluid dynamics as well as both the local and global ionization fraction.


I. INTRODUCTION
In partially ionized plasmas, the interactions between the gas and the plasma particles play a significant role in the dynamics.Since neutral atoms and ions behave differently, with only the latter subject to the magnetic Lorentz force and with their collisional operators being of different forms, a plasma model that incorporates these differences can be necessary to accurately describe a partially ionized plasma.
Interest in the inclusion of ion-neutral interactions in MHD models of partially ionized plasmas has been increasing in recent years partly due to the low ionization fraction of many astrophysical plasmas.In solar physics, there is significant interest in the damping of MHD waves in the solar chromosphere and corona, [1][2][3][4] in the emergence of magnetic flux 5 or in magnetic reconnection. 6,7[10] Single fluid models of partially ionized plasmas have been shown to be insufficient in certain regimes.Zaqarashvili et al. 11 show that a two-fluid model is necessary when studying time-scales shorter than the ion-neutral collision time.Other models have derived a two-fluid model of a partially ionized plasma from the moments of the Boltzmann equation including interspecies interactions. 12,13n a 1942 paper, Alfv en proposed an explanation for the differing compositions of planets to be a result of an, at the time, unrecognised mechanism for the ionization of a neutral gas flow that encounters a magnetized plasma, 14 this mechanism is now known as the Alfv en ionization (AI) or critical velocity ionization effect.AI is mechanism that enables the ionization of the neutral component of a partially ionized plasma by relative motion, between the ionized and neutral species, with kinetic energy greater than the ionization energy. 15This mechanism relies on the presence of a seed plasma and for the relative velocity to be of sufficient magnitude in a direction perpendicular to a magnetic field.It is thought to play a role in the atmospheric composition and the ionization fraction of stellar atmospheres. 16,17r a fluid composed of a single species, the velocity at which we would expect AI to take place is found by equating the ionization energy and kinetic energy.This defines a critical velocity, v c : the speed where the kinetic energy per gas molecule exceeds the ionization energy of that same gas 16 Experimental evidence for such a mechanism is abundant; Fahleson 18 used a "Rotating plasma device" which exploits an E Â B drift to drive ions through a background of neutral gas.Their results show a limit near the theoretical critical ionization velocity (CIV) which indicates that, in their experimental conditions, the plasma could not be driven through the neutral gas at velocities greater than the CIV.This limit can only be exceeded once the plasma is almost entirely ionized meaning total ionization is associated with a sharp rise in plasma velocity.Many other similar laboratory experiments subsequently verified the presence of a critical velocity in several different kinds of discharges (Angerth et al., 19 Lehnert et al.; 20 or reviews by Danielsson 21 and Sherman, 22 for example) and also in industrial discharges, such as thin film deposition magnetron plasmas. 23here is also a number of results from astrophysical observations that show the existence of a critical velocity ionization process.The effect is thought to have been observed in cometary comae, 24,25 in the interaction between Venus and the solar wind 26 and in the magnetosphere of Jupiter's moon, Io. 27 It is also thought to be important in understanding the solar wind and its composition 28,29 and in the atmospheres of stars, both in the Sun 16 and in the cooler atmosphere of brown dwarfs. 17either the limit of relative velocity past v c or the change in ionization fraction of a partially ionized plasma due to these conditions are replicated by any existing MHD codes.The inclusion of AI in an MHD framework allows the previously fixed ionization fraction to evolve in a nontransient fashion; the relative fraction of gas to plasma could always change dynamically as waves and bulk motions move through the gas-plasma mixture.The energetics of this non-equilibrium process are unique in that the energy used for ionization is provided by the kinetic energy of a bulk fluid motion; this is in contrast with the energy present in an electron distribution function for thermal ionization.The obvious consequence being that a partially ionized plasma where the ions are in motion w.r.t. to a gas, whatever the nature of that motion is, can have a different equilibrium ionization fraction than one at rest.Diver et al. 30 used a two-fluid MHD model where the MHD plasma constitutes one fluid and the neutral gas the other.The two fluids interact by a coupling term in the momentum equations of both fluids.The model equations of the two fluids are simultaneously solved, and the frictional drag term between the two species allows the velocity field of one fluid to influence the other.
The novelty in this paper is we have extended this numerical model with the addition of an AI term.This term allows kinetic energy from fluid motion of either species to be exploited in order to ionize neutral gas atoms via a critical velocity, AI mechanism as long as certain prerequisites conditions are met.To do so we introduce an Alfv en ionization term controlled by an Alfv en ionization efficiency factor, f, that provides a link between the microscopic, sub-Larmor radius physics that results in ionization and the macroscopic, fluid scale physics of MHD.
The addition of this term results in a two-fluid code where the gas and plasma not only interact via standard momentum coupling but also by AI which acts to reduce relative velocities in the perpendicular to the magnetic field direction whenever they exceed a critical velocity limit and controls corresponding evolution in the fluid densities of both species by ionization of the gas.The result is a feedback mechanism where fluid motions dictate ionization rates, and the density perturbations resultant from ionization drive further motion in both fluids.
The AI mechanism and the prerequisite conditions are described in Section II, the fluid approximation of this mechanism is derived in Section III, and results and discussion of the effect of this additional term in 2D finite difference simulations of a partially ionized plasma are presented in Section IV.

A. AI mechanism
Alfv en 15 attempted to explain the mechanism by suggesting that pockets of charge imbalance produced by stochastic collisions could result in accelerated electrons.The process he proposed was very simple: a neutral gas flows, at speed v, and impedes upon a stationary plasma with a perpendicular field of sufficient magnitude to magnetize the electrons.Collisions between the neutral atoms and the stationary plasma ions displace the ions from the surrounding electrons, creating a charge separation.Given the magnetic field, the electron transport is impeded by their relatively small Larmor radius compared to the ions, and so the charge imbalance persists and grows until the potential is equal to the kinetic energy of the incoming gas.At this point, collisions have insufficient energy to displace further ions against the temporary electric field.This maximum energy E al is given by, The maximum lifetime of these pockets of charge imbalance is of the order of the ioncyclotron period as this is the expected time for displaced ions to return to their initial positions.
A schematic of the formation of a charge imbalance by impinging neutrals is given in Figure 1.
Once a charge imbalance is formed in this way then electrons can be accelerated.Particle-In-Cell (PIC) simulations by MacLachlan et al. 31 show that a charge imbalance with electric potential E al /e is able to accelerate a fraction (10 À3 ) of the electrons present to E al and above.Consequently, if such a pocket of charge imbalance has an electrostatic potential equal to the ionization energy then it is capable of producing an electron distribution function which includes a significant fraction of electrons at or above the ionization energy.The presence of a magnetic field helps this acceleration process; by restricting the expansion of the electron cloud under their self-field to be mostly one-dimensional (along the field) then the energy can be preferentially imparted to a small subset of the electrons, the ones which have a net force mostly in the parallel magnetic field direction.This anisotropic acceleration produces a larger population of highly energetic electrons compared to the isotropic expansion that occurs in the absence of a magnetic field.Additional processes have been proposed since, mostly depending on the formation of equivalent pockets of charge imbalance.Lehnert 32 showed that newly created ions formed in a magnetic field at rest relative to a background flow can create spatially varying potential structures with typical size of the order of the ion gyroradius.These potentials would also be capable of accelerating electrons crossing the charge imbalance to energies comparable to the ion kinetic energy.Other mechanisms, that do not rely upon charge imbalance pockets, such as a modified two stream instability 22,33 have been proposed as alternative ways to explain the electron heating.For an overview of these mechanisms, see the review by Lai. 34

B. Criteria for AI
All proposed electron heating mechanisms share a set of similar characteristics and these must be kept in mind in deriving a fluid description of the process.
• The AI process has a minimum seed plasma number density before it proceeds. 35 The process invokes a perpendicular magnetic field of sufficient strength.AI is only observed when a component of the relative motion between the charged and neutral species is perpendicular to a magnetic field.Experiments have shown the field must be strong enough such that the electron gyro frequency is comparable to the plasma frequency.If the field either has no perpendicular component, w.r.t. the flow direction, or is sufficiently weak then the AI mechanism is no longer observed. 17,36• The plasma must initially have electron density large enough that the maximum electric potential (after all ions are removed from a Debye sphere) can accelerate electrons to the ionization energy of the neutral.• The process must transfer energy from the flowing neutrals to the electrons, via the ions.At the flow energies for which AI is observed the cross sections for ionization by a neutralneutral collision or ion-neutral collision are far too low. 21he cross-section for ionization by an ion-neutral collision does not reach a value comparable to the electron-neutral collision until the velocity of the ion is comparable to that of the electron, this only occurs at energies several orders of magnitude above the ionization energy. 37,38 Ionization proceeds only when the relative velocity between the species exceeds the threshold, critical velocity.
Since the Debye length defines the scale below which charge neutrality can be violated then, if we require a localized pocket of charge imbalance to form, we require that the Debye length, k D , to be much smaller than the Larmor radius of the ions, r li k D ( r li ; (2) This inequality is equivalent to saying v 2 a =c 2 ( 1, i.e., that the Alfv en speed is small compared to the speed of light.
The perpendicular magnetic field, B, must be of sufficient magnitude to magnetize the electrons; we require the electron cyclotron frequency to be greater than the electronneutral collision frequency eB m e ) n g r g v th ; (5) where n g , r g , and v th are the neutral density, cross section for collisions between neutrals and electrons, and the thermal speed, respectively.
If the collisional frequency is too high then electrons will not remain in the region of charge imbalance long enough before they undergo a random walk outside the potential well.
The electron density must be large enough such that the maximum potential of a charge pocket is greater than the ionization potential.If we assume our charge pocket has a radius equal to the ion Larmor radius then the solution to Poisson's equation becomes This gives a threshold of electron number density At the critical velocity of v ?¼ v c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2eU i =m g p and with m i ¼ m g this becomes which places an additional restriction on the electron density.
If a partially ionized plasma satisfies the constraints in Equations ( 4), (6), and (9) then once the relative velocity between the plasma and neutrals in the perpendicular direction exceeds the CIV threshold then ionization will take place.Electrons are stochastically accelerated up to a maximum of the kinetic energy of the incoming neutral flow which leads to a non-thermal electron distribution function.Collisions between the high energy tail of this distribution and neutral atoms result in ionizing reactions, increasing the ionization fraction.

III. AI FLUID MODEL
The AI mechanism is concerned with particle scale, kinetic processes: ion-neutral collisions, acceleration of electrons under their mutual self-repulsion and ionizing electronneutral collisions all occur at length and time-scales incompatible with the fluid description.None of these processes can be present formally in an MHD description of a plasma.
AI takes place at a small scale, not much larger than the ion gyroradius.It takes place quickly, sub-nanosecond timescales are typical of both the acceleration of electrons into the non-thermal distribution but also typical of the electronneutral collision rate.
However, the criteria laid out in Section II B all depend solely on macroscopic quantities; quantities that are contained in a two-fluid description of the plasma, B; q p ; q g ; v p ; v g .
An MHD model provides the magnetic field, plasma density, and plasma velocity, and a hydrodynamic model provides the gas density and velocity.A combined MHD-Hydrodynamic description allows the creation of a physically consistent model for AI.The numerical code described in this paper simultaneously simulates the evolution of a coexisting gas and MHD plasma and predicts the presence and the magnitude of critical velocity ionization by an appropriate choice of source terms.
The two fluid system of equations are from Diver et al. 16 which are the linearised equations of ideal MHD and Euler's equations of fluid dynamics.These were normalized using three parameters: the ratio of plasma sound speed to Alfv en speed (r), the ratio of the gas sound speed to Alfv en speed (s), and the ratio of the gas density to the plasma density (g).If the gas temperature and plasma temperature are equal then g is no longer a free parameter since, in such a case, it would be uniquely defined by the relative values of r and s.The lengths and times are normalised to the sound speed such that x 0 ¼ c s t 0 .The usual finite-difference mesh ratio k, given by the ratio of the time step to spatial step in these units (k ¼ Dt=Dx), governs the dynamical evolution on the mesh; it is necessary to satisfy k < 1 to maintain stability via the CFL condition. 39dding AI in the fluid limit presents a unique challenge.Although the fluid velocities in MHD are mean velocities that represent an underlying distribution, we shall assume that all particles (plasma or neutrals) are moving at their respective fluid velocity.Since we are operating in a regime where velocities at or above the critical velocity are typical then this assumption will be valid as long as the thermal energy of the plasma is small compared to the ionization energy (k B T e ( E i ), what we shall call cold.When in this regime the random thermal velocities of the gas molecules are small compared to the critical velocity for ionization.AI occurring at close to the theoretical critical velocity, from Equation (1), should be expected for such plasmas but as the temperature increases then this simple expression may become inappropriate.Since plasmas with a temperature that is a significant fraction of the ionization energy are already fully ionized (e.g., from Saha 40 ) then AI is already not relevant for plasmas which cannot be considered cold.
The inclusion of Alfv en ionization to an existing twodimensional, two-fluid finite difference model proceeds as follows: At each timestep, the fluid code advances giving updated velocity and density fields.These fluid parameters are used to calculate an ionization rate.First, we calculate an AI energy reservoir, E ai that comes from the difference in the kinetic energy of the gas, in the frame of reference of a stationary plasma, at the relative velocity and the kinetic energy at the critical velocity.The reason for this reservoir is that once the relative velocity drops to below the critical velocity then AI must cease as the two fluids are no longer capable of generating sufficiently high potential pockets of charge imbalance where q g is the gas density and v rel is the relative velocity in the perpendicular to the magnetic field (e.g., y) direction, given by v rel ¼ jv p;y À v g;y j.This energy reservoir represents the total amount of energy that could be utilised for AI, if the process was 100% efficient; ionization can only take place when E ai > 0. The number of particles that are ionized by this energy source per unit time is where E i is the ionization energy.We have introduced f as the ionization efficiency factor which has dimensions of time À1 and represents the fraction of the maximum possible ionization that is performed in a given fluid timestep.
There are four major factors that contribute to the value of f: (1) The number of charge-imbalance pockets formed per unit fluid cell per unit time, f 1 ; (2) The number of electrons per pocket, f 2 ; (3) The fraction of electrons within a given pocket that are accelerated to energies exceeding the ionization threshold energy, f 3 ; (4) The probability that an ionization energy electron undergoes an ionization reaction in a fluid timestep, f 4 .
The numerical value of f is then given by the product The final term, f 4 , clearly has a value close to unity (the plasma is highly collisional) and so f is determined largely by the first three factors.MacLachlan et al. 31 suggest that f 3 is of the order of 10 À3 .The second factor, f 2 , is well defined as the ambient plasma electron density multiplied by the volume of a pocket.The first factor is proportional to both the volume of a cell and the length of the fluid timestep and, as such, is determined by the fluid simulation parameters.This factor can be greater than unity because the pocket lifetime is much smaller than the fluid timestep.Given that the fluid timestep is at least an order of magnitude greater than the cyclotron period and given the desire to avoid the situation where the AI rate equations are causing non-linear changes in the fluid velocities in a single timestep; it is reasonable to chose a value for the first factor such that f % 10 À2 .
If f ¼ 1 then the entire energy reservoir would be exhausted in order to ionize an appropriate number of atoms and the relative velocity would be reduced to the critical velocity in a single fluid timestep, along with a significant perturbation to the plasma and neutral gas densities.
Since the experimental results of Fahleson, 18 Danielsson and Brenning, 42 and others show a limit, at the critical velocity v c , are unable to drive further relative velocity between the two fluids then their experiments were consistent with f being close to unity when f 1 is calculated using the entire size of their apparatus.However, we do point out that the fluid dynamics are unresolved in these experiments, they only measured the mean velocity across of the entire device.With arbitrary resolution, there must be some scale at which a hard limit of v c is violated.
Equation ( 11) is then used to calculate a rate of change of density, _ q p ai , The ionization from AI introduces some deviation in the electron distribution function from thermal equilibrium that results in a momentary excess of electrons at energies that can undergo ionizing collisions with neutral atoms.The electron distribution function then relaxes back to a new equilibrium.The electron relaxation time is shorter than the fluid timescale in an MHD framework allowing this to take place before the fluid can react. 41he relative sizes of the ionization efficiency factor, f, and C, the momentum coupling constant, tell us the relative portion of kinetic energy that goes into AI and frictional drag.f/C essentially dictates how many atoms are ionized before relative motions are damped by ion-neutral collisions.
In Section IV, to explore the principle; the choice of f is motivated partly by a desire to investigate the changes in a partially ionized plasmas dynamics caused by AI and not simply to calculate total ionization rates.For a given set of plasma parameters (temperature, density, etc.) and simulation parameters (timestep and spatial resolution) there is an appropriate choice of f.Equation ( 13) allows the fluid densities to be updated; each increase in plasma density has a corresponding decrease in gas density such that total density is preserved where the prime notation, 0 , denotes the new species density resulting from AI.In order to maintain energy conservation, the energy used for ionization is extracted from the velocity field of the gas.This can be pictured as being a result of the gas atoms that come to rest in Figure 1 lowering the mean velocity of the gas fluid as the fluid transfers some of its kinetic energy to electrostatic potential where dE is the energy required to ionize the quantity of gas dictated in Equation ( 14).The 7 indicates that a negative velocity is increased and a positive velocity is decreased such that the absolute velocity always decreases.The ionization rate tends to 0 as the relative velocity tends to the critical velocity.
The plasma velocity must also be updated; since the freshly ionized plasma used to be gas then it is born with v ¼ 0 in the rest frame of the gas.The mean plasma velocity is now a combination of the plasma velocity and the gas velocity, weighted by density The updated densities, energy density, and velocities are the fed back into the fluid code which advances to the next timestep.From here, the process repeats itself.Figure 2 is a flow chart showing the incorporation of the AI subroutine to the MHD-gas momentum coupling code.
In addition, the simulations also include a coupling term in the momentum equations of the two fluids of the form 30 These model equations are solved by method of finite difference using a Lax-Wendroff scheme.

IV. RESULTS
An obvious source of relative velocity for AI is that of a homogeneous flow of neutral gas impinging upon a stationary, magnetized seed plasma.Now critical velocities are very large, of the order of several km/s, 16 but these are velocities frequently reached in laboratory experiments, 18  dwarf atmospheres 17 or observationally detected in the solar atmosphere.For example, Rubio 43 and Vitas et al. 44 show supersonic flows in the solar photosphere that are easily in the velocity range required.Bulk flows are not the only source of relative velocity in a gas-plasma mixture.Any acoustic mode naturally comes with ion and/or gas motion in the longitudinal direction and shear Alfv en waves propagating along field lines come with ion motion in the transverse direction without corresponding motion in the gas.Diver et al., 30 among others, have previously shown that waves in partially ionized plasmas have different phase speeds compared to either species treated separately.Further, due to the anisotropic speed of magnetoacoustic waves in the plasma, there can only be one direction where waves in both species in the presence of a magnetic field are perfectly matched.We also note that, due to the presence of a momentum coupling term, a wave that is initially driven purely in one species will result in complimentary wave motion in the other species.
Waves in a partially ionized plasma must therefore result in non-zero relative velocities.Waves of sufficient amplitude should contribute to Alfv en ionization.
All of the simulation results presented in this section are conducted on a 200 Â 200 numerical grid.The boundary is transparent where the ghost points are provided by a secondorder polynomial fit through the adjacent finite-difference points.The simulations also all share the same equilibrium conditions: The background is a partially ionized plasma with an ionization fraction of 0.5 meaning the same equilibrium densities in the plasma and neutral fluids (g ¼ 1).The magnetic field is of strength 1 orientated parallel to the horizontal, x, direction.The ratios of the plasma and gas sound speeds to the Alfv en speed were taken as r ¼ s ¼ 0.25.The mesh ratio of all the following simulations was set to k ¼ 0.25.
The inequalities given in Equations ( 4), (6), and (9) require the plasma to be in a certain parameter regime in order for Alfv en ionization to occur, we check that our choice of constants place our simulations in such a regime.
As previously stated, the inequality in Equation ( 4) is always satisfied as long as v 2 a =c 2 ( 1.This relationship is already assumed true in any MHD framework as such a constraint is necessary to remove the displacement current so as long as MHD is applicable Equation ( 4) is satisfied.
Equation ( 9) is essentially a restriction of plasma beta; the choice of ionization fraction implies a further choice of temperature that depends on the ionization and recombination rates for a chosen atomic species.For example, for hydrogen to be thermally ionized to an ionization fraction of 0.5 implies a temperature of T ¼ 10 5 K. 40 In combination with the values of r and s, this temperature implies a plasma beta of 0.075.At such a beta, the inequality in Equation ( 9) works out as 1 ) 10 À8 , therefore, the condition is satisfied for our choice of simulation parameters.
Substituting the values for constants into 6 gives the relationship that B ) 10 À27 n g ; this is fairly trivially satisfied for astrophysical plasmas.For example, the atmosphere of brown dwarfs and M-dwarfs have previously been shown to satisfy this condition for the majority of the extent of their atmosphere. 17The solar photosphere, with a number density of %10 23 m À3 only requires a field greater than %10 À4 T, a number easily exceeded by magnetic features such as magnetic bright points (>0.1 T) 45 and sunspots (0.2-0.37 T). 46

A. Plane waves
To investigate the AI mechanism as a result of waves, we simulate waves propagating through a 2-dimensional domain by solving our model equations, including the AI and momentum coupling terms, by method of finite difference.We first test a plane wave that propagates through the gasplasma mixture.Since motion in an acoustic wave is confined to the direction of propagation and only relative velocity in the perpendicular to the magnetic field direction can result in Alfv en ionization then only waves which have a component of the wave vector in the perpendicular direction are of interest.
Due to the linear nature of our numerical code, we are limited to small amplitude waves, those not exceeding 10% of the background density.This also limits the maximum fluid velocity (produced by the compressional wave) to a few % of both the Alfv en speed and the sound speed.
We first look at an infinitely plane wave that is launched in the perpendicular direction.The background is a partially ionized plasma with an ionization fraction of 0.5 meaning the same equilibrium densities in the plasma and neutral fluids.The magnetic field is of strength 1 purely in the horizontal, x direction.The driver is a plane parallel wave generated at the top of the domain, this driving term is sinusoidal with time dependent amplitude of the driver given by A t ¼ A max cos½2pxt, with A max ¼ 0.05, and acts on the gas density and y-direction gas velocity in order to produce an acoustic wave in the gas.The frequency, x, of the wave is set such that there 2 complete cycles before the driver is disabled.
The threshold for CIV is set at 0.01 in normalized energy units where the total energy density of a stationary fluid cell is unity.This low threshold is chosen to be consistent with the kinetic energy restrictions central to a linear code.These limitations mean that we must avoid large fractional changes in density.As such, this threshold is close to the kinetic energy of the maximum velocity of the driver (0.01 % c0.1 2 ).
The factor f, from Equation ( 13), is set to 0.01, meaning that 1% of the maximal AI occurs on each fluid timestep.This factor constrains the two relevant timescales at work, the fluid timescale (the sound speed multiplied by the mesh ratio) with the AI timescale-defined by 1/f as the characteristic timescale for a constant velocity flow, at v ¼ 2v c , to fully ionize the plasma.Both of these time scales depend on the ion-neutral collision frequency.
The driver runs for 200 timesteps (up to t ¼ 5) and the subroutine that controls the AI is enabled after 300 timesteps have elapsed (t ¼ 6.5).This delay is to ensure that the source wave has travelled some way into the computational domain before any ionization takes place, allowing AI to be more clearly observed.
Our simple energetics argument tells us that (for an ionization fraction of 0.5) if the kinetic energy of the relative velocity exceeds twice the energy at v c then that fluid cell contains a high enough energy density that it could just reach an ionization fraction of 1.0 before the relative velocity drops below the threshold and ionization ceases.This leaves a window of between v c < v ?< 2v c where the maximal possible ionization varies between 0 < _ q p ai dt < q p;0 .Obviously, the presence of frictional drag reduces the relative velocity without enhancing the ionization fraction which reduces the maximum ionization.

Alfv en ionization in the absence of momentum coupling
In order to maximize the contrast between new plasma that is created by AI and plasma which is just perturbed due to collisions (momentum coupling) from the motion of the gas then the coupling factor is first set to 0 (C ¼ 0).
Figure 3 shows time evolution of the driver wave in the gas.It starts at the top boundary and travels through the domain at the sound speed.t ¼ 6 corresponds to just before the AI subroutine is switched on.The wavefront remains flat throughout and drops in amplitude as kinetic energy is utilised in order to ionize.
Figure 4 shows temporally matched snapshots (to Fig. 3) of the plasma density.This maps out the location of new plasma that has been created by Alfv en ionization.We can see, at t ¼ 12, the first snapshot after AI starts, that there is a large burst of ionization across the full width of the wave.The intensity of the ionization drops off until the wave leaves the domain by which point it is much reduced in amplitude.
The total ionization over time is shown in Fig. 5, the y-axis is in units of the mean fluid density such that a perturbation of 1 would be the equivalent extra mass obtained by introducing a perturbation in plasma density of 1 to a single fluid cell, i.e., a single fluid cell doubling in density.To get the mean change in density, this can be divided by the size of the simulation (200 Â 200).In other words, a perturbation of 1 indicates the mean density of plasma, over the whole domain, has increased by 0.0025%, however, the usefulness of a mean measurement is limited as the majority of this excess density is confined to relatively small regions of the simulation as can be seen in Figures 4 and 8.
There is a large amount of plasma created at t ¼ 7, which is shortly after the AI subroutine is enabled.The rate then drops over time as both the energy available above the threshold and the number of fluid cells where the critical velocity threshold is exceeded shrinks over time.The drop in density resulting in a local minima at t ¼ 30 is due to a peak in density from a fast-mode magnetosonic wave (visible as the low amplitude wave in 4) leaving the domain and not due to any recombination, which is not modelled.Time resolution of these data is limited by full data output being restricted to every 40 timesteps (to reduce both disk space requirements and runtime).Once the mass leaving the domain is subtracted from the calculation, the gas and plasma mirror each other with the total (in green) remaining constant, because each unit of plasma created results in a corresponding unit of gas being removed.The ionization fraction levels off when there is no longer any part of the domain where relative velocity exceeds the CIV threshold.The formation of this new plasma creates a divergence in density which drives fluid motion of the plasma, this is seen as the slow spread out of the overdensity through the domain at later times.The plane parallel nature on the wave means that any gradient in pressure is only in the y (or cross field) direction which means that the magnetic field inhibits the dispersion of the overdense region.The slow diffusion of plasma from the overdense region means the locations where there has been a significant amount of AI remains visible in a map of plasma density (Figure 4) for a long time (many Alfv en wave crossing times), beyond the time when the critical velocity is no longer being reached in any fluid cell.
The effect of AI, on the plasma density, is local as well as global.This has the interesting consequence that the ionization fraction of a given fluid cell depends on the history of that cell not simply its current conditions.This is in contrast with a plasma where thermal ionization is the only relevant ionization mechanism.
Once all ionization has ceased then diffusion will smooth out these locally overdense regions.The plasma is on average at a higher ionization fraction than thermal equilibrium would dictate; recombination will act to reduce total ionization fraction in order to return to some new equilibrium.Due to the different timescales of ionization and recombination, the history of AI does persist for a while; the plasma has a limited memory of the fluid velocities that have passed through it.This theoretically allows a distant observer to infer the relative motion that a given fluid cell of plasma has been subject to by observing the plasma density.

Alfv en ionization in the presence of momentum coupling
The previous simulation was conducted with the momentum coupling factor, C, set to 0. Since both frictional drag and AI depend on ion-neutral collisions, though the relative magnitudes might be large, it is difficult to imagine one FIG.4. Perturbations around the mean plasma density with C ¼ 0, f ¼ 0.01.Red indicates the areas where the most plasma has been created.The maximum deviation from the mean is 0.001.FIG. 5.The evolution of the various densities in the computational domain, obtained by integration.Steeper gradients in the plasma density (blue) indicate more AI is taking place at that particular moment.The y-axis is the total perturbation relative to the equilibrium density in normalized units.
without the other.To see how the introduction of a non-zero degree of momentum coupling effects the amount of Alfv en ionization that takes place, we repeat the same simulation as Section IV A 1 with the only change coming from the degree of momentum coupling.
As expected, the presence of frictional drag lowers the total plasma produced, this is shown in Figure 6.There are a few specific features that require explanation, the first is the spike in plasma density between t ¼ 1 and 5.This spike is caused by the driver; the gas wave exchanges momentum with the initially stationary plasma which starts a flow.This flow reduces the pressure at the boundary edge which drags in new plasma from outside the boundary until the first 1/2 of the wavecycle completes and the direction of the flow reverses.The total amount of plasma that enters through the boundary is small (we can see that the total plasma density has approximately returned to 0 before the initial burst of ionization, at t ¼ 6.5, begins) compared to the total amount of Alfv en ionization that takes place.
When this same sympathetic magnetoacoustic wave exits the bottom boundary of the simulation (at t ¼ 22) there is a corresponding reduction in plasma density, this does not indicate a true reduction in plasma density such as would occur as a result of recombination.
The presence of the coupling reduces the sum total of ionization that occurs.Since there must be a scale length, below the ion-neutral collisional mean free path, where the ionized and neutral species are not collisionally linked; it is clear that even coupled fluids will have some degree of relative motion that can be exploited to produce new plasma.Partially ionized plasmas where the gas and plasma are closely coupled will always have less opportunity for ionization than one where more independent motion between the two fluids is possible.
Since we require the ions and neutrals to interact via collisions in order for AI to take place then we require the two fluids to be collisionally coupled.AI will always be hindered somewhat by the frictional drag between the two fluids.
A key parameter is the relative size of the momentum coupling factor C and the ionization efficiency f.If f ) C then, the plasmas and gas are not closely coupled at the Alfv en ionization scales.If f ( C then the plasma and gas are strongly coupled together at the scale of Alfv en ionization and little relative velocity is ever present at fluid scale lengths.

B. Cylindrically symmetric waves
Since Alfv en ionization only occurs in the perpendicular to magnetic field direction, we wish to examine the directional dependence of the AI term in our numerical model.
The form of the driver is a 2-dimensional Gaussian with a standard deviation of 3 code cells and a hard cut off at 10 code cells beyond which the amplitude is 0, the maximum amplitude of this Gaussian is initially A max ¼ 0.1 corresponding to 10% of the sound speed.This maximum amplitude is kept constant for one wave period and then reduced exponentially with a characteristic time equal to the period of the wave such that each subsequent cycle has its mean amplitude reduced by a factor of e.The driver acts on the gas velocity only; it first drives the gas radially outwards and then inwards with the same sinusoidal amplitude dependence as the plane wave case.
The result of this driver is a cylindrically symmetric acoustic wave moving outwards from the centre of the domain.The amplitude of the driving term decreases for each consecutive cycle.The total size of the region where the driver operates is 10% of the width of the computational area, i.e., a diameter of 20 cells.Both fluids have the same equilibrium density, the magnetic field is orientated in the horizontal direction and r ¼ s ¼ 0.25.The ionization is treated the same as the plane wave case; AI remains disabled until t ¼ 8.0 and then it is enabled with the energy threshold again set at 0.01 of the total background energy density in perpendicular-to-field direction, this threshold allows the chosen driver to produce some ionization without altering the plasma density beyond the linear limit.The momentum coupling constant C is set to 10 À5 , low but not completely disabled.This leaves us with the majority of the plasma motion being from AI rather than the momentum coupling term.As before t ¼ 1 corresponds to 40 timesteps, this is approximately the duration for a wave to cross 10 units of length.
The gas density is shown in Figure 7 and the plasma density in 8.At t ¼ 5 we see only a small signal in the plasma which is purely motion induced by coupling not by ionization.When the ionization is enabled there is a sudden flash of new plasma created.This new plasma fairly rapidly spreads out in the parallel direction, along the field, as this is the direction where the plasma motion is not inhibited by a v Â B force.This is unlike the plane wave case where the gradient in the parallel to field direction was 0, limiting dispersion.
We can visually see that new plasma continues to be created as time goes on as well as the hotspots where a large FIG. 6. Integrated plasma density over time for varying C. As the coupling increases the total amount of plasma produced decreases.The y-axis is in the same units as 5. concentration of ionization took place have their density diffuse over time.
The acoustic wave in the gas starts out cylindrically symmetric but we can see that over time an asymmetry develops.The asymmetry comes as a result of AI only being concerned with relative velocities orientated perpendicular to the magnetic field.This means that even though the wave initially has the same amplitude in all directions only the directions where v ? is large has the kinetic energy wave converted to internal energy of the plasma.
Figure 8 also shows the formation of various MHD wave modes.This is due to the combination of two effects: the motion of the gas produces sympathetic motion of the plasma via the momentum coupling term, as such there is a switch from pure MHD modes to hybrid magnetoacousticacoustic modes; the new plasma created by AI constitutes local increases in plasma density and pressure that begins to drive new waves and bulk flows in the plasma fluid.This second source of waves does not initially include corresponding motion in the gas but they will, again due to the coupling of the two fluid's momenta, induce motion in the gas over time.The anisotropic ionization pattern and the formation of these waves from both processes produce the intricate pattern of plasma density seen at later times.
It may be expected that the central region, where most of the initial ionization takes place would be where the plasma density would remain the highest; similar to the plane wave case where most of the plasma remained in approximately the region it was created.However, at late times a cold region (dark blue) is seen in the centre of the computational domain.This is where the plasma density is lower than the initial equilibrium density.This means the non-zero divergence caused by the increase in plasma density set up a net outward flow in the plasma that results in the centre ending up a plasma density below the equilibrium value.This same effect is not seen in the gas where the gas motion is purely one induced by waves.The density may rise and fall dynamically as the driving wave passes by but the fluid parameters of the gas return to the equilibrium fairly rapidly; there is no net transport of neutral material.
In Figure 9, we look at the total amount of plasma created by Alfv en ionization.Approximately, the same pattern is repeated here as in the perpendicular plane wave with ionization being rapid at first and falling off as the amplitude of the wave drops off.The falloff is more rapid due to the natural dispersion and corresponding amplitude drop of a cylindrical wave compared to a planar wave.When f is increased then the difference between the planar and cylindrical case is reduced since AI is able to extract more energy quicker, before the wave amplitude falls below the critical threshold.
Since the ionization rate in a given fluid cell depends on the excess velocity available above a threshold then we expect even a small difference in wave amplitude to be capable of large changes in ionization rate.When you are close to the critical velocity, AI is very sensitive to changes in velocity.A small increase in wave amplitude might FIG. 7. The perturbations around the mean gas density showing an initially cylindrically symmetric acoustic wave in the gas as a source of relative velocity for AI with C ¼ 10 À5 , f ¼ 0.01.The portion of the wave that propagates across the field loses significant amplitude due to the influence of AI whereas the fraction of the wave in the parallel-to-field direction is relatively unimpeded.move the partially ionized plasma from a regime where there is no, or a very small amount of, ionization taking place to one where ionization is complete.This behaviour is not unique to this numerical fluid approach but is also seen in experiment. 18

V. CONCLUSIONS
We have presented simulation results showing the effect of the addition of an Alfv en ionizing mechanism to the two fluid MHD-gas momentum coupling code described in Diver et al. 30 can have on both the local and global ionization fraction and the overall gas-plasma dynamics.
The addition of this term allows the kinetic energy of the fluid motion to be harnessed in order to ionize the neutrals.The presence of this mechanism depends on prerequisite criteria (given in Equations ( 4), (6), and ( 9)) being met and on the two fluids reaching a relative velocity greater than the critical ionization velocity (1) in the direction perpendicular to the background magnetic field.The excess energy above this threshold dictates the total ionization possible and a numerical factor f dictates the rate that this energy is converted into ionization.
By removing kinetic energy from the relative motion of the fluids, we influence their dynamics.The AI term is an anisotropic amplitude and frequency dependent damping term for wave motion that acts as an energy sink, alongside the momentum coupling term which merely redistributes energy between the species.In the case of AI, energy loss is only present when the velocity exceeds the CIV and only applies in the perpendicular to magnetic field direction.
The creation of plasma is shown to produce gradients in partial pressure that drive additional flows and waves in one or both of the species.In an idealised experiment such additional dynamics have diagnostic potential; inversion may reveal the presence of and degree of Alfv en ionization by the FIG. 9.The evolution of the plasma density for the cylindrical wave case.The ionization follows a similar pattern to the plane wave case, initially rapid before falling off.We have maintained the same y-axis as Figure 5, note the lower total degree of ionization due to the more rapid falloff in wave amplitude from the radial dispersion of the wavefront.magnitude of these effects and therefore reveal fluid velocities that may not be directly observed.
Previously, there existed a clear gap in the theoretical modelling of plasmas which enter the regime where Alfv en ionization could take place.Experiments have previously shown that in these conditions, there was a process taking place that inhibited relative velocities above the critical velocity, this is not behaviour that could be replicated by any previous fluid code.The introduction of an AI mechanism that acts as a damping force if and only if the relative velocity between the ordinarily independent species exceeds a set threshold allows simulations to reproduce this behaviour.The extent to which this limit is enforced smoothly varies as f is adjusted between 0 and 1.If f approaches unity then the limit is strictly enforced and no relative velocity above v c is allowed anywhere in the domain, if f is 0, such as in a collisionless plasma then the relative velocities are limitless such as in any previous MHD simulation.
The values of f and C (the AI efficiency factor and momentum coupling constants, respectively) strongly influence the degree of AI.Several microscopic processes are encapsulated in f: the formation of pockets of charge imbalance; the fraction of the electron distribution function that is capable of ionizing; and the probability that a high energy electron will ionize.These processes are not resolvable at the fluid scale.
We have distilled a series of microscopic processes, informed by high time-resolution PIC simulations, into a source term that depends only on fluid parameters.This approach naturally can be extended to incorporating other relevant phenomena that will extend the compass and impact of MHD models.The value of f used for these simulations (10 À2 ) is appropriate for capturing both the fluid dynamics and ionization properties.This approach could be helpful for contexts in which Alfv en ionization is already implicated in the physics but has yet to be modelled, such as the cool atmosphere of brown dwarfs or the gas-plasma mixing layer in magnetic fusion refuelling.
The MHD-gas code reported here is linear (therefore is restricted to small amplitude effects) but, in general, the density changes associated with AI are not necessarily small.The approach here can be extended to even up to the case of total ionization.The authors are preparing a fully non-linear gas-MHD plasma code that will allow significant evolution of the ionizaton fraction, beyond the linear limit.

FIG. 1 .
FIG. 1.(a) A stream of neutrals with velocity, v, encounters a stationary plasma with magnetic field perpendicular to the flow.(b) Ions are displaced via collisions with the neutrals, the electrons cannot follow and a resultant potential, /, is created.

FIG. 3 .
FIG. 3. The perturbations around the mean gas density showing an acoustic wave in the gas as a source of relative velocity for AI with C ¼ 0 and f ¼ 0.01.The maximum and minimum perturbations are 60.05.

FIG. 8 .
FIG. 8. Perturbations around the mean plasma density with C ¼ 10 À5 , f ¼ 0.01.Red indicates the areas where the most plasma has been created.Plasma is created by wave motion in the perpendicular-to-field direction but due to the v Â B force it diffuses along the field lines at approximately the sound speed.
theoretically expected in Brown FIG. 2. Flowchart showing the process for incorporation of an Alfv en ionization subroutine into a finite difference code.