Methods and Challenges of Including Dust Evolution in Magnetohydrodynamic Simulations

Dust plays many important roles in the interstellar medium, affecting the heating, cooling and chemistry in all the different phases of the medium. Up until recently though, if dust was included in simulations at all, it was typically assumed that it was tightly coupled to the gas. More recent work has allowed for the dust to decouple from the gas, though this brings with it challenges for including all the dust destruction processes, especially grain-grain collisions. We discuss the various methods of simulating dust evolution in the context of magneto-hydrodynamical simulations of the interstellar medium and their strengths and weaknesses. In particular we present a method for improving on the tight coupling approach that shows promise for simulating the evolution of small grains subject to shocks.


Introduction
Dust plays a crucial role in determining the thermal state of the interstellar medium (ISM) via a variety of mechanisms: • photoelectric heating • depletion of important coolants (e.g.ions of C and O) from the gas phase • cooling of hot gas by conversion of thermal energy into infrared radiation.
In addition dust in dense ISM is important for star and planet formation because of its role in catalyzing molecule formation.Molecules are important for cooling gas in those dense clouds that are optically thick to UV and optical radiation, which in turns allows the dense clouds to cool further and condense on their way to star and planet formation.By locking up metals, dust is also important for the transport of metals, both within galaxies and from galactic disks to halos.Dust formed in supernovae carry the metals synthesized in the progenitor and in the explosion into the ISM.Similarly, dust formed in the dense winds of AGB stars enriches the ISM with metals.Finally dust absorbs and scatters starlight, especially in the far UV, which is important for the ionization of low first ionization potential elements such as Mg, Fe and Si.

Why Simulating Dust is Hard
Given the importance of interstellar dust as outlined above, there has been increasing interest in including it in simulations of the evolution of the ISM of galaxies.However, there are a number of characteristics of dust that make it intrinsically difficult to include in numerical magnetohydrodynamical calculations.These include: The big peak in gyroradius occurs at the start of the shock cooling zone where the grains are betatron accelerated as the magnetic field is compressed along with the gas.We see the extremely large range of scales over which grain trajectories need to be followed.The grain charges range from ∼ −350 − +1600 with the highest charges in the hot shock region (on left) and the largest negative charge in cooling zone at ∼ 0.55 pc into the shock.
• dust grains are charged with the charging depending on a variety of factors including the electron temperature, the electron density, the interstellar radiation field and the speed of the dust relative to the surrounding plasma, • the charge-to-mass ratio is quite low for grains compared with ions, but also spans a wide range, leading to a wide range in gyroradii (see Figure 1), • for grains on the small end of the grain size distribution (which are the most numerous, though they carry less mass) their gyro-periods are short, demanding short timesteps to calculate their motions, • grain-grain interactions such as shattering and vaporization can be important for grain evolution and destruction and mainly involve collisions of small grains with large grains • dust instabilities, so called "resonant drag instabilites" can cause a high degree of dust clumping that doesn't coincide with gas clumping [1] These factors especially come into play in highly dynamical situations in which shocks and outflows occur.

Grain destruction in shocks
A primary motivation for including dust in simulations of the ISM is to study grain destruction, especially as it occurs when shocks are present.The mechanisms of grain destruction are sputtering of atoms off of grains by ion impacts and vaporization of grains in grain-grain collisions.Thermal sputtering occurs behind fast shocks that heat the gas to temperatures above about ∼ 10 5 K. So-called non-thermal or inertial sputtering occurs when grains are moving at a substantial speed relative to the gas.Vaporization in grain-grain collisions occurs when such collisions are at speeds high enough that one of the grains is completely destroyed [see 3].When a grain is moving with the gas and enters a fast shock, the gas is quickly slowed, in the shock frame, to ∼ 1 4 v shock .The grain is not slowed quickly and thus attains a gyro velocity of ∼ 3 4 v shock , the difference between the pre-shock speed in the frame of the shock and the post-shock speed of the grain relative to the gas (for high mach number shocks).In radiative shocks, some distance behind the shock, the gas cools and condenses and the component of the magnetic field perpendicular to the shock direction gets amplified (in addition to the factor of about four amplification at the shock front).This in turn leads to betatron acceleration of the grains and higher gas-grain relative velocities.In some cases grains can be reflected back upstream and undergo Fermi acceleration leading very large grain velocities, possibly helping with the injection problem for cosmic rays.Figure 2 shows some example trajectories calculated by [2] for various grain types and shock speeds.The trajectory in the upper right shows a grain that has been Fermi accelerated.shock with only grain-grain collisions or sputtering turned on.Note that the overall level of grain destruction for the two cases is not greatly different, despite the large differences in size distributions.
Shocks not only destroy grains but also change the grain size distribution.Observations of dust in the ISM, especially extinction of starlight, have lead to the inference that grains have a distribution of sizes ranging from ∼ 0.005 − 0.25 µm in the diffuse ISM, with roughly a power law shape in size, n(a) ∝ a −3.5 , where n(a) is the number of grains per grain size and a is the grain radius [e.g., 4].This implies that most of the surface area of the grains is in small grains and most of the mass is in large grains.Sputtering reduces the grain size by the same amount in radius, ∆a, per unit time regardless of grain size, so the fractional decrease in mass is larger for small grains.Thus sputtering tends to shift the grain size distribution to smaller sizes removing more grains from the small grain end.Grain-grain collisions leads to shattering of grains and vaporization.Shattering doesn't destroy grains, i.e. return grain material back to the gas phase, but does reshape the grain size distribution by shifting mass from large grains to small grains.Small grains are stopped more quickly by gas drag, which plays a part in promoting grain-grain collisions because of the relative velocity of the small and large grains.The differing effects of sputtering and grain-grain collisions are illustrated in Figure 3.

Methods of Calculating Grain Evolution in Hydro/MHD Simulations
The interest in understanding grain evolution in the ISM has lead to a variety of methods of including grains in hydrodynamical and MHD simulations.The three main types of calculation are: • tight coupling assumption: grains are assumed to move with the gas because of their small gyroradii and gas drag [e.g. 6, 7], • pressureless fluid method: dust is treated as a fluid (or several fluids) that can move relative to the gas and interacts with it via drag and possibly Lorentz forces [e.g. 8, 9], • free particle/kinetic approach: grains are treated as particles that are allowed to move relative to the gas subject to drag and (sometimes) Lorentz forces [e.g. 5, 10, see Figure 4 for an example].
Note that within these broad categories there are several variations among the approaches used.Each of the approaches listed above has its advantages and disadvantages.Though a thorough exploration of these would require more space than we have here, we briefly describe these pros and cons of each method.
For the tight coupling method the advantages are: 1) it is easy to include a range of grain sizes, for example using bins of grain sizes [3], and grain-grain interactions, 2) gyro speed can be included as an internal variable, 3) there are no stringent time constraints -one only needs to resolve the grain destruction and slowing timescales.The main disadvantage of the method is that large grains, i.e. grains that decouple from the gas on scales of the variations in gas properties, cannot be modeled accurately.Calculations of grain destruction using steady shocks [e.g.3], give the grains an initial gyro speed equal to 3  4 v shock , however in hydro/MHD calculations one doesn't know a priori when the grains will enter a shock.Thus coupled grains, which are essentially tracer particles with extra characteristics, do not naturally experience small scale decoupling at shock fronts and thus will not get the boost to their gyro speeds (carried as an internal variable) that they should, We discuss a method for overcoming this limitation below.The multi-fluid approach treats dust as a separate fluid that couples to the gas.This approach has the advantages that: 1) dust can slip relative to the gas, 2) using multiple grain size bins/fluids allows for inclusion of grain-grain interactions, 3) small scale decoupling occurs naturally.The main disadvantages are: 1) it assumes that all grains of a given size in the same region have the same velocity, 2) if grain gyro radii are large, grain fluids could overlap in space raising the question of how to define the fluid velocity, 3) very large decoupling, such as Fermi acceleration at a shock front does not seem possible to accommodate.More generally dust grains are not coupled tightly to each other; collisions between grains are infrequent so a fluid approximation is not justified.
The free particle/kinetic approach treats grains as free particles subject to interactions with the gas and, potentially, the magnetic field.The advantages are that: 1) large grains can be handled accurately without any assumptions about a uniform grain velocity, 2) particles are correctly accelerated at shock fronts, including any possible Fermi acceleration.The disadvantages mainly concern grain-grain interactions since those require knowing the characteristics of the other grains in the vicinity.Also, the accuracy of such calculations are dependent on the number of grains used and so could become computationally expensive for a large simulation.An additional issue is that in an MHD calculation, one needs to resolve the gyro periods which could be very short for small grains.

Improving the Tight Coupling Approximation
As detailed above, the tight coupling approach could, at least for small grains, be a good approximation that would allow inclusion of all the main processes that affect grain evolution, including grain-grain collisions.The lack of small-scale decoupling at the shock front, however is a stumbling block.In general the gyro speed of grains evolves via electro-magnetic acceleration and gas drag: where v ⊥ is the gyro speed (grain speed perpendicular to the magnetic field direction), B is the total magnetic field strength, F d is the drag force (∝ v 2 ⊥ ), and m gr is the grain mass.An implicit assumption is that except for gyrating around the magnetic field, the grain does not move relative to the gas.Said another way, the guiding center of the grain motion moves with the gas.One thing to note is that the betatron acceleration from the first term on the right is only effective once the grain has a significant gyro velocity.We note that we are only concerned here with the velocity perpendicular to the magnetic field since, motions along the field cannot increase the gyro velocity.Thus this method is incapable of following grain motion along the magnetic field that is decoupled from the gas motion.The drag force then should only use the gyro velocity of the grains.
It is clear that grains should experience an effective boost to their gyro velocity any time that the gas speed perpendicular to the magnetic field changes quickly over distances smaller the grain's gyro radius.Using this insight and via experimentation we have found that an effective method is to 1) add ∆v gyro = ∆v gas to the v gyro iff ∆v gas > 0.1c sound , 2) when this boost is applied, betatron and drag terms are turned off.(c sound is the local sound speed.)Note that ∆v gas here is only the component of the gas speed perpendicular to the magnetic field.Thus during time steps that the criterion is met the grains increase their gyro speeds by an amount ∆v gas , but this only occurs near shocks.If shocks were infinitely sharp in the simulation, one might want the factor of 0.1 above to be closer to 0.75, however with the shock spread over a few zones, the factor of 0.1 has worked well so far.
Our development of this approach was done using the MHD code FLASH, version 4.6.2 with some additional physics that we have added as enhancements to the code.We add optically thin radiative cooling and heating via a table look-up.With regards to dust processes, we have implemented, for tightly coupled grains, time advancement of position and velocity using mapping of gas properties from the grid to the grain positions, which is already included in the code, but with the addition of a calculation of the grain charge using a table look-up scheme and the calculation of the gyro velocity in the way discussed above.We generated the grain charge table using code we have written which incorporates code generously provided to us by J. Weingartner [see 11], extended to include relative gas-grain motions.Gas (collisional) drag is calculated and leads to slowing of the gyro speed.Thermo-kinetic sputtering of the grain is also included.We have implemented both the tight coupling scheme and a kinetic scheme in which we calculate the grain trajectories via a Boris pusher method.The shock profile is set up initially via a steady shock calculation.Our results presented here are for a 200 km s −1 shock, which we found to be stable and not subject to the cooling overstability that affects somewhat slower shocks [12].The shock profile is shown in Figure 5.The grains are spread over several zones upstream of the shock and moving with the gas.The calculation is in 1D, though for the Figure 6.Comparison of the grain gyro speed (top) and grain mass (bottom) for our tight coupling and kinetic grain calculations.As can be seen there is excellent agreement, with the gyro velocity boost method performing very well in giving the tight coupling grains the appropriate gyro velocity after entering the shock.free particle scheme, the grains are allowed to move in 3D.Since the magnetic field is assumed to be in the z direction and the gas motion is in the x direction, the grain motion is only in x and y.
Figure 6 shows a comparison of our results for the grain gyro speed and grain mass for coupled and kinetic assumptions for 0.1 µm silicate grains.As can be seen, the agreement is excellent for both the gyro speed and the grain mass.Very little tuning was done to achieve these results, which leads us to believe that they are robust.Further testing will be necessary to confirm that, especially in 2D and 3D and for oblique shocks.

Summary
Including dust grains in magneto-hydrodynamical calculations of the evolution of the ISM is complicated by the large range of size scales involved and because of their self-interaction via grain-grain collisions.Various methods for including grains have their strengths and weaknesses.Some methods are fairly accurate for small grains but not for large grains and vice versa.A kinetic approach, which could be accurate for large grains is expensive for small grains because of their small gyro radii.
A tight coupling approach is appropriate for small grains and can allow for relatively easy inclusion of grain-grain collisions, however, for use in an MHD code, a mechanism is needed to allow for the boost to the gyro speed imparted by small scale decoupling of the grains from the gas at shock fronts.We have found a method that allows for the correct behavior at a shock front without explicit shock tracking and depends only on the change in gas speed and local sound speed.For a radiative shock that is stable, we have shown that the grain gyro speed and mass match closely those found for a kinetic grain treatment.
The ultimate goal of this work is to find a method that is accurate and not too expensive for calculating the evolution of grains over the whole observed size range and including all the physical processes expected to be important in such environments.One possible avenue is to include small grains with a tight coupling approximation and large grains with a kinetic treatment in the same simulation.Then the erosion of the large grains via collisions with the small grains could be calculated.Such calculations could improve our understanding of grain destruction in clouds that are overrun by supernova remnant shocks, which is a key component of dust evolution in the ISM.

Figure 1 .
Figure1.Comparison of the gyroradii for 0.1 µm and 0.25 µm silicate grains as they move through a steady 200 km s −1 radiative shock.The big peak in gyroradius occurs at the start of the shock cooling zone where the grains are betatron accelerated as the magnetic field is compressed along with the gas.We see the extremely large range of scales over which grain trajectories need to be followed.The grain charges range from ∼ −350 − +1600 with the highest charges in the hot shock region (on left) and the largest negative charge in cooling zone at ∼ 0.55 pc into the shock.

Figure 2 .
Figure 2. Grain trajectories behind steady radiative shocks for various shock speeds and grain types [from 2].A range of behaviours is seen including relatively tight coupling (upper left), reflection upstream leading to Fermi acceleration (upper right), and moderate decoupling (lower right) and finally almost full decoupling of a very large grain that passes through the shock only slightly deflected.Note that the x-axis limits for the inset plots are the same as for the corresponding larger plots.

Figure 3 .
Figure 3.Comparison of the effects of sputtering and grain-grain collisions on the grain size distribution.The size distributions are for grains that have passed through a 200 km s −1shock with only grain-grain collisions or sputtering turned on.Note that the overall level of grain destruction for the two cases is not greatly different, despite the large differences in size distributions.

Figure 4 .
Figure 4. Evolution of dust in an HD simulation of the reverse shock overrunning an ejecta clump in an expanding supernova remnant [from 5].The green dots are representative dust grains.There is no magnetic field in the simulation.Note that the grains, which start out close together, are widely scattered by the last frame.In this simulation grains are formed in dense ejecta clumps and escape the remnant into the ISM.

Figure 5 .
Figure 5. Shock profile used for testing grain evolution code.The calculation takes place in the shock frame with the upstream gas flowing in from the left and the shock at 0.5 pc.The gas flows off the grid on the right with outflow boundary conditions.