Dust Dynamics in Kelvin-Helmholtz Instabilities

The Kelvin-Helmholtz instability (KHI) is a fluid instability which arises when two contacting flows have different tangential velocities. As shearing flows are very common in all sorts of (astro)physical fluid setups, the KHI is frequently encountered. In many astrophysical fluids the gas fluid in loaded with additional dust particles. Here we study the influence of these dust particles on the initiation of the KHI, as well as the effect the KHI has on the density distribution of dust species in a range of different particle sizes. This redistribution by the instability is of importance in the formation of dust structures in astrophysical fluids. To study the effect of dust on the linear and nonlinear phase of the KHI, we use the multi-fluid dust+gas module of the MPI-AMRVAC [1] code to perform 2D and 3D simulations of KHI in setups with physical quantities relevant to astrophysical fluids. A clear dependency on dust sizes is seen, with larger dust particles displaying significantly more clumping than smaller ones. Dust resides in vastly different locations in space, from tori around supermassive black holes to outflows of AGB stars and tails of comets. However, the dust particles have a large (and often unknown) variety of sizes and compositions. Furthermore, in most cases the density of the dust is several orders of magnitude lower than the local gas density, and as a consequence the contribution of dust on the dynamics has often been ignored. Here we gain insight in the consequence of having dust particles embedded in an astrophysical fluid by investigating the effect of dust in the KHI. 1 Problem setup and model To simulate the dynamics of the gas and dust in the KHI we use the MPI-AMRVAC [1] code, which has a multi-fluid gas+dust module. In the module the gas is treated as an Eulerian fluid. Multiple dust fluids can be added to the simulation, each representing different dust particle sizes and/or dust material densities. Each dust fluid follows the equations of a pressureless gas. The dust fluids interact with the gas by a drag force fd that couples the momentum equations of gas and dust fluids: fd = −(1 − α)πndρad∆u √ ∆u2 + vt , (1) with nd the dust particle density, ad the grain radius of fluid d, ∆u the difference between the gas and dust velocities (∆u = u − ud ) and vt ∝ √ p/ρ the thermal speed of the gas, with density ρ and pressure p. This drag force fd is a combination of the Epstein drag law for the subsonic regime and the Stokes law for the supersonic regime [2]. No dust-dust interaction is included, which is a reasonable assumption as in most typical astrophysical fluids the mean free path of dust-dust collisions is much larger than that of gas-dust and gas-gas interactions. The initial domain we set up for the 2D simulations has a rectangular shape and initially consists of three regions: an upper and a lower region with equal velocities, v0, in opposite directions along the x-axis, separated by a thin middle layer with (total) thickness D where the velocity of the flow varies ae-mail: Tom.Hendrix@wis.kuleuven.be EPJ Web of Conferences 46, 06003 (2013) DOI: 10.1051/epjconf/20134606003 © Owned by the authors, published by EDP Sciences, 2013 This is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Article available at http://www.epj-conferences.org or http://dx.doi.org/10.1051/epjconf/20134606003 EPJ Web of Conferences linearly from the upper to the lower velocity. The existence of the middle layer inhibits instabilities with wavelengths λmin . 5D for gas-only KHI [3]. The boundaries perpendicular to the x-axis are periodic. Perpendicular to the flow direction we use open boundaries. Depending on the simulation different initial perturbations have been used. Unless mentioned otherwise, we use a “smooth" perturbation of the perpendicular velocity near the middle layer: vy = U0 exp ( − (y − M) 2 2σ2 ) sin (kxx). (2) U0 sets amplitude of the perturbation and is set as U0 = 10v0. M is the y-coordinate of the middle of the separation layer. To make the perturbation smooth, we take σ = 5D, making the region influenced by the perturbation several times larger than the middle region. In a typical setup we use seven levels of adaptive mesh refinement (AMR) on a 16×32 base grid, resulting in an effective resolution of 1024×2048. When setting up the simulation to study the growth of an instability of a certain wavelength λ, we choose the physical dimensions of the box to be 20λ in the x-direction while keeping the size in the y-direction fixed. Doing so,we assure that in every simulation every wavelength of the perturbation is resolved with the same number of cells (viz. ∼ 52 cells per wavelength), independent of the actual wavelength. The aspect ratio of the axes is thus different for each simulation; in terms of code units, the lengths of the x and y axes vary between 0.17×0.3 up to 1.26×0.3. In these units, the middle layer has thickness D = 2 × 10−3. Initially the gas and dust fluids are uniformly distributed. We use a gas density ρ = 10−20 g cm−3, or ∼ 6× 103 mH cm−3. The temperature of the gas, used for the pressure calculation, is set to T = 100 K. In the upper and lower layer the gas has a velocity v0 = 5 × 104 cm s−1, in opposite directions. This is ≈ 0.43 times the speed of sound of the medium. We use an adiabatic index of γ = 5/3. The total dust density is chosen by setting the dust to gas ratio δ. We use ratios between 0.01 and 1, with δ = 0.01 the canonical value often used for the ISM. The velocities of the dust fluids are initially taken identical to those of the gas fluid. For each dust species one has to specify the material density of the dust grains, which is used to calculate the drag force. Here we take the dust grains to be composed of silicates, resulting in an internal grain density of ρp = 3.3 g cm−3 [4]. When setting up the dust fluids, a range in dust grain sizes is chosen. We typically choose to include grains of sizes between 5 nm and 250 nm. We then choose the number of bins to represent a certain part of this range in sizes. Each bin will be represented as a separate dust fluid during the simulation. In most simulations discussed here we use 4 dust fluids, unless mentioned explicitly. Once the number of dust fluids has been chosen, the size range for each bin is chosen by setting an equal density in each bin. This is done by integration over the dust density as computed from the size distribution, the grain volume and grain density: ∫ Rmax Rmin n(r) [ 4 3 πr3 ] ρp dr = 1

Dust resides in vastly different locations in space, from tori around supermassive black holes to outflows of AGB stars and tails of comets.However, the dust particles have a large (and often unknown) variety of sizes and compositions.Furthermore, in most cases the density of the dust is several orders of magnitude lower than the local gas density, and as a consequence the contribution of dust on the dynamics has often been ignored.Here we gain insight in the consequence of having dust particles embedded in an astrophysical fluid by investigating the effect of dust in the KHI.

Problem setup and model
To simulate the dynamics of the gas and dust in the KHI we use the MPI-AMRVAC [1] code, which has a multi-fluid gas+dust module.In the module the gas is treated as an Eulerian fluid.Multiple dust fluids can be added to the simulation, each representing different dust particle sizes and/or dust material densities.Each dust fluid follows the equations of a pressureless gas.The dust fluids interact with the gas by a drag force f d that couples the momentum equations of gas and dust fluids: with n d the dust particle density, a d the grain radius of fluid d, ∆u the difference between the gas and dust velocities (∆u = u − u d ) and v t ∝ p/ρ the thermal speed of the gas, with density ρ and pressure p.This drag force f d is a combination of the Epstein drag law for the subsonic regime and the Stokes law for the supersonic regime [2].No dust-dust interaction is included, which is a reasonable assumption as in most typical astrophysical fluids the mean free path of dust-dust collisions is much larger than that of gas-dust and gas-gas interactions.The initial domain we set up for the 2D simulations has a rectangular shape and initially consists of three regions: an upper and a lower region with equal velocities, v 0 , in opposite directions along the x-axis, separated by a thin middle layer with (total) thickness D where the velocity of the flow varies EPJ Web of Conferences linearly from the upper to the lower velocity.The existence of the middle layer inhibits instabilities with wavelengths λ min 5D for gas-only KHI [3].The boundaries perpendicular to the x-axis are periodic.Perpendicular to the flow direction we use open boundaries.Depending on the simulation different initial perturbations have been used.Unless mentioned otherwise, we use a "smooth" perturbation of the perpendicular velocity near the middle layer: U 0 sets amplitude of the perturbation and is set as U 0 = 10 −3 v 0 .M is the y-coordinate of the middle of the separation layer.To make the perturbation smooth, we take σ = 5D, making the region influenced by the perturbation several times larger than the middle region.In a typical setup we use seven levels of adaptive mesh refinement (AMR) on a 16×32 base grid, resulting in an effective resolution of 1024×2048.When setting up the simulation to study the growth of an instability of a certain wavelength λ, we choose the physical dimensions of the box to be 20λ in the x-direction while keeping the size in the y-direction fixed.Doing so,we assure that in every simulation every wavelength of the perturbation is resolved with the same number of cells (viz.∼ 52 cells per wavelength), independent of the actual wavelength.The aspect ratio of the axes is thus different for each simulation; in terms of code units, the lengths of the x and y axes vary between 0.17×0.3 up to 1.26×0.3.In these units, the middle layer has thickness D = 2 × 10 −3 .Initially the gas and dust fluids are uniformly distributed.We use a gas density ρ = 10 −20 g cm −3 , or ∼ 6 × 10 3 m H cm −3 .The temperature of the gas, used for the pressure calculation, is set to T = 100 K.
In the upper and lower layer the gas has a velocity v 0 = 5 × 10 4 cm s −1 , in opposite directions.This is ≈ 0.43 times the speed of sound of the medium.We use an adiabatic index of γ = 5/3.The total dust density is chosen by setting the dust to gas ratio δ.We use ratios between 0.01 and 1, with δ = 0.01 the canonical value often used for the ISM.The velocities of the dust fluids are initially taken identical to those of the gas fluid.For each dust species one has to specify the material density of the dust grains, which is used to calculate the drag force.Here we take the dust grains to be composed of silicates, resulting in an internal grain density of ρ p = 3.3 g cm −3 [4].When setting up the dust fluids, a range in dust grain sizes is chosen.We typically choose to include grains of sizes between 5 nm and 250 nm.We then choose the number of bins to represent a certain part of this range in sizes.Each bin will be represented as a separate dust fluid during the simulation.In most simulations discussed here we use 4 dust fluids, unless mentioned explicitly.Once the number of dust fluids has been chosen, the size range for each bin is chosen by setting an equal density in each bin.This is done by integration over the dust density as computed from the size distribution, the grain volume and grain density: with R min and R max the delimiters of a certain bin, N the number of bins and a min and a max the minimum and maximum size of the total range, respectively.The function n(r) represents the distribution function for dust particles of radius r.We use n(r) ∝ r −3.5 , which is a typical size distribution for dust particles in the ISM.Then we choose a single grain size r for each bin, which is representative for the bin, by taking the weighted mean over the bin in question, using the drag force (equation 1) as a weighting function.

Linear phase of the KHI
To study the growth of the dusty KHI, we compare the growth rates with those of gas-only KHI.The gas-only setup with two shear layers and a boundary layer with total thickness D is unstable between 0 < κ < 1.2785, with κ = k x D a dimensionless number and k x the wavenumber of the perturbation [3].The setup is most unstable around κ = 0.7968.We simulate three different setups, each setup has a different dust to gas ratio δ, namely δ = 0 (no dust, to compare the simulations with the theoretical result), δ = 0.01, δ = 0.1 and δ = 1.For each setup we investigate 11 wavelengths.
To quantify the growth of the instability through time, we use the total kinetic energy of the gas in the y-direction (viz.normal to the flow boundary in 2D and 3D) as an indicator.The velocity perturbation of the boundary layer initiates the KHI for the gas fluids, resulting in an exponential growth of the transverse kinetic energy of the gas.Even though the dust fluids are also perturbed, they are not prone to the KHI themselves.In figure 1 we compare the analytic solution of the dispersion relation with growth rates from simulations.We see that the δ = 0 simulations closely resemble the theoretical curve.The growth rates for wavelengths shorter than the most unstable wavelength (viz.κ > 0.7968) can be seen to be slightly higher than the theoretical growth.Adding dust with dust to gas ratio of δ = 0.01 only has minor consequences on the growth rates, which are ∼ 2% lower than for the gasonly simulation.Increasing δ to 0.1, the reduction in growth rate ranges between 5% and 20%, with a stronger decrease for larger κ.By interpolating the data we derive an approximation of the maximum value of the growth rates.As noted earlier the growth rate for the gas-only is higher than theoretical for shorter wavelengths, resulting in a maximum at κ = 0.81.For δ = 0.1 the maximum shifts slightly to smaller κ, having its maximum at κ = 0.79.Further raising the dust to gas ratio to δ = 1.0, we see a strong decrease in growth rates.As in the δ = 0.1 case, the decrease is strongest for the shorter wavelengths.Wavenumbers κ > 1.0 (or λ > 2πD) are stabilized.The most unstable wavelengths can be found around κ = 0.51 and growth rates are down by around 50%.During the linear phase the typical KH vortices are formed.The gas density distribution is similar to what is typically seen for the gas-only KHI, with alternating overdense and underdense islands where the pressure has maxima and minima, respectively.In the underdense vortices, the gas rotates forming closed elliptical streamlines.
Initially the rotational velocity of the gas in the vortices increases.After the linear phase the rotation period of the vortices stays constant, with a rotation period P ≈ 2λ/v 0 .Due to the rotational motion of the gas the dust is accelerated outwards.As the dust fluids are pressureless, only the interaction with the gas fluid counteracts the outward motion.Smaller dust particles (in our simulation the dust fluid with r ∼ 8 × 10 −7 cm) are more easily stopped as they have lower momentum.This difference in stopping time makes that the lightest dust is coupled to the gas fluid (the density distribution closely resembles the vorticity of the gas velocity), while heavier dust species are propagated further outward.The heaviest dust species are not stopped and form a layer of increased density around the vortices.
During the linear phase an accumulation of dust is observed in the overdense islands.

Non-linear phase and dust redistribution
In the non-linear phase, the separation of dust from the vortices continues.If we look at the global minima of the dust fluid densities (corresponding to centres of vortices) we notice that, for all the dust sizes, this quantity decreases exponentially after the end of the linear phase (see left figure 2).The decrease is stronger for heavier particles as they are more easily removed from the vortex due to their higher particle momentum.At the same time, we see that the global maximum of the dust density increases strongly after the linear phase (right in figure 2).First the accumulation is in the regions where the gas is overdense, later these regions break up and the maxima of the dust density is located in layers around the vortices.In later stages, the vortices become unstable and start to merge.Mergers can also increase local dust density by pushing together dust layers.
If we compare setups with different dust to gas ratios, we see that setups with lower initial values of δ will have a (slightly) stronger increase in maximum density, as compared to the initial density.If we look at the percentage of the volume in which the density has decreased by at least a factor of 2.5, we see that this quantity increases linearly during the non-linear phase for both δ = 0.01 and δ = 0.1.The linear increase is strongest for δ = 0.01, and slightly less strong for δ = 0.1.This trend of linear increase is observed until the end of our simulations at t = 120, long after the end of the linear phase.For δ = 1 the linear trend is less clear, and the increase rate is also slower.We have also performed simulations of the KHI in 3D, using 2 dust species in a square box setup with effective resolution of 256×1024×256, using random perpendicular velocities at the boundary layer.The size of the box is 2 times the most unstable wavelength in all directions.The boundaries are periodic in the x and z direction, and open in the y direction.In the linear stage the 3D case is quantitatively comparable with the 2D simulations.Later in the non-linear phase however, the vortex "tubes" start bending, and eventually become unstable in the extra third dimension, leading to a breakup of the vortices in smaller parts.This is in contrast to the 2D case, where the small vortices over time amalgamate into one large vortex.This additional instability compresses the dense dust layers around the vortex tubes, forming filament structures in which the dust density is enhanced up to 40 times, which makes the clumping of dust even stronger in 3D as compared to the similar 2D setup.

Figure 1 .
Figure1.Dependence of the growth rate of the KHI on the wavelength of the initial perturbation.The continuous line is the analytic solution of the dispersion relation of a gas-only fluid with the setup as described in 1.Each growth rate is calculated from a different simulation by fitting the linear growth phase of the transverse kinetic energy with an exponential.

Figure 2 .
Figure 2. Left: evolution of the global minimum in dust density with times, for all 4 dust fluids.The size of the assumed particles in the fluid increases from fluid 1 to 4. After the end of the linear phase (t ∼ 10) the global minimum decreases exponentially.The decrease is faster for heavier dust particles.Right: global maximum for the four dust species.Density increases are stronger for heavier dust species.