© ESO, 2014

1. Introduction

The Kelvin-Helmholtz instability (KHI) occurs at the interface of shearing fluids and taps energy from the velocity difference to create rotating structures in-between. It has been studied repeatedly, on the one hand, because it has important effects on the surrounding environment, such as mixing and redistribution of matter and energy on various scales; on the other hand, it is of relevance in many different branches of fluid studies, such as oceanic circulation (van Haren & Gostiaux 2010), winds on planet surfaces, magnetic reconnection in the solar corona (Lapenta & Knoll 2003), interaction between comet tails and the solar wind (Ershkovich 1980), astrophysical jets (Baty & Keppens 2006), and many others. In many fluids, however, the fluid does not consist of only one smooth mixture, but one can think of the flow as loaded with extra “impurities”, having different physical properties as the particles make up the flow itself. An example of such impurities are dust grains. Even though the number of dust particles per volume is typically low compared to that of the gas because of the large size and mass of the dust particles compared to typical fluid particles, they can still be of importance on the dynamics of the system. Furthermore, the dynamics of the dust itself can be of major interest, as for example in protoplanetary disks where the initially dilute dust particles embedded in the gas can separate itself from the gas and eventually clump to form planetesimals and later planets. However, how the dust can actually grow to large macroscopic sizes is not fully understood. To explain the dust growth in protoplanetary disks, high dust densities are required, which can be obtained through (the interplay of) several instabilities. Often an instability can dynamically raise the local dust density to values high enough to initiate a gravitational collapse. Meheut et al. (2012) used numerical simulations to show that the Rossby wave instability in disks can trap dust. Jacquet & Balbus (2012) discuss how dust in disks is influenced by the magnetorotational instability. In protoplanetary disks, the KHI can be of major importance. As the dust, which moves closer to the Keplerian velocity than the gas, settles gravitationally towards the midplane, it can accelerate the gas in the midplane to higher velocities, potentially causing a shear instability between the faster moving gas in the midplane and the slower gas above and below. In this way, the KHI can potentially disrupt the growths to higher density in the midplane, or dynamically cause higher densities in some regions. This has been studied analytically as well as numerically in Michikoshi & Inutsuka (2006), Johansen et al. (2006), and Barranco (2009).

In systems where dust is typically found, particles have a large (and often unknown) variety of sizes and compositions; therefore, including the effect of dust in the analysis of the dynamics of these systems causes a significant increase in numerical and analytical complexity. Furthermore, the density of the dust is several orders of magnitude lower than the local gas density in most cases, and as a consequence, the contribution of dust on the dynamics has often been ignored. In the last few years, several works have been presented with different numerical approaches towards the incorporation of dust. Examples range from fully particle-based two-fluid smoothed-particle hydrodynamics (SPH) methods (Cha & Nayakshin 2011; Laibe & Price 2012) to hybrid particle-fluid methods (Johansen et al. 2006; Miniati 2010) and two-fluid hydrodynamic models (Paardekooper & Mellema 2006; van Marle et al. 2011; Meheut et al. 2012).

Here, we gain insight in the consequence of having dust particles embedded in a fluid by investigating the effect of dust on the KHI. We do this by making 2D and 3D simulations of the dusty KHI with the fluid code MPI-AMRVAC (Keppens et al. 2012). While many previous studies that included dust dynamics treat the system as being composed of a two-fluid gas+dust, MPI-AMRVAC simulates an entire range of dust sizes simultaneously, allowing us to compare the effect of the KHI on dust of different size and vice versa. Additionally, while many previous gas+dust simulations are focused on a specific astrophysical setup, the results we obtain here are generally valid for many different systems with dust that form shear instabilities. How this is done is outlined in Sects. 2 and 3. In Sect. 4.1, we study how the addition of dust changes the linear phase of the KHI as compared to the analytic gas-only solution and the effect of different dust-to-gas density ratios on the growth rate and wavelength dependence of the instability. Further we discuss how different flow velocities change the growth rate of the dusty KHI as compared to the situation with only gas in Sect. 4.1. In Sect. 4.2, we discuss the dynamics after the linear phase and look at how the distribution of the gas and dust densities are changed due to vortex formation, mixing, and vortex merging. We outline how the KHI strongly alters the distribution of the dust and forms regions devoid of dust and layers where the dust density is increased significantly. We also discuss how the amount of dust in the fluid changes these effects. In Sect. 4.3, we discuss in which way 3D effects change the situation in the previous sections. Finally, we compare the results we found for the structure formation of dust in KHI with some findings on filaments observed in molecular clouds and highlight some interesting similarities in Sect. 5.

2. Problem setup

2.1. Numerical method

To simulate the KHI, we use the MPI-AMRVAC code (Keppens et al. 2012), a finite volume code suitable for solving any system of hyperbolic partial differential equations. The inclusion of a dust module in the code is discussed in van Marle et al. (2011). Both gas and dust are treated as separate fluids. To capture species-dependent dust dynamics, multiple dust fluids with different characteristics (e.g. different grain size or material densities) can be used simultaneously. The gas dynamics are governed by the following equations: $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }\mathrm{\left(}\mathit{\rho }{v}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\rho }{v}{v}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathit{p}\mathrm{=}\sum _{\mathit{d}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}{{f}}_{\mathit{d}}\mathit{,}\\ & & \frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\left[\mathrm{(}\mathit{p}\mathrm{+}\mathit{e}\mathrm{)}{v}\right]\mathrm{=}\sum _{\mathit{d}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}{v}\mathrm{·}{{f}}_{\mathit{d}}\mathit{,}\\ & & \mathit{e}\mathrm{=}\frac{\mathit{p}}{\mathit{\gamma }\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}\mathrm{+}\frac{\mathit{\rho }{\mathit{v}}^{\mathrm{2}}}{\mathrm{2}}\mathit{,}\end{array}$with ρ as the gas density, v the gas velocity, p the thermal pressure, and f_d the drag force per volume as applied on the gas by dust type d. The parameter N is the number of dust fluids, e is the total energy density of the gas. Equation (4) describes the two parts contributing to the energy density e - namely, the internal energy with γ the adiabatic index and the bulk kinetic energy of the gas. Equation (2) can be seen as the Euler equation for momentum conservation, with an added term on the right-hand side that represents the total force per volume acting on the gas. Similarly, Eq. (3) is the Euler equation for energy conservation supplemented with an extra term on the right-hand side which represents the work done by the dust on the gas due to the dragforce.

The dust fluids are treated as a pressureless gas: the internal energy of a dust grain only influences the surface temperature of the grain and has no influence on its movement (see LeVeque 2004, for a discussion on the dynamics of the pressureless gas equations). Furthermore, dust-dust collisions are typically rare as dust densities are often low in astrophysical fluids (the typical dust particle density in the ISM is ≈1 m^-3, or 10⁶ times lower than that of the gas), and the collisional cross section of grains is small. Even in regions where the dust density is enhanced or comparable to that of the gas density, the particles themselves are still rare compared to the amount of gas particles. Furthermore, dust-dust collisions at typical astrophysical velocities are often inelastic (see e.g. Hirashita & Yan 2009; and Ormel et al. 2009) and would not contribute to an internal pressure of a dust fluid. For the treatment of the dust fluids, the governing equations are thus $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{\rho }}_{\mathit{d}}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathit{d}}{{v}}_{\mathit{d}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \frac{\mathit{\partial }\mathrm{\left(}{\mathit{\rho }}_{\mathit{d}}{{v}}_{\mathit{d}}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathit{d}}{{v}}_{\mathit{d}}{{v}}_{\mathit{d}}\mathrm{\right)}\mathrm{=}\mathrm{-}{{f}}_{\mathit{d}}\mathit{.}\end{array}$Similar to the variables introduced for the gas fluid, ρ_d, v_d, and f_d are the dust density, bulk dust velocity and drag force of dust fluid d, respectively. We have compared our implementation of the gas+dust hydrodynamics with the benchmarks proposed in Laibe & Price (2012) and find our solutions in agreement with theirs.

For the drag force f_d, we use a combination of the Epstein drag law for the subsonic regime and the Stokes law for the supersonic regime (Kwok 1975), namely ${{f}}_{\mathit{d}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\mathit{\pi }{\mathit{n}}_{\mathit{d}}\mathit{\rho }{\mathit{a}}_{\mathit{d}}^{\mathrm{2}}\mathrm{\Delta }{v}\sqrt{\mathrm{\Delta }{{v}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{t}}^{\mathrm{2}}}\mathit{,}$(7)with n_d as the dust particle density, a_d the grain radius of fluid d, Δv the difference between the gas and dust velocities (Δv = v − v_d), and ${\mathit{v}}_{\mathrm{t}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\sqrt{\mathrm{3}\mathit{p}\mathit{/}\mathit{\rho }}$ (Kwok 1975) the thermal speed of the gas. The force is evaluated separately in each spatial direction and the same expression is used independent of the dimensionality of the problem. In this way, the 2D simulations represent a 3D simulation with embedded 3D grains in the limit of vanishing third dimensional velocities and not a truly physical 2D setup with disk-like dust particles, which would represent infinitely long cylinders if expanded to 3D. This way we can directly compare 2D simulations with cut of 3D simulations in the same plane. The variable α is the sticking coefficient and expresses the amount of gas particles that stick to the dust grain after collision. Thus, (1 − α) is a measure of the effectiveness of momentum transfer from gas to dust and vice versa. Following Kwok (1975), we use α = 0.25.

Sometimes, it is useful to represent regions as being free of dust, or in some cases, the dust density can become very low due the dynamics. To be able to handle these cases numerically, we introduce a threshold value for minimum density. If the density in a fluid would become lower than the threshold value, the dust from this fluid is removed by setting both the density and the velocity of the fluid to zero. In all simulations described here, the threshold is 10 orders of magnitude below the initial dust density.

Several numerical methods for the time advance are implemented in the code, as well as multiple flux limiters. In the simulations described in this paper we used the second order TVDLF (total variation diminishing Lax-Friedrich) scheme in combination with a monotonized central (MC) type limiter (van Leer 1977).

2.2. Numerical setup

The initial domain we set up for the 2D simulations has a rectangular shape and initially consists of three regions: an upper and lower region with equal velocities, v₀, in opposite directions along the x-axis, and a thin middle layer that separates them with (total) thickness D, where the velocity of the flow varies linearly from the upper to the lower velocity. The existence of the middle layer inhibits instabilities with wavelengths λ_min ≲ 5D for incompressible gas-only KHI (Chandrasekhar 1961). This guarantees that the growth of the KHI is not limited by the smallest resolution scale, which is contrary to a setup of KHI without a middle layer, gravity, and surface tension in which the initial growth increases as ω ∝ k_x(=2π/λ) for all wavenumbers (Chandrasekhar 1961). The boundaries perpendicular to the x-axis are periodic, limiting the development of instabilities of wavelengths longer than the box size; however, as we will typically take a box several times larger than the most unstable wavelength, this will not prove a limitation. Perpendicular to the flow direction, we use open boundaries. Care is taken to assure that these boundaries are always far enough from the instability dominated region during the entire simulation.

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: ${\mathit{v}}_{\mathit{y}}\mathrm{=}{\mathit{U}}_{\mathrm{0}}\exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{y}\mathrm{-}\mathit{M}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\sin \mathrm{\left(}{\mathit{k}}_{\mathit{x}}\mathit{x}\mathrm{\right)}$(8)for both the gas and dust fluids. The parameter U₀ sets the amplitude of the perturbation and is set as U₀ = 10^-3v₀ and M is the y-coordinate of the middle of the separation layer. To make the perturbation smooth, we take σ = 5D, which makes the region influenced by the perturbation several times larger than the middle region. Another approach we used to set up the perturbations is to fill the middle layer with different random v_y values for each cell: ${\mathit{v}}_{\mathit{y}}\mathrm{=}{\mathit{U}}_{\mathrm{0}}\left[\mathrm{2}\hspace{0.17em}\mathrm{rand}\mathrm{\left(}\mathrm{\right)}\mathrm{-}\mathrm{1}\right]$(9)with rand() being a random number between 0 and 1 and U₀ having the same value as before. All fluids have the same initial velocities and the same perturbation velocities.

In a typical setup, we use seven levels of adaptive mesh refinement (AMR) on a 16 × 32 base grid, which results in an effective resolution of 1024 × 2048. The AMR refinements are triggered by a Löhner type estimator (Löhner 1987), which computes normalized second derivatives to locate strong variations in specified variables. In all simulations in this paper, the AMR is triggered by a combination of the gas velocity perpendicular to the initial flow (caused initially by the KHI), the gas density ρ, and the density of the dust species representing the largest particles (see Sect. 3.2). Figure 1 gives an example of the refinement in a typical 2D simulation. By tracking the time evolution of the actual number of cells as compared to the number of cells needed in a similar simulation without AMR (see Fig. 2), we can compute that the simulations in Fig. 2 would take roughly 7.12 and 5.83 times the current time for the 2D and 3D run, respectively. A typical 2D simulation up to t = 120 uses ~200 CPU h, and the 3D simulation in Fig. 2 needs about 15 000 CPU h to complete. More details about the 3D simulation can be found in Sect. 4.3.

Fig. 1 Full domain of a 2D simulation at t = 40. The background colour represents the refinement level with lighter colours being higher levels of refinement. While seven levels are used in the simulation, the coarsed level is no longer used at this time. In the middle, the gas density is overlaid to illustrate how the AMR chooses maximal refinement around the vortices.

Fig. 2 Evolution of the ratio between the actual amount of cells in a simulation and the amount needed to fill a uniform grid with the same resolution. The difference between the 2D and 3D setup are due to differences in the extent of the simulated domain.

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), which is 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.

3. Fluid properties

The simulations of the dusty KHI are performed with an astrophysical fluid in mind, and the fluid properties have been chosen accordingly. However, the physical quantities in the code itself are scaled to values close to unity, as the Euler equations used to describe the fluid can be written in dimensionless form. Therefore, these simulations are thus also equally relevant for other environments.

3.1. Gas fluid

Initially, the gas fluid is uniformly distributed over the entire grid. We use a gas density ρ = 10^-20 g cm^-3, or ~6 × 10³ m_H cm^-3. The temperature of the gas, which is used for the initial pressure calculation (assuming an ideal gas EOS), is initially set to T = 100 K. In the upper and lower layer, the gas has a velocity v₀ = 5 × 10⁴ cm s^-1 in opposite directions (different values are used in Sect. 4.1.2). This is ~0.43 times the speed of sound of the medium. We use an adiabatic index of 5/3, which is appropriate for the used values of gas density and temperature (Tomida et al. 2013). A summary of the physical values is given in Table 1.

Fig. 3 Evolution of the transverse kinetic energy of the gas in six different simulations with the same dimensionless wavenumber κ = 0.7968 (the most unstable wavenumber in the gas-only case). The vertical axis has been set to logarithmic scale to clearly show the exponential growth in the initial linear growth phase. Three simulations start with a lower initial transverse kinetic energy. This is because they start with small random perturbations in the middle layer, instead of being excited with the perturbation described in Eq. (8). These three simulations are discussed in the Sect. 4.3. The behavior of the simulation with δ = 1 is different; this is explained in Sect. 4.1.1.

3.2. Dust fluid

The total dust density (the sum of the densities of the dust fluids) is chosen by setting the dust-to-gas (mass density) ratio δ. Typically, we use ratios between 0.01 and 1. A canonical value for the dust-to-gas ratio in the ISM of 0.01 (Spitzer 1954) is often used, so the ratios we consider are standard mixtures up to blends which are enriched with dust. The velocities of the dust fluids are initially taken identical to those of the gas fluid, including in the perturbed regions. For each dust species, one has to specify the material density of the dust grains ρ_p (the density of the matter inside a spherical dust grain), 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 (Draine & Lee 1984). 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 is represented as a separate dust fluid during the simulation. In most simulations discussed here, we use four dust fluids; however, this is mentioned explicitly for simulations with another number of subdivisions. 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: $\begin{array}{ccc}\frac{\mathrm{1}}{\mathit{N}}{\mathrm{\int }}_{{\mathit{a}}_{\min }}^{{\mathit{a}}_{\max }}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\left[\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\right]{\mathit{\rho }}_{\mathit{p}}\hspace{0.17em}\mathrm{d}\mathit{r}\mathrm{=}{\mathrm{\int }}_{{\mathit{R}}_{\min }}^{{\mathit{R}}_{\max }}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\left[\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}\right]{\mathit{\rho }}_{\mathit{p}}\hspace{0.17em}\mathrm{d}\mathit{r}& & \end{array}$(10)with R_min and R_max as 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 (Draine & Lee 1984), which is a typical size distribution for dust particles in the ISM in Eq. (10) to calculate the upper limit R_max of a bin with lower limit R_min: $\begin{array}{ccc}\mathrm{\Rightarrow }\hspace{0.17em}{\mathit{R}}_{\max }\mathrm{=}{\left[\sqrt{{\mathit{R}}_{\min }}\mathrm{+}\frac{\mathrm{1}}{\mathit{N}}\mathrm{(}\sqrt{{\mathit{a}}_{\max }}\mathrm{-}{\sqrt{{\mathit{a}}_{\min }}}^{\mathrm{)}}\right]}^{\mathrm{2}}\mathit{.}& & \end{array}$(11)The R_max of the first bin can be calculated immediately, further bin widths can be calculated by using Eq. (11) with the width of the first bin. Once the range has been further divided in N bins with for with each bin a (different) width, we need to choose a single grain size for each bin, which is representative for the bin during the simulation. This is done by taking a weighted mean over the bin in question. As a weighting function, one can take the particle size or particle mass; however, we choose a function, which corresponds to the drag force (Eq. (7), with a_d = r): $\begin{array}{ccc}{{f}}_{\mathit{d,}{\Delta }{r}}\mathrm{=}\underset{{\mathit{R}}_{\min }}{\overset{{\mathit{R}}_{\max }}{\mathrm{\int }}}{{f}}_{\mathit{d}}\hspace{0.17em}\mathrm{d}\mathit{r}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\mathit{\pi \rho }\mathit{r̅}\mathrm{2}\mathrm{\Delta }{v}\sqrt{\mathrm{\Delta }{{v}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{t}}^{\mathrm{2}}}\underset{{\mathit{R}}_{\min }}{\overset{{\mathit{R}}_{\max }}{\mathrm{\int }}}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}& & \end{array}$(12)viz. the total force on the bin with range from R_min to R_max, called f_d,Δr, is the integral over all the grain sizes in the bin. This total force is supposed to be equal to an equivalent force in which the particle size is replaced by the representative size . Thus, isolating in Eq. (12), canceling all parameters outside the integrals, and assuming the same particle size distribution as before gives us $\begin{array}{ccc}& & \mathit{r̅}\mathrm{2}\mathrm{=}\frac{\underset{{\mathit{R}}_{\min }}{\overset{{\mathit{R}}_{\max }}{\mathrm{\int }}}{\mathit{r}}^{\mathrm{2}}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}{\underset{{\mathit{R}}_{\min }}{\overset{{\mathit{R}}_{\max }}{\mathrm{\int }}}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}\\ & & \mathrm{\Rightarrow }\hspace{0.17em}\mathit{r̅}\mathrm{=}\sqrt{\mathrm{5}\frac{{\mathit{R}}_{\max }^{-0.5}\mathrm{-}{\mathit{R}}_{\min }^{-0.5}}{{\mathit{R}}_{\max }^{-2.5}\mathrm{-}{\mathit{R}}_{\min }^{-2.5}}}\mathrm{·}\end{array}$This method assures us that the drag force on all grain sizes is correctly represented by a force on the representative grain size in the bin. An overview of the resulting representative grain sizes for different amounts of dust fluids is given in Table A.1 in the appendix.

All dust fluids interact with the gas fluid through the drag force (Eq. (7)); however, dust-dust interactions are not included. One can demonstrate that the mean free path of a dust grain, l, can be written as $\mathit{l}\mathrm{=}\frac{{\mathit{\rho }}_{\mathit{p}}{\mathit{a}}_{\mathit{d}}}{\mathrm{3}\sqrt{\mathrm{2}}{\mathit{\rho }}_{\mathit{d}}}$(15)with ρ_p the internal grain density. The mean free path increases with grain size for a certain dust density ρ_d. For the values of δ and a_d mentioned above, the mean free path ranges between 5.9 × 10¹³ cm and 2.9 × 10¹⁷ cm, while dust-gas collisions only have a mean free path between 2.4 × 10⁵ cm and 6.0 × 10⁸ cm. Thus, the exclusion of dust-dust interactions is a valid assumption. As mentioned earlier, the inclusion of dust-dust interactions would involve taking into account physical phenomena, such as dust coagulation and shattering, which are not of the essence here but will be included in the code in the future.

In Appendix A, a demonstration of the validity of the approach described in the previous section is given by comparing the structures and properties of setups with different amounts of dust fluids.

4. Results

4.1. Linear phase

To quantify the growth of the instability through time in the simulations, we use the total kinetic energy of the gas in the y-direction (viz. normal to the initial background flow in 2D and 3D) as an indicator. From this measurement, we see how the instability in the boundary layer grows as free energy from the unstable setup is converted in movement perpendicular to the background flow. Several examples of the evolution of the transverse kinetic energy can be seen in Fig. 3, where six evolutions are compared. In the figure, we see how the kinetic energy starts to increase exponentially after a short settling period. This phase, called the linear phase, ends around t ~ 10−25, depending on the simulation. Each growth rate is calculated from a different simulation by fitting the linear growth phase of the transverse kinetic energy with an exponential. We then adjust the values of the growth rate found by the analytic analysis (see Sect. 4.1.1) so that it can directly be compared with the growth of the kinetic energy by converting the analytic growth rate of velocity to one of kinetic energy (i.e., there is an extra factor of two as the analytic growth rate is valid for v_y and the kinetic energy is dependent on ${\mathit{v}}_{\mathit{y}}^{\mathrm{2}}$).

During the linear phase, the spatial evolution of the velocity perpendicular to the flow, v_y, can be found as the corresponding eigenfunction. The two-dimensional form of v_y is described by the Taylor-Goldstein equation (Taylor 1931; Goldstein 1931), $\mathrm{(}{\mathit{v}}_{\mathrm{0}}\mathrm{-}\mathit{c}\mathrm{)}\left[\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\mathrm{-}{\mathit{k}}^{\mathrm{2}}\mathit{\psi }\right]\mathrm{+}\left[\frac{{\mathit{N}}^{\mathrm{2}}}{{\mathit{v}}_{\mathrm{0}}\mathrm{-}\mathit{c}}\mathrm{-}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{v}}_{\mathrm{0}}}{\mathit{\partial }{\mathit{y}}^{\mathrm{2}}}\right]\mathrm{=}\mathrm{0}\mathit{,}$(16)with v₀ as the initial flow as described above; ψ a two-dimensional stream function with $\mathit{\psi }\mathrm{=}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial y}}$ and $\mathit{\psi }\mathrm{=}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{y}}}{\mathit{\partial x}}$; c the phase speed; and $\mathit{N}\mathrm{=}\sqrt{\mathrm{-}\frac{\mathit{g}}{\mathit{\rho }}\frac{\mathit{\partial \rho }}{\mathit{\partial y}}}$ the Brunt-Väisälä frequency (g is a (gravitational) acceleration in the y-direction). It can be seen that Eq. (16) can easily be solved for the three parts of the domain and coupled at the interface, as both N and the second derivative of v₀ vanish everywhere in our setup. Figure 4 shows the comparison between the values of v_y in a simulation with δ = 0.01 at t = 3 and the eigenfunctions in the three parts, as calculated from the gas-only version of the Taylor-Goldstein equation.

Fig. 4 Comparison of the values of vy of a simulation (+ symbols) and the derived eigenfunctions (dotted lines) on a cut perpendicular to the flow direction. For clarity, only a limited amount of simulated points is shown, The values of vy are given in cm s-1. The cut, taken at t = 3, only spans a small portion of the range along the y-axis, and is centred on the middle of the boundary layer. The distance from the centre is given in units of D. The dotted lines each represent a solution of the eigenfunction, the red line in the upper flow, teal in the middle layer, and blue in the lower flow. From the correspondence with the simulated values, it is clear that during the perpendicular velocity the linear phase evolves from the perturbation in Eq. (8) to the distribution expected from the analytical derivation of the eigenfunction.

4.1.1. Wavelength dependence of the growth rate

To study the growth of the dusty KHI, we compare the growth rates with those of gas-only KHI. As derived in Chandrasekhar (1961) for an incompressible gas, the gas-only setup with two shear layers at ± D/2 and a boundary layer with total thickness D is unstable between 0 < κ < 1.2785 with κ = k_xD, a dimensionless number, and k_x the wavenumber of the perturbation. By solving the analytically known dispersion relation, one can find that the setup is most unstable around κ = 0.7968. As described in Sect. 2.2, each wavelength is investigated by setting up a simulation with length 20λ in the x-direction. For each wavelength, we simulate four 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 different wavelengths. As mentioned in Sect. 4.1, we make quantifications of the growth rate of the transverse kinetic energy from the simulations.

The initial transverse velocity perturbation of the middle layer initiates the KHI for the gas fluids, resulting in an exponential growth of the initial perturbation for the gas. Even though the dust fluids are also initially perturbed, they are not prone to the KHI themselves, as the fluids are pressureless. In Fig. 5, we compare the analytic solution of the dispersion relation with growth rates derived from simulations. It can be seen that the gas-only simulations closely resemble the theoretical curve. The growth rates for wavelengths shorter than the most unstable wavelength (viz. κ larger than 0.7968) can be seen to be slightly higher than the theoretical growth for an incompressible gas, which we attribute to finite compressibility effects; the deviation from uniform gas density, Δρ/ρ, is typically up to 5−10% near the vortices at the end of the linear phase. Adding dust with a total dust-to-gas ratio of δ = 0.01 only has minor consequences for the growth of the KHI: growth rates are ~2% lower than for the gas-only simulation. We also note that the simulations with dust follow the trend of having raised growth rates compared to the theoretical gas-only curve for higher values of κ. Increasing δ to 0.1, we see that the effect of the added dust increases. The reduction in growth rate ranges between 5% and 20% with a stronger decrease for shorter wavelengths. By making a fourth order spline interpolation of the data, we can derive an approximation of the maximum value of the growth rate from the simulations. As noted earlier the growth rate for the gas-only is higher than theoretical for shorter wavelengths, resulting in a derived maximum at κ = 0.81. For the δ = 0.1 simulations, this maximum shifts slightly to longer wavelengths, 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 seen before for the δ = 0.1 case, the decrease is strongest for the shorter wavelengths. Wavenumbers larger than κ = 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%. An example of the slower growth of the δ = 1 simulations can be seen in Fig. 3. The simulation shown there is excited with a κ = 0.7968 perturbation, and between t = 6 and t = 16, the perturbation grows exponentially. The figure also shows clearly that the growth is significantly slower than in the other simulations. Between t = 18 and t = 28, a second settling occurs in which a mode with longer wavelength starts to grow. After t = 28, a new linear phase starts; this time with a longer wavelength, where we initially set up the simulation to include 20 perturbations with wavenumber κ = 0.7968. We now see only 12 growing perturbations. This means the wavenumber effectively shifts from κ = 0.7968 to κ = 0.48 at which the growth rate is faster, as can be seen in Fig. 5. Note that this is not a physical merger of vortices but the growth of the larger wavelength as it overtakes the growth of the shorter one.

The reason for the shift towards smaller wavenumbers can be interpreted by noting that the total inertia of the fluid caused by both dust and gas particles, increases by increasing δ. In contrast, the amount of gas, which causes the KHI, remains the same. Not only is the development of the perturbation velocity in the y-direction slowed by the added dust that needs to be dragged along, the added inertia of the dust moving in the background flow means that a longer time is needed to deviate a volume element from a movement in the x-direction to a rotational movement in a vortex, and therefore smaller vortices are less effective. This is similar to the effect of a faster background flow, which we discuss in the next section.

Fig. 5 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 Sect. 2.2.

4.1.2. Velocity dependence of the growth rate

In Sect. 4.1.1 we discussed how the growth rate of the KHI is influenced by the wavelength of the imposed perturbation. Another important quantity that affects the growth rate of the instability is the velocity difference of the two fluids. A comparison between the setup without dust and with dust (δ = 0.1) is given in Fig. 6. In Sect. 4.1.1, we saw that the growth rate for the δ = 0.1 is always lower than for the gas-only case. When we keep the perturbation wavelength fixed (to κ = 0.7968) but change the initial flow velocity v₀ we see that the growth rate of the KHI for all values of v₀ is once again decreased as compared to the gas-only simulations. For low velocities, growth rates are almost the same with and without dust, but differences increase for larger velocities. The growth rates reach a maximum around Δv₀ = 2v₀ = v_s (with v_s the speed of sound). After the maximum, the growth rate decreases, and we observe that the KHI is stabilized when the flow velocity v₀ reaches the transonic regime. For supersonic values of Δv₀, the nonlinear effect ultimately allows the formation of shocks (see Fig. 7), which can be seen to emanate from the vortices (see e.g. Baty et al. 2003). We find that these shocks are only seen for fluids with small grains for the dust fluids. This is because small grains are coupled tightly to the dynamics of the gas fluid, while the higher inertia of large grains allows them to cross shocks more easily (van Marle et al. 2011).

Fig. 6 Comparison of the growth rates of simulations with different fluid velocities. The Mach number indicated is the Mach number of the velocity difference of the upper and lower flow. Fourth order spline fits have been added to guide the eye. Dust tends to stabilize the KHI for all velocities with an increasing effect at higher velocities. The growth rate increases until the velocity difference Δv becomes supersonic.

Fig. 7 Vortices that have formed shock structures at t = 35 in a simulation with Δv0 = 1.51vs, κ = 0.7968 and δ = 0.1. Shocks are clearly visible for the gas (left panel) and the smallest dust species (middle panel). In the second smallest dust species (right panel), the shocks can still be discerned; however, they are less evident. In larger dust species, these shocks are not seen.

Fig. 8 Density distribution of the largest dust grains out of four species at t = 10 in a 2D simulation with δ = 0.01 and κ = 0.7968. Only part of the simulated domain is shown. In the centre of the middle vortex, the black ellipse is a streamline of the gas fluid. The gas trajectory is (almost) closed as the gas fluid is pressure supported by the lower pressure in the centre of the vortex. By contrast, the streamlines of the large dust particles (shown in red) can be seen to spiral out of the vortex. Ultimately, the dust ends up in two layers below and above the vortex structures.

Fig. 9 Left: evolution of the global minimum in dust density with time for all 4 dust fluids in a simulation with δ = 0.01 and κ = 0.7968. The size of the assumed particles in the fluid increases from dust 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. Only the part of the simulation is shown, where all minima are above the allowed minimum density (see Sect. 2.1). Right: global maximum for the four dust species. Density increases are stronger for heavier dust species.

4.1.3. Dust separation

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, corresponding to regions where the pressure has maxima and minima, respectively. In the underdense vortices, the gas rotates forming closed elliptical streamlines (see Fig. 8). For each vortex, we can study the rotational speed, ${\mathit{\omega }}_{\mathrm{rot}}\mathrm{=}{\mathit{v}}_{\mathrm{rot}}\mathit{/}\mathit{r}$(17)by looking at the velocities along several of these streamlines (v_rot) at multiple distances r from the centre of the vortex. Initially, the rotational speed of the gas in the vortices increases, and the fluids in the vortex have a differential motion: at t = 5, the rotational speed drops about 50% from around the centre to the edge of the vortex (where the density is the same as the ambient density). Near the end of the linear phase, the rotational speed reaches a uniform value over the entire vortex. After the linear phase the rotation period of the vortices stays constant with a rotation time T ≈ 2λ/v₀. 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. The smallest dust particles (in our simulation the dust fluid with cm, see Table A.1) are more easily stopped (The stopping time of grains due to dust-gas collisions scales linearly with the grain size.) as the particles have lower inertia. Indeed, the stopping time in a simulation with four dust fluids is ~24 times larger as that of the smallest ones. This difference in stopping time ensures 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. During the linear phase, an accumulation of dust is observed in the overdense islands.

4.2. Non-linear phase

4.2.1. Local dust density increase

In the non-linear phase, which begins around t ~ 10 in Fig. 3, the separation of dust from the vortices continues. The heaviest dust species are not stopped during their outward motion and form a layer of increased density around the vortices (see Fig. 8) by strongly lowering the dust density inside the vortices. If we look at the global minima of the dust fluids (corresponding to centres of vortices) we notice that this quantity decreases exponentially for all the dust sizes after the end of the linear phase (see left Fig. 9). The decrease is stronger for heavier particles as they are more easily removed from the vortex due to their higher particle inertia. At the same time, we see that the global maximum of the dust density increases strongly after the linear phase (right in Fig. 9). In the beginning, the accumulation is in the regions where the gas is overdense. Later, these regions break up and the maxima of the dust density are located in layers around the vortices. In later stages, vortices become unstable and start to merge. Mergers can also increase local dust density by pushing together dust layers. The dynamical breaking and merging of overdense layers is responsible for the minima and maxima in the right image of Fig. 9. Figure 10 gives a comparison of the gas density with the dust densities of the heaviest and lightest species. The vorticity of the gas is also shown to demonstrate the strong connection between the gas dynamics and the density distribution of the lightest dust species.

Fig. 10 Density distribution of the gas (top left), the lightest dust species (top right), and the heaviest dust species (bottom left), as well as the vorticity of the gas fluid (bottom right) at t = 60 in a simulation with δ = 0.01 and κ = 0.7968. Only part of the total domain is shown. In the image, the merging of two vortices is ongoing. A clear connection can be seen between the vorticity of the gas and the distribution of the smallest dust species.

If we compare setups with different dust-to-gas ratios, we see that setups with lower initial values of δ have a (slightly) stronger increase in maximum dust density, as compared to the initial density. The maxima and minima of the fluids can, of course, only describe the local behavior in the fluids. For the global density dynamics of the dust fluids we see that this quantity increases linearly during the non-linear phase for both δ = 0.01 and δ = 0.1 if we look at the percentage of the volume in which the density has increased by at least a factor of 2.5 (see Fig. 11). This trend of linear increase is observed until the end of our simulations at t = 120 long after the end of the linear phase. The increase is strongest for δ = 0.01 and slightly less strong for δ = 0.1. For δ = 1, the increase rate is slower. The reason for the observed increase is the merging of vortices. As the vortices merge, the surrounding dust is pushed together, leading to further dust enhancement. This also explains the decreasing growth of volume percentage with dust-to-gas ratio because the dust slows the merging process: the added mass of the dust decelerates the merging motion, which is initiated by the growth of the subharmonic instability of the gas, because it needs to be pulled along by the gas through the drag force, which transfers kinetic energy of the acceleration from the gas to the dust. At the end of the simulation shown in Fig. 11, the initial 20 vortices have merged to 2, 3, and 7 vortices for δ = 0.01, 0.1, and 1.0, respectively.

4.3. 3D effects

To expand and compare the previous results, we have performed simulations of the dusty KHI in 3D. The setup is comparable to ones in 2D with physical quantities as described in Sect. 3. The simulation described here has a dust-to-gas ratio δ = 0.1. While this value is higher than that typially assumed for the ISM, most conclusions also hold for δ = 0.01, as the dynamics in 2D is seen to be largely comparable. Instead of starting with a setup with a fixed wavelength perturbation and a domain with length 20λ, we adopt a simulated domain of (2λ_c)³ with λ_c as the most unstable wavelength (κ = 0.7968, see Sect. 4.1.1). In physical units, this translates to a domain of ~(0.01 pc)³. The smooth perturbation in Eq. (8) is replaced with a random perturbation of v_y in the middle plane. The resolution is set to 256 × 1024 × 256, which is four times higher than the resolution perpendicular to the contact plane. One can see that the spatial resolution is actually higher than in the 2D case. To reduce the computational cost, the number of dust species has been reduced to two. The boundaries in the x- and z-directions are periodic. The boundaries in the y-direction (parallel to the contact plane) have open outflow conditions. AMR is used; in the middle region however the grid is forced to the highest refinement during the entire simulation to allow small scale perturbations in the y- and z-directions to develop from the initial random perturbations. In the surrounding regions, the AMR follows variations in the density of the two dust fluids. While variations in the gas fluid is not traced actively, (because the dust fluid with the smallest dust particles is strongly coupled to the gas fluid) a reasonable grid coverage of the domain of the gas fluid is expected.

We choose to end the simulation at t = 40, as this value allows us to explore a significant timescale of the non-linear phase. On the other hand, turbulent motions time start filling the simulated domain at this time and increasing refinement (see Fig. 2) causes the simulation time to increase strongly. No stationary state has been reached at this point.

Fig. 11 Volume percentage of the simulated domain where the dust density has at least increased to 2.5 times the initial dust density. All simulations have v0 = 0.42 and are perturbed with wavelengths close to their most unstable value (see Fig. 5). The portion of the domain where the dust density is significantly enhanced is seen to increase linearly with time after the linear phase (t ~ 10 for δ = 0.01 and δ = 0.1 and t ~ 20 for δ ~ 1.0). The rate of increase decreases with dust-to-gas ratio δ.

Fig. 12 Volume rendering after the end of the linear phase (t = 22.5). The entire domain of the simulation is shown. The volume rendering shows the gas distribution of the two straight vortex-tubes that have formed. To visualize the vortices in the volume rendering, we only show gas densities below 9.2 × 10-21 g cm-3. At the back and left side of the cube (xy- and yz-planes, respectively), the dust distribution of the heaviest dust species is shown. The dust forms tubes around the gas vortices, and sheets of increased density between the tubes.

Fig. 13 Volume rendering at t = 40. To visualize the vortices in the volume rendering, we only show gas densities below 8.6 × 10-21 g cm-3. The straight vortex tubes seen at the end of the linear phase in Fig. 12 have become severely twisted and are breaking apart. The dust distribution at the same time is shown in Fig. 17.

Fig. 14 Comparison of the dust enhancement in four different simulations with different values of δ. The graph shows the maximum of the dust density in the simulation, ascompared to the initial value of the dust density. Dust enhancements are (initially) slightly higher for lower values of δ. The enhancements in the 3D case are significantly stronger due to the additional instability of the vortex tubes in the added third dimension. A log-scale is used for the maximum density. Note that the maximum density starts to increase at a later time in the 3D simulation as compared to the 2D cases because a random initial perturbation is used in 3D.

Fig. 15 Power per lengthscale derived from the 3D Fourier transform of the densities of the three fluids at t = 40. Formed structures are smaller for the dust species, with dust 1 being the dust species representing the smaller dust particles, and dust 2 representing the larger particles. The ripples at high wavenumber for the gas are due to the AMR, as we choose the active refinement to follow the density variations of the two dust fluids.

Fig. 16 Comparison of the strength of each mode of vy along the z-direction of the middle plane. Modes 0 to 13 are shown. At each time, the strength of each mode is compared to the strength of the dominant mode at that time. Three distinct regimes can be discerned in the figure. Initially, the initial random perturbations evolve up to t ~ 8 from small-scale structures to a large scale feature. In the linear phase and the early non-linear phase (between t ~ 8 and t ~ 30), the vortex tubes are elongated along the z-direction and the 0-mode is dominant. After t ~ 30, the vortex tubes become unstable, resulting in the formation of small-scale structures.

To be able to analyze the results in 3D, we have made a 2D simulation with the same setup as described above and for the same simulation in 3D without dust. In the linear part of the instability, the 3D KHI behaves similar to what we have seen in the 2D cases in the previous sections. The observed growth rate in 3D is almost exactly the same as in the 2D case, as can be seen in Fig. 3 (i.e. it is 1% slower). Also in the same line as found earlier in 2D, the growth in 3D without dust is 8% faster than with the added dust. As is expected from the setup, initially two vortices are formed. In the linear phase and the early non-linear phase, these vortices are straight, tube-like structures in the added z-direction, which adds no 3D-effects to the simulation (see Fig. 12). Similar to the heavy dust layers in 2D, sheets of increased dust density are formed around the vortices and between the vortices in the braids. At around t = 30, a secondary instability occurs (see Fig. 16). A discussion of secondary KH instabilities which are inherently 3D effects in the KHI can be found in Klaassen & Peltier (1991); Caulfield & Peltier (2000) for stratified and homogeneous flows. Here, the three-dimensional instability suppresses the merging mode of the vortices and results in a twisting motion of the vortices, lifting the tubes out of the middle plane and causing a turbulent-like state as discussed in Metcalfe et al. (1987). We also find the secondary instability in the gas-only 3D simulation, as is expected because it is not triggered by the dust itself. The twisting motion causes the dust sheets to be folded and pushed together, leading to the formation of elongated, filamentary dust structures with enhanced dust densities. If we compare the dust density enhancement of simulations with different values of δ with the enhancement in the 3D simulation (Fig. 14), a clear difference can be seen. Where in the 2D cases, enhancements up to a factor of 10 are typical; the folding of the dust sheets due to the bending of the vortex tubes creates a significantly stronger increase. Dust clumping in 3D is up to almost an order of magnitude more efficient as compared to the 2D cases. In the 2D simulation with the same configuration as that of the 3D one, the enhancement factors are similar to those in other 2D simulations, reaching a maximum enhancement of ten times the initial density dust. This confirms that the increased maximum increase is inherently due to 3D effects.

Eventually, the twisting of the vortex tubes leads to a disruption of the vortex structures, as can be seen in Fig. 13. This is in contrast to the 2D case, where an evolution from small to large scales is typically seen. The formation of small-scale dust structures with enhanced densities can also clearly be seen in the power spectrum of the simulation. If we take a 3D FFT (fast Fourier transform) of the density distributions of the gas and dust fluids at t = 40 (the end of the simulation) and convert this to give the power per lengthscale as is shown in Fig. 15, we clearly see how the maximum in the spectrum shifts to smaller scales if we go from the gas fluid (wavenumber κ ~ 4, tracing the size of the vortices) to the smallest dust species (κ ~ 17) and to even smaller scales if we go to the largest dust species (κ ~ 24). Figure 16 gives a quantitative comparison of the strength of the 13 first modes in the middle plane as a function of time. It clearly shows the evolution to straight vortex tubes during the linear phase, followed by an early non-linear phase where there are no important 3D effects and ultimately the disruption to a turbulent phase in which the vortices break up into small-scale structures.

5. Discussion: filaments in molecular clouds

The simulations described above all have physical parameters, which are similar to values typically observed in molecular clouds. In this section, we highlight some similarities between the simulated KHI in 3D and structures observed in molecular clouds. Molecular clouds are regions in interstellar space where the local density is much higher than in the surrounding regions. Also temperatures are typically much lower. The mean diameter of a giant molecular cloud is 45 pc (Blitz 1993). Structure formation is observed on all sizes in molecular clouds from 0.003 to 30 pc (Genzel 1991). Typical substructures are called “clumps”, “sheets”, “bubbles” and “filaments”. We have seen in Sect. 4.3 that these types of structures are also formed in the 3D KHI (as can be seen in the left panel of Fig. 17). Due to the high densities and low temperatures, observations of molecular clouds are often made in infrared wavelengths, in which the emission traces the dust distribution. To compare our 3D simulation with observation, we make a synthetic observation by summing the total dust density along the line of sight in a certain direction, giving us the dust mass column density σ_d: ${\mathit{\sigma }}_{\mathit{d}}\mathrm{=}\mathrm{\int }{\mathit{\rho }}_{\mathit{d}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\kappa }}_{\mathit{\nu }}}\mathrm{\int }{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{\rho }}_{\mathit{d}}\hspace{0.17em}\mathrm{d}\mathit{s}\mathrm{=}\frac{{\mathit{\tau }}_{\mathit{\nu }}}{{\mathit{\kappa }}_{\mathit{\nu }}}\mathit{,}$(18)with κ_ν as the dust opacity and τ_ν the optical depth at frequency ν. If we assume that we look at the system in a wavelength at which it is optically thin (τ_ν ≪ 1), the radiation is caused by thermal dust emission and the dust temperature is uniform, then σ_d can indeed be used as an estimator for the observed flux of the system. The radiative transfer then simplifies to ${\mathit{I}}_{\mathit{\nu }}\mathrm{\approx }{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{\sigma }}_{\mathit{d}}\mathit{,}$(19)and thus the emitted intensity scales linearly with the mass column density. The σ_d is calculated in the right panel of Fig. 17 from the 3D data of the left panel. This total density along the line of sight can then be compared to column densities, which are derived from the observations. One of the important measurements derived from the column density is the probability density function (PDF), which is the probability distribution of the column density values in the observation. It has been shown that the PDF of isothermal structures with turbulence and negligible magnetic and gravitational forces is taking the shape of a log-normal distribution (Vazquez-Semadeni 1994). This result has also been confirmed from magnetohydrodynamics-simulations with supersonic turbulence and self-gravity, (Ostriker et al. 2001) and isothermal hydrodynamical simulations with supersonic turbulence without self-gravity (Federrath et al. 2010). This log-normal shape of the PDF is also seen in observation; however, often a power-law excess is seen on the side of the higher column densities (Kainulainen et al. 2009). Kainulainen et al. (2009) find evidence that links this possible excess to molecular clouds with active star-formation.

Fig. 17 Left: volume plot of a simulation of the KHI at t = 40. The heaviest dust species (of two species) is shown with the opacity in the plot scaled with the dust density. One can clearly see the elongated, filamentary structures with several high-density clumps where the density has increased to more than 10 times its value. In between, the structures are vacuum bubbles where the density has decreased by several orders of magnitude. Right: line of sight surface density of the heaviest dust species, as derived from the 3D simulation on the left. The density scale is in code units. Elongated structures with embedded high-density clumps are clearly seen.

Significant deviations from log-normal behavior is also seen in the low-density regions of many clouds. While the log-normal distribution is generally seen as a consequence of supersonic turbulence, Tassis et al. (2010) show that such distributions are obtained without supersonic turbulence in three very different numerical setups of molecular clouds. Furthermore, they claim that the development of power-law tails is not due to intermittency or gravity taking over the dynamics.

Fig. 18 Column density distribution of the synthetic observation of the KHI on the right side of Fig. 17 at t = 40, as shown in black. The column density is shown in code units. The central part of the distribution is fitted with a log-normal distribution, as shown as the dashed red line. A power-law tail is seen on the high-density side, as well as an additional excess on the low-density side, which is comparable to what is seen in some observations of molecular clouds.

If we make a PDF from our synthetic dust observation as shown in the right side of Fig. 17, we get the distribution shown in Fig. 18. We thus see that the PDF in our simulation of the KHI has evolved from one initial value of the column density to a distribution spanning two orders of magnitude. The central part of the PDF is fitted by a log-normal distribution; a power-law-like excess is seen on the high density side, as well as an excess on the low density side. This is very similar to some of the observations in Kainulainen et al. (2009). The excesses can be explained in the case of the KHI; values on the low density side correspond to the low-density regions created by the vortex formation in the KHI. The higher density regions are the structures of enhanced dust density formed first by the dust separation in the non-linear phase of the KHI, and later in the non-linear phase by the instability in the third dimension.

While the other mechanisms that we discussed above, which have been studied previously, may indeed be important for the formation of the structures observed in molecular clouds, it seems reasonable that the omniprecense of shear flows in astrophysical fluids can also contribute to the formation of the structures and to the formation of the log-normal distribution in the PDFs. This may explain some of the deviations from the log-normal distribution. In future work, we explore this further by performing dedicated simulations with properties derived from observations and link these with the dust radiative transfer code SKIRT (Baes et al. 2011) to be able to conduct a more detailed comparison with observations.

6. Conclusions

We have presented the wavelength dependence of the growth rate of the KHI with dust-to-gas ratios which range between δ = 0 and δ = 1. To do so, we performed 44 individual simulations with different values of δ and κ. We see that the growth rate reduced slightly (2%) for δ = 0.01. Adding more dust reduces the growth rate more significantly. The reduction is stronger for short wavelengths; for δ = 1, short wavelengths are completely stabilized. This leads to a shift of the most unstable wavelength to longer wavelengths. For the linear phase of the KHI, we have also investigated the dependence of the growth rate on the shearing velocity from the case δ = 0.1. We see that the difference with gas-only simulations is stronger for increasing shearing velocities. Shocks are seen in the gas when the velocity difference becomes supersonic; for the dust, these shocks are only seen for small grains, as the larger grains are more weakly coupled to the gas.

Fig. A.1 The summation of the densities for all the different dust species gives us the total dust density, as shown here for a simulation with two dust species (left panel), four dust species (central panel), and eight dust species (right panel). Only part of the larger domain of (1.58 × 2.4) × 1017 cm is shown. While the structures in the simulations with four and eight species are similar, the simulations with only two dust species display a less diffusive structure. The rarefaction rate is also underestimated in the simulation with two species, leading to higher minimum densities than in the two other cases. This is due to lower values of the largest represented dust species.

The formation of the vortices evacuates dust from the low-pressure vortices at a fast rate, with the minimum density decreasing exponentially for all dust species after the linear phase. High density dust layers are formed around the vortices and in the braids between vortices. Small dust particles form structures closer to the centre of the vortex. The formation of both low density and high density dust regions is stronger for increasing dust particle sizes in the range between 5−250 nm investigated here. The volume percentage of the regions, which have strongly increased dust densities, increases linearly after the end of the linear phase for the cases δ = 0.01 and δ = 0.1. For δ = 1, the trend is less clear. In 2D simulations, where we start with 20 vortices, merging leads to an increase in the size of the structures. In our 3D simulations, the vortex tubes become unstable in the additional third direction, which leads to a breakup of the tubes and eventually to the formation of small-scale structures. These structures are smaller for larger particles than for smaller particles or gas. The additional instability also increases the enhancement of dust densities even further with peaks in density enhancement up to seven times stronger than in comparable 2D simulations.

We see that a dusty KHI with physical values comparable in molecular clouds is able to form structures such as filaments, which are prevalent in molecular clouds, in relative short amounts of time (the end of the 3D simulation at t = 40 corresponds to 0.13 Myr). Furthermore, the KHI reshapes the PDF of the column densities from a single value in the initial uniform distribution to a central log-normal distribution with excesses on both sides. This is comparable to what is seen in observations without needing supersonic turbulence or self-gravity to form the log-normal distribution and the power-law tail.

Acknowledgments

We acknowledge financial support from the EC FP7/2007-2013 grant agreement SWIFF (No. 263340) and from project GOA/2009/009 (KU Leuven). This research has been funded by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office (IAP P7/08 CHARM). Part of the simulations used the infrastructure of the VSC − Flemish Supercomputer Center, funded by the Hercules Foundation and the Flemish Government − Department EWI. The authors would like to thank the anonymous referee for the valuable comments.

References

Appendix A: Method validation

To test the effectiveness of the dust treatment described in Sect. 3.2, we compare between three simulations with the same physical parameters (see Table 1 with δ = 0.01 and a total domain of (1.58 × 2.4) × 10¹⁷ cm), however with different amounts of dust fluids; namely two, four and eight different species with sizes between 5 nm and 250 nm. The effective sizes for each fluid, calculated using Eqs. (10) and (11), are given in Table A.1. For all three setups, simulations are conducted up to t = 60, where non-linear effects, such as vortex merging, are observed.

Most importantly, no large differences are seen between the behavior of these three simulations. Their linear phase is equally long, and in the non-linear phase, the vortices merge at the same rate and form structures of comparable sizes. However, small differences can be seen: if we look at the dust distribution in the vortices (see Fig. A.1), we see that the total dust density is more spread out for the simulation with four and eight species compared to the other setup. In simulations with more fluids, the outer bins have smaller/larger representative particle radii, leading to more extreme behavior. As different dust sizes are evacuated from the vortex at different rates, this allows simulations with more fluids to have more distinct behaviors, leading to more diffuse distributions of the total dust density. The minimum values are also slightly lower for simulations with more fluids, an effect of having larger values for the maximum effective grain radius, as discussed in Sect. 4.2.1. More quantitatively, we compare the evolution of the transverse kinetic energy of the three simulations. An example of the typical evolution of the transverse kinetic energy is shown in Fig. 3. Here, we look at the relative error of the simulations, as compared to the setup with eight species (Fig. A.2). For the simulation with four species, the difference is less than 2% until t ≈ 40, which is long after the end of the linear phase (t = 10). Even in the non-linear phase differences are never bigger than 50%. The differences for the simulation with two species are larger, with a difference of ~20% at the beginning of the simulation and a relative difference up to three times the value of the simulation with eight species in the non-linear phase. Nonetheless, differences between slightly different setups are expected due to the non-linear nature later in the simulations.

Fig. A.2 Relative error of the transverse kinetic energy in simulations with two and four dust species as compared to the simulation with eight species.

While we can conclude that the overall physics is captured in all these simulations, we note that minor small-scale differences are observed in the dust distribution. While the growth rate of the linear phase is comparable for all three simulations, differences in evolution during the non-linear phase are also observed between the simulation with two species on the one hand and the simulations with four and eight species on the other hand. Choosing for a setup with four dust species is therefore a good choice with the treatment that we described in Sect. 3.2.

All Tables

All Figures

	Fig. 1 Full domain of a 2D simulation at t = 40. The background colour represents the refinement level with lighter colours being higher levels of refinement. While seven levels are used in the simulation, the coarsed level is no longer used at this time. In the middle, the gas density is overlaid to illustrate how the AMR chooses maximal refinement around the vortices.
In the text

	Fig. 2 Evolution of the ratio between the actual amount of cells in a simulation and the amount needed to fill a uniform grid with the same resolution. The difference between the 2D and 3D setup are due to differences in the extent of the simulated domain.
In the text

Fig. 3 Evolution of the transverse kinetic energy of the gas in six different simulations with the same dimensionless wavenumber κ = 0.7968 (the most unstable wavenumber in the gas-only case). The vertical axis has been set to logarithmic scale to clearly show the exponential growth in the initial linear growth phase. Three simulations start with a lower initial transverse kinetic energy. This is because they start with small random perturbations in the middle layer, instead of being excited with the perturbation described in Eq. (8). These three simulations are discussed in the Sect. 4.3. The behavior of the simulation with δ = 1 is different; this is explained in Sect. 4.1.1.

In the text

Fig. 4 Comparison of the values of vy of a simulation (+ symbols) and the derived eigenfunctions (dotted lines) on a cut perpendicular to the flow direction. For clarity, only a limited amount of simulated points is shown, The values of vy are given in cm s-1. The cut, taken at t = 3, only spans a small portion of the range along the y-axis, and is centred on the middle of the boundary layer. The distance from the centre is given in units of D. The dotted lines each represent a solution of the eigenfunction, the red line in the upper flow, teal in the middle layer, and blue in the lower flow. From the correspondence with the simulated values, it is clear that during the perpendicular velocity the linear phase evolves from the perturbation in Eq. (8) to the distribution expected from the analytical derivation of the eigenfunction.

In the text

	Fig. 5 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 Sect. 2.2.
In the text

Fig. 6 Comparison of the growth rates of simulations with different fluid velocities. The Mach number indicated is the Mach number of the velocity difference of the upper and lower flow. Fourth order spline fits have been added to guide the eye. Dust tends to stabilize the KHI for all velocities with an increasing effect at higher velocities. The growth rate increases until the velocity difference Δv becomes supersonic.

In the text

	Fig. 7 Vortices that have formed shock structures at t = 35 in a simulation with Δv0 = 1.51vs, κ = 0.7968 and δ = 0.1. Shocks are clearly visible for the gas (left panel) and the smallest dust species (middle panel). In the second smallest dust species (right panel), the shocks can still be discerned; however, they are less evident. In larger dust species, these shocks are not seen.
In the text

Fig. 8 Density distribution of the largest dust grains out of four species at t = 10 in a 2D simulation with δ = 0.01 and κ = 0.7968. Only part of the simulated domain is shown. In the centre of the middle vortex, the black ellipse is a streamline of the gas fluid. The gas trajectory is (almost) closed as the gas fluid is pressure supported by the lower pressure in the centre of the vortex. By contrast, the streamlines of the large dust particles (shown in red) can be seen to spiral out of the vortex. Ultimately, the dust ends up in two layers below and above the vortex structures.

In the text

Fig. 9 Left: evolution of the global minimum in dust density with time for all 4 dust fluids in a simulation with δ = 0.01 and κ = 0.7968. The size of the assumed particles in the fluid increases from dust 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. Only the part of the simulation is shown, where all minima are above the allowed minimum density (see Sect. 2.1). Right: global maximum for the four dust species. Density increases are stronger for heavier dust species.

In the text

Fig. 10 Density distribution of the gas (top left), the lightest dust species (top right), and the heaviest dust species (bottom left), as well as the vorticity of the gas fluid (bottom right) at t = 60 in a simulation with δ = 0.01 and κ = 0.7968. Only part of the total domain is shown. In the image, the merging of two vortices is ongoing. A clear connection can be seen between the vorticity of the gas and the distribution of the smallest dust species.

In the text

Fig. 11 Volume percentage of the simulated domain where the dust density has at least increased to 2.5 times the initial dust density. All simulations have v0 = 0.42 and are perturbed with wavelengths close to their most unstable value (see Fig. 5). The portion of the domain where the dust density is significantly enhanced is seen to increase linearly with time after the linear phase (t ~ 10 for δ = 0.01 and δ = 0.1 and t ~ 20 for δ ~ 1.0). The rate of increase decreases with dust-to-gas ratio δ.

In the text

Fig. 12 Volume rendering after the end of the linear phase (t = 22.5). The entire domain of the simulation is shown. The volume rendering shows the gas distribution of the two straight vortex-tubes that have formed. To visualize the vortices in the volume rendering, we only show gas densities below 9.2 × 10-21 g cm-3. At the back and left side of the cube (xy- and yz-planes, respectively), the dust distribution of the heaviest dust species is shown. The dust forms tubes around the gas vortices, and sheets of increased density between the tubes.

In the text

	Fig. 13 Volume rendering at t = 40. To visualize the vortices in the volume rendering, we only show gas densities below 8.6 × 10-21 g cm-3. The straight vortex tubes seen at the end of the linear phase in Fig. 12 have become severely twisted and are breaking apart. The dust distribution at the same time is shown in Fig. 17.
In the text

Fig. 14 Comparison of the dust enhancement in four different simulations with different values of δ. The graph shows the maximum of the dust density in the simulation, ascompared to the initial value of the dust density. Dust enhancements are (initially) slightly higher for lower values of δ. The enhancements in the 3D case are significantly stronger due to the additional instability of the vortex tubes in the added third dimension. A log-scale is used for the maximum density. Note that the maximum density starts to increase at a later time in the 3D simulation as compared to the 2D cases because a random initial perturbation is used in 3D.

In the text

Fig. 15 Power per lengthscale derived from the 3D Fourier transform of the densities of the three fluids at t = 40. Formed structures are smaller for the dust species, with dust 1 being the dust species representing the smaller dust particles, and dust 2 representing the larger particles. The ripples at high wavenumber for the gas are due to the AMR, as we choose the active refinement to follow the density variations of the two dust fluids.

In the text

Fig. 16 Comparison of the strength of each mode of vy along the z-direction of the middle plane. Modes 0 to 13 are shown. At each time, the strength of each mode is compared to the strength of the dominant mode at that time. Three distinct regimes can be discerned in the figure. Initially, the initial random perturbations evolve up to t ~ 8 from small-scale structures to a large scale feature. In the linear phase and the early non-linear phase (between t ~ 8 and t ~ 30), the vortex tubes are elongated along the z-direction and the 0-mode is dominant. After t ~ 30, the vortex tubes become unstable, resulting in the formation of small-scale structures.

In the text

Fig. 17 Left: volume plot of a simulation of the KHI at t = 40. The heaviest dust species (of two species) is shown with the opacity in the plot scaled with the dust density. One can clearly see the elongated, filamentary structures with several high-density clumps where the density has increased to more than 10 times its value. In between, the structures are vacuum bubbles where the density has decreased by several orders of magnitude. Right: line of sight surface density of the heaviest dust species, as derived from the 3D simulation on the left. The density scale is in code units. Elongated structures with embedded high-density clumps are clearly seen.

In the text

Fig. 18 Column density distribution of the synthetic observation of the KHI on the right side of Fig. 17 at t = 40, as shown in black. The column density is shown in code units. The central part of the distribution is fitted with a log-normal distribution, as shown as the dashed red line. A power-law tail is seen on the high-density side, as well as an additional excess on the low-density side, which is comparable to what is seen in some observations of molecular clouds.

In the text

Fig. A.1 The summation of the densities for all the different dust species gives us the total dust density, as shown here for a simulation with two dust species (left panel), four dust species (central panel), and eight dust species (right panel). Only part of the larger domain of (1.58 × 2.4) × 1017 cm is shown. While the structures in the simulations with four and eight species are similar, the simulations with only two dust species display a less diffusive structure. The rarefaction rate is also underestimated in the simulation with two species, leading to higher minimum densities than in the two other cases. This is due to lower values of the largest represented dust species.

In the text

	Fig. A.2 Relative error of the transverse kinetic energy in simulations with two and four dust species as compared to the simulation with eight species.
In the text