Modelling grain-size distributions in C-type shocks using a discrete power-law model

In this paper we discuss the implementation of a discrete, piecewise power-law grain-size distribution method into a numerical multifluid MHD code as described in Sumpter (2020). Such a description allows to capture the full size range of dust grains and their dynamical effects. The only assumptions are that grains within a single discrete bin have the same velocity and charge. We test the implementation by modelling plane-parallel C-type shocks and compare the results with shock models of multispecies grain models. We find that both the discrete and multispecies grain models converge to the same shock profile. However, the convergence for the discrete models is faster than for the multispecies grain models. For the pure advection models a single discrete bin is sufficient, while the multispecies grain models need a minimum of 8 grain species. When including grain sputtering the necessary number of discrete bins increases to 4, as the grain distribution cannot be described by a single power-law as in the advection models. The multispecies grain models still need more grain species to model the distribution, but the number does not increase compared to the pure advection models. Our results show that modelling the grain distribution function using a discrete distribution reduces the computational cost needed to capture the grain physics significantly.


Introduction
Although the typical grain-size distribution in the Interstellar Medium (ISM) follows the Mathis et al. (1977) (MRN) distribution, local variations are expected as the grain-size distribution is subjected to different grain processes such as dust production (Maercker et al., 2018), grain growth by coagulation and mantle accretion (Jones and Williams, 1985;Liffman and Clayton, 1989;Ossenkopf, 1993;Inoue, 2003;Ormel et al., 2009;Asano et al., 2013;Ysard et al., 2016;Jones et al., 2017), and destruction processes like sputtering, shattering and vaporisation (Tielens et al., 1994;Jones et al., 1996;Flower and Pineau des Forêts, 2003;Hirashita and Yan, 2009;Guillet et al., 2007Guillet et al., , 2009Guillet et al., , 2011Anderl et al., 2013). These grain processes do not only affect the overall grain distribution, but also have significant effects on the dynamics of the ISM. In the outflows of Young Stellar Objects (YSOs), dust grains are important charge and current carriers and therefore determine the structure of C-type shocks (Van Loo et al., 2009). Furthermore, grain-processing leads to observational signatures. In typical ISM conditions silicon is adsorbed onto dust grains. However, gas phase SiO is detected in the clumpy structure of YSO outflows due to the shock-induced sputtering releasing silicon into the gas phase (e.g. Martin-Pintado et al., 1992;Mikami et al., 1992). Thus, it is essential to accurately follow the evolution of the dust grain distribution to model both the dynamics and emission signatures.
A number of different approaches have been used to model the dynamics and processing of dust grains in Ctype shocks. In the simplest method, i.e. the multispecies grain model, only a few dust species with specified radii, i.e. typically, one or two grain species representing small and/or large dust grains, are evolved (e.g. Draine et al., 1983;Van Loo et al., 2009. This approach captures not only the dynamical importance of dust grains to the shock structure, but also the sputtering process of the grains. However, the limited number of dust species fails to represent the grain-size distribution adequately and is furthermore restricted in use as it cannot model grain-grain processes such as shattering and vaporisation (Sumpter and Van Loo, 2020). A more rigorous approach uses a discrete distribution function in which the grain-size distribution is modelled using different size bins spanning the full radial range (e.g. Jones et al., 1996;Guillet et al., 2007Guillet et al., , 2009Guillet et al., , 2011Anderl et al., 2013). While both sputtering and grain-grain processing are followed to provide a more realistic shock and emission profile, these models only follow the total dust mass and approximate the dust number density. McKinnon et al. (2018) improved this method by using a linear piece-wise discretisation in a grain size bin to conserve both mass and number density and showed it to be second-order accurate in number of size bins. However, this still requires 50-100 bins to achieve an accuracy of the order of 1%. A further improvement was suggested by Sumpter and Van Loo (2020, hereafter Paper I): they show that a power-law discretisation of the distribution function preserves the global distribution properties to a high degree even with small bin numbers.
In this paper we will implement this discrete power-law method in a multifluid magnetohydrodynamics (MHD) code (Sect. 2). Then, in Sect. 3, we will compare the profiles of oblique C-type shocks obtained with both the powerlaw discretisation and multispecies approaches, while, in Sect. 4, we will present the results for the sputtering of grain cores with the aim of determining the number of size bins required to accurately model this destructive grain sputtering process. Finally, Sect. 5 contains a discussion of the results and conclusions.

Dynamics of discrete grain distributions
To model a weakly-ionised dusty plasma, we use the multifluid MHD code MG with the chemical network described in Van Loo et al. (2009. The numerical scheme underlying the code solves the continuity, momentum and energy equation for the neutral particles. In the limit of small mass densities for the charged fluids, their inertia can be neglected and, also, the charged fluids are in thermal balance as their heat capacities are small. Therefore we can restrict ourselves to solving the continuity and reduced momentum equations for all charged fluids and the reduced energy equations for ions and electrons only (as the dust grains can be treated as a pressureless fluid). Furthermore, the neutral and charged fluids are coupled by mass, momentum and energy transfer coefficients. These equations are solved using a secondorder hydrodynamic Godunov solver for the neutral fluid equations. The charged fluid densities are calculated using an explicit upwind approximation to the mass conservation equation, while the velocities and temperatures are calculated iteratively from the reduced momentum and energy equations. The magnetic field is calculated using the induction equation with the resistivities obtained from Ohm's law. The magnetic field is advanced explicitly, even though this implies a restriction on the stable time step at high numerical resolution (Falle, 2003).
However, we make three modifications. While Van Loo et al. (2013) modelled grain sputtering and thus grain mass loss, the mass loss was not taken into account when determining the grain dynamics. Our first modification is to relax this assumption. The second change concerns the dust grain's charge. While in previous models the grains are always assumed to be negatively charged, this is also relaxed and, in specific conditions, grains can attain a positive charge. The final change relates to the calculation of the velocities of the charged particles from the reduced momentum equations and provides a more robust derivation of these. The latter two modifications are described in Appendix A and Appendix B.
When incorporating the grain-size distribution in the multifluid MHD scheme the equations governing the (pressureless) grain fluids need to be modified. As the description of the implementation is detailed in Paper I we only give a short summary here. We assume that the grain number density distribution in bin i, given by the size interval [a i , a i+1 ], has a power-law shape, i.e. ∂n(a, t) ∂a where A i is the power-law coefficient and α i the power-law index which both implicitly depend on time. Then, by integrating the continuity equations and reduced momentum equation over the grain number density distribution in bin i, the governing equations for bin i become (4) where n i is the total number density in bin i, ρ i the total mass density, and Z i e the average grain charge. Furthermore, is the mean specific collision coefficient between neutrals and grains in bin i (and is determined by integrating the collision coefficient K gn ). Here, the average value of a variable f is given by The only assumption made in the derivation of these equations is that all the grains within a bin have the same velocity (v i ). While this is not necessarily valid, it is expected to only have a minor effect. The grain velocity depends effectively on the grain radius through the Hall parameter, which is the ratio of the gas-grain collision frequency to the gyrofrequency, β = ZeB/(mρ n K gn ) ∝ a −1 . For small Hall parameters, i.e. when |β| < 0.1, the grains move with the neutrals, while, for |β| > 2 they move with the electrons and ions. Thus, there is only a small range of β values, or grain radii, for which grains have a velocity in between these limits. S i,sputt and S i,sputt represent the sputtering losses in bin i and are calculated as described in Paper I. For this, we specifically need to calculate the rate of change of the grain radius, i.e.ȧ, in each bin due to sputtering. The sputtering yields due to H 2 , CO, O, Mg, H 2 O and SiO projectiles are calculated using the same expressions as in Van Loo et al. (2013). It is important to realise thatȧ is not explicitly dependent on the grain size, only implicitly via the grain-neutral relative velocity. Given that the velocity is assumed to be the same for all grains in a bin, this relative velocity, and thus alsoȧ is also constant in each bin (which is implicitly assumed in Paper I).

Initial conditions
We follow Van Loo et al. (2009) andVan Loo et al. (2013) for the initial set-up of the shock models. This assumes that the initial flow profile corresponds to that of a discontinuity similar to a J-type shock. The upstream magnetic field is B 0 = 1µG(n H /1 cm −3 ) 1/2 , where n H is the upstream value of the number density of hydrogen nuclei. A shock propagates with a speed v s at an angle θ = 45 • with respect to the upstream magnetic field. We will consider two different shock set-ups corresponding to models with or without grain sputtering. For the advection only models (no sputtering), the pre-shock density is either n H = 10 4 cm −3 or n H = 10 6 cm −3 with a shock speed v s = 25 km s −1 , while, for the sputtering models, the shock propagates with a higher speed of 40 km s −1 . We adopt this larger shock speed for sputtering as 25 kms −1 is the threshold speed for core sputtering and would have little effect on the dust distribution. Furthermore we adopt the same fractional abundances for the upstream chemical species, that is the upstream gas phase fractional abundances of O, Mg, H 2 O and SiO are 4.25 × 10 −4 , 10 −7 , 0 and 0, respectively. The fractional abundance of CO we keep constant at 5 × 10 −5 throughout the numerical domain.
The dust grains are modelled using an initial powerlaw size distribution of grain cores with an index of -3.5 (Mathis et al., 1977, referred to as the MRN distribution function) and set the lower and upper radial limit of the distribution to a min = 5×10 −7 cm and a max = 3×10 −5 cm, respectively. We do not include grain mantles in this paper as we consider a shock speed above ≈ 20km s −1 for which grain mantles are completely, but also rapidly, eroded (Guillet et al., 2011;Van Loo et al., 2013). In the upstream region, one percent of the mass is contained in grains. The full radial range of the initial distribution is divided up logarithmically with a spacing determined by where N is the number of bins. Then the edges of bin i are effectively a i = a min e i∆a and a i+1 = a min e (i+1)∆a . The number density n i in each bin i can then be calculated as and the mass density ρ i as where ρ g is the density of the grain material. For models with the multispecies grains, an appropriate radiusā i needs to chosen for each bin. As the average mass per grain is given by ρ i /n i , it is clear thatā i = a 3 i 1/3 .

Advection of the grain distribution
While Paper I studied a numerical approach for accurately evolving a dust grain-size distribution undergoing grain processing, it only suggested an implementation of the technique in a (magneto)-hydrodynamical code. Therefore, it is yet unclear whether a discrete distribution function can adequately be used for modelling the dynamics. Here we compare the shock structures using both single-species and discrete distributions with the number of species (or bins) ranging from 1 to 16. Figures 1 and 2 show characteristic measures for the shock structure: the shock width measured as the part of the shock with a neutral temperatures above 100 K (as in Anderl et al., 2013) and the maximum neutral temperature as hotter temperatures are achieved in thinner shocks with a narrower radiating layer. We find that the shock is thinnest (2-2.5 times) and hottest (25%) for the model with one grain species. As more multispecies grains are included, the shock width and maximum temperature seem to converge. This convergence is more apparent if we compare the results with the ones from the discrete distribution function. Both models tend to the same values. However, contrary to the multispecies grain models, the shocks widths and maximum temperatures for the discrete distribution models have quasi-constant values even for a low number of bins. For the n H = 10 4 cm −3 runs the relative variation in the shock width is less than 1%, while it is slightly higher for the n H = 10 6 cm −3 model at 6%, compared to 52% and 56% for the multispecies grain models respectively.
To understand this behaviour, we need to understand the physical process governing the shock structure. In particular, the shock width is set by balancing the Lorentz force with the drag forces. From summing the reduced momentum equations (Eq. 4) for the charged particles together, we find Although dust grains are not necessarily the dominant charge carriers (they are for n H = 10 6 cm −3 , but not for n H = 10 4 cm −3 ), they do contribute significantly to the drag force. In regions with low ionisation fraction χ, i.e. χ < 3 × 10 −9 (v s /km s −1 ), as in these models, the drag force of the neutral-grain collisions dominates the ionneutral drag force (Van Loo et al., 2009). As the neutralgrain drag force is proportional to ρ g K gn , where ρ g is the grain mass density, the smaller grains in a MRN size distribution contribute the most to the drag. In the discrete distribution models, the collision coefficient is averaged over all of the grain radii within a bin (see Eq. 5) so that the evaluation of the drag force only depends on the number of bins needed to adequately describe the distribution function. For n H = 10 4 cm −3 all grains move with the electrons and ions as seen in Fig. 3. This means that the dust distribution keeps the same initial power-law distribution as the grains move from upstream to downstream and can, thus, be described by single discrete bin. However, the n H = 10 6 cm −3 models show a more complicated situation. The higher neutral density implies an increased collision frequency and the larger grains become decoupled from the electrons and ions (see Fig. 4) and even move along with the neutrals for a short distance. This is because, for n H = 10 6 cm −3 , the grains become the dominant charge carriers depleting the available free electrons that attach to the grains causing the largest grains to become quasi-neutral. Consequently, their Hall parameter is smaller than unity and the grains start moving with the neutrals so that these grains do not contribute to the drag force. The distribution function then needs to be modelled using more than 1 discrete bin. However, as the contribution of the large grains to the dust drag is limited, the overestimation of the drag force is small as can be seen in Fig. 2.
For the single-size models the collision coefficient is not evaluated across a range of grain radii, but just at a i . Remember that this value is determined from the initial distribution function (see Sect. 2). When considering only one grain species or discrete bin, it is possible to compare the drag force exerted by the dust, i.e. the ratio of the single size to discrete bin model is given by ρ g K g n/n g K * gn = a 3 2/3 / a 2 where the average is determined across the full grain size range [a min , a max ]. For an initial MRN distribution this ratio is then ≈ 2.4 which also is the relative difference observed in the shock widths between the single size and discrete bin models in Figs. 1 and 2. Including more grain species will reduce this overestimation in the drag force until convergence. We find that at least 8 to 16 grains species are required for the results to converge.
Both the discrete models and the multispecies grain models converge to the same values giving confidence that the discrete distribution models reproduce the correct shock structures. The importance of using a discrete distribution of grain sizes is demonstrated, with multispecies fluids only accurately achieving the same shock profile as the distribution for upwards of about 8 fluids. While discrete models with a single bin produce good results, it is useful to use a minimum of two bins to model the different dynamics of small and large grains. This becomes especially important when including grain processing physics such as sputtering in the next section.

Sputtering of the grain distribution
Here we explore the effect of sputtering as the initial MRN distribution is advected through the shock front. In the previous section we found that, for advection-only models, only two discrete bins are necessary to correctly simulate the shock dynamics. However, grain sputtering depends on the relative velocity between the grains and the impacting particles which can either be ions or neutrals. Therefore it is not only important to describe the grain size distribution but also the velocity distribution. As the grain velocity in a size bin is assumed to be constant, the velocity distribution of the grains is better described by more size bins as will be the sputtering of those grains. Figures 5 and 6 show the shock properties for a 40 kms −1 shock propagating in a n H = 10 4 cm −3 and 10 6 cm −3 medium as function of the number of discrete bins or multispecies while including grain sputtering. Similarly to the advection models ( Figs. 1 and 2) the shock width and the maximum neutral temperature converge as more discrete bins or multispecies are used to model the grain distribution. This is more so for the multispecies models than for the discrete bin models. For n H = 10 4 cm −3 , the multispecies models have a relative variation in the shock width from 19% for N = 2 and dropping below 5% for N = 8, while the discrete bin models are all within 2%.     For n H = 10 6 cm −3 , the relative differences are larger as in the pure-advection models: 34% and 12% for the multispecies and discrete distribution models respectively when N = 2, but fall below 5% for N = 4 for the discrete models and N = 8 for the multispecies models. Similarly the maximum neutral temperature exhibits the same convergent behaviour although the variation is smaller, i.e. all relative differences are smaller than 5% except for the multispecies grains for N = 2 when the variation is 6% and 10% for n H = 10 4 cm −3 and 10 6 cm −3 respectively. Note that the maximum neutral temperature for the n H = 10 6 cm −3 models is above 4000 K which is when H 2 dissociation becomes important (Draine et al., 1983). Then molecular cooling is significantly reduced and it is unlikely that a Ctype shock can exist in these conditions. However, as we do not include H 2 dissociation, the transition from C-type to J-type shock does not naturally occur in our models. When including sputtering, more discrete bins are needed before convergence of the shock properties. Furthermore, the multispecies models still need to include more multispecies grains than bins needed in the discrete distribution models, i.e. N = 8 compared to N = 4. When including sputtering the shocks widths are larger than for the advection-only models. This can be readily explained by the neutral-grain drag force: as the grains experience sputtering the grain distribution alters. At the small grain radii the distribution drops below the typical MRN distribution function (see e.g. Fig. 7). As the smallest grains contribute most to the neutral-grain drag, the total drag force decreases. To maintain the force balance in Eq. 10 the shock width needs to increase (as j × B = (∇ × B) × B ∝ B 2 /L where L is the shock width). This can be also inferred from Fig. 8 which shows the normal velocity for the ion and neutral species for both advection and sputtering models. While the ion velocity does not change between the models, the neutral velocity shows a slower deceleration in the shock frame. The deviation starts when the relative velocity between the ions and neutrals gets above 25 kms −1 . This corresponds to the threshold for grain core sputtering (most grains move with the ions and electrons). As the distribution changes due to sputtering and the grain-neutral drag reduces, the neutrals are not decelerated as quickly as in the advection only model.
Furthermore the neutral-grain drag force also explains the relative difference of 1-2.5% in the shock width between the discrete distribution and multispecies grain models contrary to the convergence seen in the advection only models. We see in Fig. 7 that the multispecies distribution function has the same overall shape as the discrete powerlaw distribution function, but lies a bit higher. This is because, in our models, grains are assumed to be completely destroyed if their radius drops below a min = 5 × 10 −7 . As the discrete distribution describes grains with a power-law size-distribution in each bin, dust grains with radii close to a min will lose enough mass due to sputtering to reduce their radius below a min . In the multispecies grain model with N = 8 grain species, the smallest grains, with radius a 0 , do not experience enough sputtering for their radius to be reduced below a min . Hence, the number of grains remains constant in the multispecies models and decreases in the discrete model. Note that, if N becomes large enough in the multispecies models, the smallest grains will be also removed from the distribution. This in return changes the amount of drag experienced between neutral and charged particles and changes the shock width.
While we already discussed the fact that the dust grains are sputtered and that the grain-size distribution alters, we have not discussed how the sputtering process affects dif-ferent grain sizes. This differs between the low and high density models. For the low density models, the grains contribute to the drag force, but do not carry much of the charge. As a result, all dust grains have a large Hall parameter and move with the other charged particles, i.e. the electrons and ions. As the relative velocity difference with the neutrals is then the same for all grains, the rate of change of the grain radius da/dt is independent of grain radius. However, although all grains experience the same rate of change ∆a, its effect is more readily seen at small grain radii as it is a larger fraction of their initial size. Figure 7 shows this effect by normalising the distribution with the initial MRN distribution. While the large grains still follow the initial MRN distribution, the distribution at small grain radii has a smaller power-law index. As seen before, both the multispecies models and discrete bin models produce similar results although the multispecies models lie slightly higher. Sputtering actually removes about 13% of the total number of grains in the discrete model, while the multispecies models do not lose any. Thus the difference is mainly due to the normalisation with different total grain numbers. However, we note that the mass lost in both the single-size and discrete distribution models is similar at about 4% as the small grains that are lost do not contribute much to the total grain mass. For the higher density models, the grains do carry a significant amount of the total charge. In fact, they dominate the free electrons in a large region of the shock front. Due to the lack of free electrons the large grains are close to zero or even positive charge and, thus, have a small Hall parameter. The large grains then move predominantly with the neutrals. As the dominant sputtering projectile is H 2 , the sputtering yield from the larger grains is non-existent. Only the small grains (up to 6×10 −6 cm) contribute to the sputtering and release of grain material such as silicon. However, the resulting downstream grain distribution function has a similar shape as the one for n H = 10 4 cm −3 (Fig. 7) with the large grains following the MRN distribution and a turnover to lower power-law indices at small grain radii. So, although sputtering acts differently at the two densities, this cannot be inferred from the downstream distribution function.
As the dust grains loose mass due to sputtering, silicon returns to the gas phase as SiO. Figure 9 shows the post-shock abundance of SiO for both background densities. Again we notice a convergence of the values as the number of bins or species increases. However, the relative difference is quite small: the only relative difference above 5% is when using a multispecies model with N = 2 and even that is only 8%. All other models have relative differences below 2%. Note that the shape of the SiO abundance as function of number of bins or fluids for n H = 10 4 cm −3 follows that of the shock width, but does not for the n H = 10 6 cm 3 models. This is because, as discussed above, in the lower density model, all grains move together with the other charged particles and the sputtering region within the shock front is a fraction of the shock width. On the other hand, for the higher density, the larger grains decouple from the magnetic field and move with the neutrals. With lower number of bins or fluids, it is not possible to correctly capture the velocity distribution across the entire grain size distribution and the sputtering yield is either over-or underestimated and is no longer a constant fraction of the shock width as in the lower density models.

Discussion and Conclusions
In this paper we implemented a discrete grain-size distribution method into a numerical multifluid MHD code as described in Paper I. Such a description allows to capture the full size range of dust grains and their dynamical effects. The only assumptions made are that grains within a single discrete bin have the same velocity and charge. We tested the implementation by modelling plane-parallel Ctype shocks and compared the results with shock models of multispecies grain models.
For models without grain processing the advection of the dust distribution function is to very good accuracy modelled by a single discrete size bin even for models where the large and small grains have a different velocity. The error in the characteristic measures is 1% for n H = 10 4 cm −3 and only 6% for n H = 10 6 cm −3 . This is not surprising in the low density case as all grains have large Hall parameters and move with the same velocity as the ions and electrons throughout the shock front. However, for the high density case, the large grains move not with the electrons and ions. Then a single discrete bin still produces very accurate shock profiles because the shock structure of C-type shocks is dominated by the dust-grain drag exerted by the small grains. Furthermore, the grain-size distribution does not vary much from the initial distribution as it is advected through the shock front. Contrary to the discrete bin models, shock profiles produced with multispecies grains cannot reproduce the shock structure with a single grain species. They actually produce much hotter and thinner shock structures. The models require the use of minimal 8 grain species for the shock profile to converge to the same result as the discrete bin models.
Following the pure-advection models we considered the effect of grain processing. In this paper we only take into account grain sputtering as the discrete bin models are directly comparable with multispecies grain models. We use the method described in Paper I for number-conserving processes. Similarly to the pure-advection models, the discrete bin models converge quicker than the multispecies grain models. As the grain sputtering changes the size distribution, especially at the smaller size range, the distribution no longer can be described with a single bin, but has to be increased to about a minimum of 4 bins. The main reason for this is that sputtering changes the grain size distribution throughout the shock as material is removed from the grains. However, this is still less than for the multispecies models which need a minimum of 8 species. The converged results of the multispecies and discrete bin models are not identical but similar to a relative error of a few percent. This is because the sputtering process is treated slightly different in the two models, i.e. in the discrete bin models grains with radii below a min are removed from distribution.
We find that sputtering increases the shock widths of the C-type shocks as this grain process removes a fraction of the small grains and, thus, reduces the dust-grain drag balancing the Lorentz force. Furthermore, as material is removed from the grains, SiO is released into the gas phase. While sputtering has a large effect on the shock properties, especially for the multispecies grain models, the effect on the SiO abundance is much less. The relative error when modelling the distribution with N = 2 discrete bins or multispecies grains is less than 10%. Furthermore the fraction of Si removed from the grains is in agreement with the results of Van Loo et al. (2013) whose models did not include the mass loss of the grains in determining the shock structure. Because only a small part of the shock front contributes to sputtering and thus to releasing SiO, the dynamical effects do not have a significant impact on the sputtering result.
Our results show that the implementation of the dust distribution function within a multifluid MHD code modelling C-type shocks is successful as the piecewise powerlaw and multispecies grain models converge. The only assumption that remains is that the velocity of the grains in each bin is equal. Therefore this method can be readily applied to other problems in which single or multispecies dust models have been used such as in e.g. dust transport and evolution in galaxies and protoplanetary discs. Most studies in these research areas either follow the dust distribution evolution (Asano et al., 2013;Aoyama et al., 2018;Granato et al., 2021) or the transport of a single dust species (Kanagawa et al., 2018;Hu and Bai, 2021).
However it is obvious that the dust distribution affects the dust dynamics and vice versa. Recent studies of planets in protoplanetary discs (e.g. Drażkowska et al., 2019;Karlin et al., submitted) reveal the interplay of multiple grain species on each others dynamics. Also, simulations of the streaming instability which plays a role in the early stages of planet formation show that, when multiple dust species are considered, the resulting instability growth rate depends strongly on the number of dust species considered (Krapp et al., 2019). Zhu and Yang (2021) show that, in the regime of slow growth rate, the number of species required to achieve convergence of the growth rate is very large (of the order of 1000 species). It is clear that the discrete power-law method would be highly beneficial in this case.
The only other method incorporating dust dynamics and grain size evolution in hydrodynamical simulations is by McKinnon et al. (2018) who use a piecewise linear method in the moving-mesh code AREPO. They demonstrate its applicability on determining the contribution of different grain processes, such as sputtering, shattering, coagulation and dust growth on setting the dust grain size distribution. However, they do not investigate the resolution of the dust distribution on the results. In Paper I we showed that the piecewise power-law method needs significantly less bins than the piecewise linear method to achieve the same accuracy. With a minimum of 4 discrete bins for just modelling grain sputtering using the power-law method, this means that it is unclear whether the piece-wise linear method is useful in practise. In a subsequent paper we will investigate the additional effect of grain shattering and vaporisation on the dust-distribution throughout a C-type shock and compare it to the piecewise linear method. n e . If we ignore the relative velocity of the ions and grains or v gi = 0, we find that this corresponds to When this situation occurs, the small grains in the distribution are negatively charged, while the larger grains become neutral or positively charged, as can be seen, for example, in Fig. 4.

Appendix B. Charged particle velocities
In a weakly-ionised plasma the momentum equations for a charged species i reduce to where ρ i is the density of the charged particles, ρ n the neutral density, α i the charge-to-mass ratio of the charged species, K in the collision coefficient between the neutrals and the charged species, v i the velocity of the charged species, v n the neutral velocity, E the electric field and B the magnetic field. Combined with Ohm's law, with r 0 the resistivity along the magnetic field, r 1 the Hall resistivity, r 2 the ambipolar resistivity and J the current, we can derive the charged velocities. It is often easier to work in the frame co-moving with the neutral fluid. Using u i = v i − v n we can rewrite Eq. B.1 as where we also substituted in the Hall parameter of the charged species β i = α i B/ρ n K in and E * = E + v n × B is the electric field in the neutral frame. This equation can also be expressed in matrix form as It is then straightforward to find the charged velocities from as long as the matrix is invertible. As the determinant of the matrix equals B 3 βi 1 + 1 β 2 i , we find that this is satisfied as long as β i is not too large. However, this cannot be guaranteed in simulations making the numerical scheme unstable.
In order to avoid this from occurring, a different choice of coordinate frame is necessary. We set the z-axis along the magnetic field and the x-axis along the current perpendicular to the magnetic field, i.e. B = (0, 0, B) and J = (J ⊥ , 0, J || ). This simplifies the expression for the electric field to E * = r 2 J ⊥ , −r 1 J ⊥ , r 0 J || and Equation B.4 to We can see that the velocity along the magnetic field can be calculated easily as while the velocities perpendicular to the magnetic field can be found from This solution shows that, for β i → ∞, u x,i = r 1 J ⊥ /B and u y,i = r 2 J ⊥ /B, while, for β i → 0, both u x,i and u y,i are zero and the charged particles move with the neutral, as expected. Furthermore, using J = i α i ρ i u i and substituting the charged velocities from Eqs. B.7-B.10, we find that J = (J ⊥ , 0, J || ), as required. Thus, this way of deriving the charged velocities is robust although it does requires extra calculations to rotate the frame of reference, first, along the local magnetic field and then back to the original frame.