Formation of Dust Filaments in the Diffuse Envelopes of Molecular Clouds

Molecular clouds (MCs) are the galactic reservoirs of molecular material where star formation takes place. It is therefore necessary to understand their formation and evolution if one aspires to comprehend the physical processes involved in star formation.

It is known that MCs are formed from a more diffuse medium that becomes thermally unstable (see the recent reviews by, e.g., Ballesteros-Paredes et al. 2020 and Chevance et al. 2020). Their turbulent nature (Elmegreen & Scalo 2004) is inherited from the parental gas (Meidt et al. 2018; Kruijssen et al. 2019), and is maintained by several processes both at large and small scales, such as supernova explosions (Padoan et al. 2016; Bacchini et al. 2020) or stellar jets (Federrath 2015; Offner & Liu 2018). This turbulence, together with the pervasive galactic magnetic field, plays a fundamental role in the fragmentation of the cloud and the formation of dense gas filaments where stars are actually formed (Myers & Goodman 1988; Fiege & Pudritz 2000; Nakamura & Li 2008; Inutsuka et al. 2015).

MCs are formed from the accretion of warm, diffuse gas from the warm neutral medium (WNM), characterized by its atomic composition (mainly H I), temperatures of approximately 6000 K, densities ranging from 0.2 to 0.5 cm⁻³ (Ferrière 2001) and ionization fractions of χ ∼ 0.1 (Klessen & Glover 2016); these conditions are sustained in time due to the constant bombardment of the interstellar radiation field over the molecular gas. Therefore, MC envelopes share properties with the WNM, and they constitute the transition layers between the hot ionized gas of the diffuse interstellar medium (ISM) and the interiors of the clouds. The net charge of the species in the outskirts of these clouds favors their coupling with the galactic magnetic field, which acts as a carrier of hydromagnetic waves (Balsara 1996); they may also facilitate the penetration of cosmic rays inside the cloud (Ivlev et al. 2018).

A key ingredient in the transition from atomic to molecular gas is interstellar dust, for which the distribution seems to be correlated with that of the gas (Bohlin et al. 1978; Boulanger et al. 1996; Roman-Duval et al. 2010; Lenz et al. 2017). Despite representing only 1% of the mass of the ISM, interstellar dust grains play a fundamental role in several physical phenomena. They act as catalysts in chemical reactions (Hollenbach & Salpeter 1971; Cazaux & Tielens 2004; Shingledecker et al. 2020), shield the interior of the clouds from the high energy ultraviolet (UV) radiation, and contribute to the thermal regulation of the cloud (Bakes & Tielens 1994; Wolfire et al. 1995; Weingartner & Draine 2001; Hocuk & Spaans 2010).

Due to their composition, dust grains acquire a net charge in the presence of a radiation field (Weingartner & Draine 2001); even in the densest regions of MCs, grains may be charged due to the action of cosmic rays (Ivlev et al. 2015). In MC envelopes, dust grains contribute to the magnetic field coupling and affect the propagation of hydromagnetic waves (Pilipp et al. 1987; Cramer & Vladimirov 1997) and, in the last instance, star formation (Nakano 1998). It is then crucial to consider interstellar dust grains when studying the formation and evolution of MCs.

Over the last decade, our knowledge of the cold ISM has grown exponentially thanks to the results of dedicated missions such as Herschel (see, e.g., André et al. 2010; Arzoumanian et al. 2011; Palmeirim et al. 2013) and Planck (Planck Collaboration et al. 2020a, 2020b). However, the characterization of the diffuse ISM is more challenging, since it requires tracers sensitive to low column densities, such as ¹²CO (Falgarone & Phillips 1996; Falgarone et al. 2005) or very high-resolution H I maps (Clark et al. 2014; Martin et al. 2015). An alternative approach, which is very time-consuming but provides a higher resolution, is to perform a spectroscopic survey of the neutral and ionized components of the ISM at optical and UV wavelengths (Savage & Sembach 1991; Lallement et al. 2003; Redfield et al. 2004; Welsh et al. 2010).

With these ingredients, it seems that a full characterization of MCs may be within reach. However, observations give at most a line-of-sight averaged view of the processes involved in MC evolution and thus, to understand the underlying physics, numerical simulations of these processes are required.

It has been predicted that negatively charged grains may affect the star formation efficiency due to gyroresonance effects with MHD waves propagating inside MCs (Pilipp et al. 1987; Nakano 1998; Falceta-Gonçalves et al. 2003; Hopkins et al. 2020). For instance, Chapman & Wardle (2006) demonstrated that C-type MHD shock profiles are deeply affected by the presence of a population of charged grains; Pandey & Vladimirov (2019) have recently shown that the Kelvin–Helmholtz instability in MCs depends on the dust charge density and size distribution; and Wardle (2007) and Pinto et al. (2008) derived the modified coefficients for nonideal MHD effects (ohmic resistivity, and ambipolar and Hall diffusion) in the presence of charged dust grains. However, none of the abovementioned studies followed the numerical evolution of an MC in the presence of a population of dust.

The main trend has been to consider the aerodynamic drag of dust in the framework of protoplanetary disks (Bai & Stone 2010; Price et al. 2018; Mignone et al. 2019). To our knowledge, the only authors that have included the physics of charged grains in their MHD simulations are Lee et al. (2017), who have devoted a series of papers to the exhaustive study of resonant instabilities that arise between gas and charged dust over different phases of the ISM when dust motions are driven by an external force (Hopkins & Squire 2018; Seligman et al. 2019; Hopkins et al. 2020).

In this article, we study the dynamics of both gas and charged dust under conditions typical of a molecular cloud envelope with a modified version of the MHD code Athena ³ (Stone et al. 2008). In Section 2 we present our implementation of the dynamics of charged dust. In Section 3, we present the numerical 2D model for the MC envelope (Section 3.1) and analyze the local variations of the dust-to-gas ratio (Section 3.2) as well as the dust filamentary structure (Section 3.3). In Section 4, a discussion of the results is provided and the main results are highlighted in Section 5.

We worked with the latest stable version of Athena (v4.2). This C-version has already implemented an aerodynamic particles module (Bai & Stone 2010, hereafter BS2010). In this section, we detail how we have expanded this module by adding the necessary equations that take into account the dynamics of charged dust. We have also carried out different tests for checking this new implementation that are detailed in the Appendix.

Apart from the traditional drag term characterized by the stopping time t_s (Weidenschilling 1977), we added two terms to the equation of motion of a charged particle that account for the Coulomb drag from charged gas species, and the Lorentz force:

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=-\displaystyle \frac{{\boldsymbol{v}}-{\boldsymbol{u}}}{{t}_{s}}-{\nu }_{C}({\boldsymbol{v}}-{\boldsymbol{u}})+\displaystyle \frac{{Z}_{d}e}{{m}_{d}c}({\boldsymbol{v}}-{\boldsymbol{u}})\times {\boldsymbol{B}}.\end{eqnarray} \tag{ 1 }$$

In the above equation, v is the velocity of the dust grain, u is the gas velocity, t_s is the aerodynamic stopping time, Z_d is the grain charge, m_d is the dust mass, c is the speed of light, and ν_C is an additional drag term that accounts for Coulomb interactions between charged grains and ions/electrons in the plasma:

$$\begin{eqnarray}&&{\nu }_{C}=\displaystyle \frac{4}{3}\displaystyle \frac{\sqrt{2\pi }{Z}_{d}^{2}{e}^{4}\mathrm{ln}{\rm{\Lambda }}{n}_{d}}{{({kT})}^{3/2}}\left(\displaystyle \frac{\delta }{\sqrt{{m}_{\mathrm{ion}}}}+\displaystyle \frac{(1-\delta )}{\sqrt{{m}_{e}}}\right).\end{eqnarray} \tag{ 2 }$$

In this equation, $\mathrm{ln}{\rm{\Lambda }}=\mathrm{ln}[{(3{kT}/2{e}^{3})({{kTm}}_{e}/{\rho }_{e}\pi )}^{1/2}]$ is the Coulomb logarithm, n_d is the density of charged dust grains, and δ = 1 if Z < 0, δ = 0 otherwise; this means that negatively (positively) charged dust grains interact with the ions (electrons) in the gas. We disregarded the effect of the electric field derived from the magnetic field in the Lorentz force because we are mainly interested in studying the effects of first-order terms.

We coupled the integration of the Lorentz force and the Coulomb drag in the semi-implicit integrator of BS2010, which is very similar to the Boris integrator (Boris 1970) implemented by Lehe et al. (2009). Following the notation by BS2010, the equations to solve for the particle dynamics are:

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\boldsymbol{v}},\,\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}={\boldsymbol{a}}[{\boldsymbol{v}},{\boldsymbol{x}},{{\boldsymbol{u}}}^{n+1/2}({\boldsymbol{x}})],\end{eqnarray} \tag{ 3 }$$

where we have grouped all the acceleration terms in Equation (1) under the term a . Then, the velocity update can be computed using the predicted particle position x ' as:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}^{(n+1)} & = & {{\boldsymbol{v}}}^{(n)}+h{{\rm{\Lambda }}}^{-1}{\boldsymbol{a}}({{\boldsymbol{v}}}^{(n)},{\boldsymbol{x}}^{\prime} ,{\boldsymbol{u}}({\boldsymbol{x}}^{\prime} )),\\ {\rm{\Lambda }} & = & 1-\displaystyle \frac{h}{2}\displaystyle \frac{\partial {\boldsymbol{a}}}{\partial {\boldsymbol{v}}}\end{array}\end{eqnarray} \tag{ 4 }$$

where the Jacobian matrix for Equation (1) is:

$$\begin{eqnarray}\displaystyle \frac{\partial {\boldsymbol{a}}}{\partial {\boldsymbol{v}}}=\left(\begin{array}{ccc}-\tfrac{1}{{t}_{s}}-{\nu }_{C} & {Q}_{m}{B}_{z} & -{Q}_{m}{B}_{y}\\ -{Q}_{m}{B}_{z} & -\tfrac{1}{{t}_{s}}-{\nu }_{C} & {Q}_{m}{B}_{x}\\ {Q}_{m}{B}_{y} & -{Q}_{m}{B}_{x} & -\tfrac{1}{{t}_{s}}-{\nu }_{C}\end{array}\right).\end{eqnarray} \tag{ 5 }$$

In the above matrix, we have renamed the term Z_d e/m_d c as Q_m for simplicity.

The last consideration for the implementation of charged particles is the numerical timestep. In order to integrate accurately the oscillation of a charged particle, the timestep has to be considerably smaller than the particle's Larmor time. As a compromise between numerical accuracy and computational cost, we have chosen the same value as that by Lee et al. (2017):

$$\begin{eqnarray}&&{\rm{\Delta }}t=\displaystyle \frac{\mathrm{CFL}}{10}\min \left(\displaystyle \frac{{m}_{d}c}{{Z}_{d}e| {\boldsymbol{B}}| },{\rm{\Delta }}{t}_{\mathrm{MHD}}\right),\end{eqnarray} \tag{ 6 }$$

where CFL is the Courant–Friedrichs–Lewy number and t_MHD is the timestep computed without considering the dust particles.

3.1. Numerical Model

We solve with Athena the equations for isothermal, ideal MHD in 2D together with Equation (1) for the dust dynamics in a box of length L with periodic boundary conditions.

As a first step in the study of the dynamics of charged dust in MCs, we developed a very simple model of an MC envelope that responds to the common physical parameters of the WNM. The ambient gas is composed by H I at T = 6000 K with an ionization fraction χ = 0.1; the magnetic field strength is chosen to be B = 10⁻⁶ G, in concordance with the measured values in the ISM of a few μm (Ferrière 2001). The gas density is n = 10 cm⁻³, a selection justified by observational evidence (Blitz et al. 2007; Fukui et al. 2009, 2017).

The dust population is represented by 80,000 test particles that share the same properties: a radius of a_d  = 5 × 10⁻⁶ cm and an internal density of ${\rho }_{d}^{\mathrm{int}}=1\,{\rm{g}}$ cm⁻³. They are randomly distributed over the whole domain with a velocity taken from a Gaussian distribution with deviation v_th, the isothermal sound speed. For these grains, we computed an average electric charge under the assumption of statistical equilibrium of Z = −17 following the procedure described in Weingartner & Draine (2001); the intensity of the interstellar radiation field (Mathis et al. 1983) is considered to be G = 1.

Initially, the medium is uniform in density, with the magnetic field parallel to the x-axis. The box size is chosen to be L = 1 pc, a value that allows us to resolve the apparent characteristic molecular gas filament widths of 0.1 pc (Arzoumanian et al. 2011) and, of more relevance for the present study, the width of H I fibers observed in the diffuse ISM (sizes below 0.04 pc, Clark et al. 2014): at the chosen resolution (1024 × 1024, see below), one pixel size corresponds to 1/1024 ∼ 9.76 × 10⁻⁴ pc. We introduced a turbulent spectrum in velocity as the superposition of four waves with wavelengths λ₀ = L/24, λ₁ = 26L/72, λ₂ = 49L/72, and λ₃ = L and wave amplitude u₀ = v_th(k_n /k₀)^3/2:

$$\begin{eqnarray}\begin{array}{rcl}{u}_{x} & = & {u}_{0}\displaystyle \sum _{n=0}^{3}\cos ({k}_{n}y)\\ {u}_{y} & = & {u}_{0}\displaystyle \sum _{n=0}^{3}\sin ({k}_{n}x)\\ {u}_{z} & = & 0.\end{array}\end{eqnarray} \tag{ 7 }$$

The wave spectrum has been chosen such that the number of waves (N_waves = 4) is large enough to produce a turbulent state in the gas, but small enough to follow the individual effects of each wave. The longest wavelength is chosen to be λ_max = L, while the shortest one is λ_min = 1/24, much larger than the resolution of the domain; the other wavelengths have been chosen from the following relationship: λ_n  = λ_min + n ∗ (λ_max − λ_min)/(N_waves − 1).

The simulation has been run with a resolution of 1024 × 1024 pixels² and up to t = 5 (code units, ∼0.4 Myr). The reason behind this choice is that at t = 5, the initial perturbations have been damped by a factor of 2 and they are bounded for t > 5, so we consider that the system has been stabilized.

3.2. Decoupled Motions of Gas and Dust

Motivated by the apparent linear relationship between gas and dust column density reported by several authors (Savage et al. 1977; Bohlin et al. 1978; Predehl & Schmitt 1995; Lenz et al. 2017), we studied the distribution of gas and dust in the simulation.

We built dust distribution maps from the particle's positions assuming an initial dust-to-gas ratio ρ_d /ρ_g  = 0.01 (Spitzer 1954). For the initial conditions, we estimated the dust mass contained in the total volume of 1 pc² (×1 pc) and distributed it equally over the 80,000 test particles. We considered that the test particles trace the distribution of real particles, too large to be computationally resolved, and that those real particles are uniformly distributed around the test ones. This post-processing is done after the main simulation, so it does not affect the formation of dust filaments. We distributed the dust mass over a matrix with the same size as the resolution of the simulation (1024 × 1024) and normalized it such that a value of 1 corresponds to the assumed dust-to-gas ratio of 0.01. The dust and gas maps at the final stage of the simulation are shown in Figure 1.

Zoom In Zoom Out Reset image size

If we look at the dust-to-gas ratio maps in Figure 2, we observe that there are certain regions where the values are significantly larger than the initial average of 0.01 (value equal to unity in the normalized maps). On the other hand, there are some regions of the simulation where the dust-to-gas ratio is very small or that are completely devoid of dust. This is an expected result of the submission of charged dust to a spectrum of Alfvén waves. Finally, it is worth noting that the dust overdensities show a filamentary structure that will be studied in Section 3.3.

Zoom In Zoom Out Reset image size

The decoupling of gas and dust could arise due to the large differences between the particle stopping time t_s and the Larmor time, t_L  = m_d c/∣Z_d ∣e∣ B ∣. For all particles in the simulation at the final stage, the ratio t_s /t_L is always greater than one (see Figure 3) indicating that the coupling with the magnetic field is dominant. We selected some regions in the gas density map in order to study in detail the relationship between the distribution of dust, the gas density and the magnetic field strength. The regions are displayed in Figure 2 and were selected such that a region of high (low) gas density shows values greater (lower) than 1σ, since the gas distribution can be approximated by a Gaussian.

Zoom In Zoom Out Reset image size

The analysis of the selected regions shown in Figure 2 may be summarized as follows:

• 1.  
  At high gas density regions (Figure 4) we observe large values of the dust-to-gas ratio forming filaments, parallel to the magnetic field.
• 2.  
  Except for region H3, where the dust density is extremely poor, the dust mean velocity is similar to the Alfvén velocity, pointing toward a strong coupling of dust with the magnetic field.
• 3.  
  At low gas density regions (Figure 5), the dust filamentary structure is still observed, but there is more granulation and the filament widths are narrower. The modulus of the magnetic field in these regions is approximately 2.5 times greater than in regions H1-H4, and the dust velocity is again similar (and for regions L1 and L4 almost equal) to the Alfvén velocity.
• 4.  
  Finally, for the mid gas density region (Figure 6) we observe a hybrid behavior. Dust filaments are also formed parallel to the magnetic field, but the overdensities are broader and, generally, perpendicular to the field. In this case, the mean dust motion seems more coupled to the gas dynamics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Filamentary Dust Structure

In order to better understand the properties of the dust filaments generally decoupled from the gas (see Section 3.2), we analyzed in detail their properties using FilFinder (Koch & Rosolowsky 2015) over the dust mass map at the final stage of the simulation. For the analysis, we created a custom mask preprocessing the map with the python packages scikit-image (van der Walt et al. 2014) and scikit-learn (Pedregosa et al. 2011). First, we binarized the map considering only pixels with a dust mass greater than the median; then, we identified connected regions and removed those with fewer than 10 pixels. The resulting mask is then passed as an argument to FilFinder together with the original dust map; we then skipped the flattening and creation of a mask steps and performed directly the identification and pruning of filaments; the results are shown in Figure 7.

Zoom In Zoom Out Reset image size

We found it useful to study the filamentary structure over a degraded version of the map (256 × 256 px²) as well as over the original one (1024 × 1024 px²). While the latter represents well the fragmented structure of the dust clumps, the former gives an overall picture of the morphology of the dust filaments. The typical length of the filaments in the original image is 0.04 pc, while for the degraded image is 0.21 pc. We also computed the typical widths of the filaments as the mean width perpendicular to the longest path and found a value of 0.0035 pc for the high-resolution image and 0.0087 pc for the degraded one; note that the image size imposes a resolution limit of 9.8 × 10⁻⁴ pc for the 1024 × 1024 px² map, and of 3.9 × 10⁻³ pc for the 256 × 256 px² one, but both mean widths are above that threshold. Table 1 summarizes all the parameters of the analyzed regions. It is evident that, although the resolution of the dust map critically affects the filament length, the widths are similar and an order of magnitude smaller than the typical width of 0.01 pc for gas filaments in molecular clouds found by Herschel (Arzoumanian et al. 2011; but see Panopoulou et al. 2017 for a recent discussion on the validity of such a characteristic width); since at large scales the correlation between dust and gas still holds, our results are in concordance with the observational report by Clark et al. (2014) on diffuse H I filaments at all resolved scales, up to a lower limit of 0.04 pc imposed by the instrumental resolution. However, the dust filament width should not depend on the resolution of the global simulation, because dust dynamics is resolved based on the gas properties on the adjacent cells. Hence, as dust grains are coupled to the magnetic field, an increase of the resolution only contributes to a better definition of the shape of the gas filaments, but the overall dynamics of dust grains should remain unchanged. What does affect the morphology of the dust filaments is the wave spectrum, as is discussed later in the following section. The dust structures will be more clumpy or filamentary depending on the propagating waves. In that sense, the box size may play a role since it sets an upper limit on the longest wavelength.

Table 1. Properties of the Analyzed Regions

Region	Dimensions	〈DGR〉	〈B〉	〈vA 〉	〈vgas〉	〈vdust〉	w1024	w256
	pc × pc		10−6 G	105 cm s−1	105 cm s−1	105 cm s−1	10−3 pc	10−3 pc
H1	0.145 × 0.138	1.403	0.737	0.476	2.997	0.205	4.385	10.449
H2	0.121 × 0.089	0.798	0.829	0.534	2.402	0.135	2.981	9.975
H3	0.128 × 0.133	0.280	1.262	0.806	3.069	0.361	3.818	6.331
H4	0.131 × 0.163	1.090	1.015	0.649	2.629	0.558	4.109	9.827
L1	0.233 × 0.103	1.597	2.703	2.124	6.293	0.479	3.046	8.523
L2	0.443 × 0.109	1.767	2.437	1.841	5.935	0.597	3.300	11.728
L3	0.087 × 0.032	1.921	1.862	1.368	1.982	0.162	4.112	8.464
L4	0.238 × 0.140	1.459	2.514	1.910	5.010	0.280	3.362	10.919
M1	0.062 × 0.060	2.514	0.777	0.532	6.104	1.534	4.939	17.578

Note. Mean quantities for the dust-to-gas ratio (denoted by DGR), magnetic field, Alfvén velocity, gas velocity, and dust velocity. Also included are the mean dust filament widths on each region for the original map (w₁₀₂₄) and for a degraded version (w₂₅₆).

Download table as:  ASCIITypeset image

We have also studied the mean dust-to-gas ratio along the longest paths of the filaments, and the results are shown in Figure 8. We find that, in general, the ratio reaches values six times greater than the average value for the ISM and spans from ∼3 to ∼17, another indication that dust filaments are not tightly correlated with gas filaments. However, if those mean values are normalized by the filament's length (number of pixels, L_pix, Figure 8, lower panel) one can see that the dust-to-gas ratio inside each filament is approximately constant and proportional to its length by a factor that oscillates between 0.1 and 0.3. There does not seem to be any apparent relationship between the dust-to-gas ratio of a dust filament and its width, which suggests that such a width may be characteristic of the dust filaments. However, to confirm that hypothesis, more realistic simulations that consider a whole representative dust population are needed.

Zoom In Zoom Out Reset image size

With respect to the relative orientation between the dust filaments and the magnetic field, we find that in general dust filaments tend to be aligned with the magnetic field, see Figure 9; in the high-resolution map, ∼15% of the filaments have relative orientations between ±0.05 rad (parallel), while the vast majority (98%) have relative angles inside the range ±π/4 (dotted lines in Figure 9). These angles have been computed along the skeleton of each dust filament as follows: every three pixels of the skeleton, the dot product between the magnetic field vector at the intermediate pixel and the skeleton vector is used to retrieve the relative angle between them, and is restricted to be inside the range ([−π/2, π/2)).

Zoom In Zoom Out Reset image size

It is a common practice to assume a standard relationship between the dust and gas abundance. However, local variations of the dust-to-gas ratio are to be expected (Padoan et al. 2006; Hopkins 2014), especially when the dust particles are charged (Lazarian & Yan 2002; Lee et al. 2017), and they are actually found in the diffuse medium and even inside a cloud (Chen et al. 2015; Reach et al. 2015; Siebenmorgen et al. 2017). There were some early attempts to determine the scale of the fluctuations from extinction measurements (Thoraval et al. 1997, 1999) and variations of the dust size distribution have been found in MCs (Flagey et al. 2009; Gómez de Castro et al. 2015; Köhler et al. 2015; Beitia-Antero & Gómez de Castro 2017). However, most of the studies are focused on the cool interiors of MCs since the main interest is to determine the possible influence of this phenomenon over the star formation processes. But if one aims to understand how dust may affect star formation and the formation of molecular clouds themselves, it is necessary to begin analyzing its behavior in MC envelopes, where the coupling to the magnetic field favors the propagation of hydromagnetic waves into the cloud and increase the turbulent state of the molecular gas.

Despite the simplicity of our model and having considered an initial homogeneous distribution of dust and gas, we have shown that the propagation of velocity waves produce local variations of the dust-to-gas ratio, forming narrow dust filaments aligned with the magnetic field. Due to the random distribution of dust grains at the beginning of the simulation, the ratio is not constant even at the initial stage of the simulation, although the mean scarcely varies with time at large scales, as shown in Figure 10. This fact could support the argument that, at large scales, the dust-to-gas ratio is nearly constant. However, at small scales there are appreciable deviations of the ratio that are far from a linear relationship. In Figure 11, the histogram of the dust-to-gas ratio map presented in Figure 2 is shown; there are clear variations of the ratio, and the peak of the distribution is far from the assumed dust-to-gas ratio (dashed vertical line).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recently, Lee et al. (2017) reported deviations of the dust-to-gas ratio in the cold neutral medium and Hopkins & Squire (2018) studied in detail the spectrum of instabilities that may arise when charged dust moves in a magnetized gas. In particular, they report that for the WNM the relative velocity between dust particles and gas should be subsonic and that the dominant resonances tend to form structures parallel to the magnetic field. Our results are compatible up to a point, since the majority of the structures are indeed aligned with the magnetic field and the mean dust and gas velocities differ by a few tens of km s⁻¹, see Table 1. These results are also in good agreement with observational evidence of dust filaments with low column densities aligned with the magnetic field (Planck Collaboration et al. 2016; Fissel et al. 2019). There seems to be only an exception with the mid-density region (see Figure 6), where a dust structure is formed perpendicular to the magnetic field direction. We attribute this localized phenomenon to the clear decrease of the magnetic field strength along the structure.

Additionally, we studied separately the influence of the largest and shortest waves of the initial velocity spectrum in the formation of dust structures. As can be seen in Figure 12, the shortest wave does not inject enough energy into the medium to produce conspicuous filamentary structures. Nevertheless, tiny clumps of a few (≤10) pixels are formed, and the histogram of the dust-to-gas ratio is very similar to the one for the complete simulation shown in Figure 11. For the longest wave, the energy input is so strong that relevant fluctuations are appreciable in density, momentum, and magnetic field distributions. These variations enable the formation of dust filaments with a dust-to-gas ratio distribution that again resembles to the one in Figure 11. We draw two conclusions from this analysis: first, that variations in the dust-to-gas ratio are expected to arise whenever a wave propagates in a molecular cloud envelope, no matter the energy of the wave; and second, that the morphology of the structures that can be formed in such a medium depends on the wave spectrum. It is of interest to explore in a future work if other properties of the dust filaments, such as their widths, depend on the wave spectrum or if they have a characteristic size. Finally, we want to highlight again that the choice of the box size sets an upper limit to the longest wavelength that we are able to resolve (L₀), so we are always talking about structures at scales lower than L₀ = 1 pc.

Zoom In Zoom Out Reset image size

4.1. Limitations of the Model

In this work, we have described a new module developed for Athena to deal with the interaction of charged dust particles in astronomical magnetized fluids. The code is applied to the study of dust dynamics in diffuse molecular cloud envelopes, and its performance is illustrated through the formation of dust filaments decoupled from the gas. The model for a molecular cloud envelope considers a density of n = 10 cm⁻³, with a temperature of T = 6000 K and an ionization fraction χ = 0.1, and includes a single population of charged dust silicate particles of radius 0.05 μm.

We show that charged dust grains follow the magnetic field and form filamentary structures with a typical width of 0.0035–0.0087 pc. Although the mean dust-to-gas ratio is globally preserved, local variations are observed that are not always correlated with gas density; these variations are in concordance with previous theoretical (Padoan et al. 2006; Hopkins 2014) and observational (Siebenmorgen et al. 2017; Chen et al. 2015; Reach et al. 2015) results.

Although a more realistic model (nonideal MHD effects and a full dust distribution) is needed, our results highlight the importance of considering the dynamics of charged dust for star formation studies even at a low-density regime characteristic of a molecular cloud envelope.

We want to thank an anonymous referee for many helpful comments and suggestions that have helped to improve this manuscript. L. B.-A. acknowledges Universidad Complutense de Madrid and Banco Santander for the grant "Personal Investigador en Formación CT17/17-CT18/17". This work has been partially funded by the Ministry of Economy and Competitiveness of Spain through grants MINECO-ESP2015-68908-R and MINECO-ESP2017-87813-R. L. B.-A. also wants to acknowledge I. Prada for a useful discussion about how to build the dust maps and with the usage of the python image analysis packages. This research has made use of NASA's Astrophysics Data System.

We designed two very simple tests to check that the movement of charged particles is correctly solved in 2D. Both tests consider a positively charged grain, initially at rest, embedded in a uniform medium moving at constant speed parallel to the x-axis. The constant magnetic field is parallel to the z-axis, perpendicular to the domain we are solving.

A.1. Lorentz Force Term

The first test is intended to check the integration accuracy of the Lorentz force term. The problem to be solved analytically is:

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{d{\boldsymbol{x}}}{{dt}} & = & {\boldsymbol{v}},\,{\boldsymbol{x}}(0)=({x}_{0},{x}_{1})\\ \displaystyle \frac{d{\boldsymbol{v}}}{{dt}} & = & -\displaystyle \frac{{\boldsymbol{v}}-{\boldsymbol{u}}}{{t}_{s}}+Q({\boldsymbol{v}}-{\boldsymbol{u}})\times {\boldsymbol{B}},\,{\boldsymbol{v}}(0)=(0,0).\end{array}\end{eqnarray*}$$

Depending on the ratio between t_s and Q∣ B ∣ ≐ QB, the particle motion will be dominated by gas drag (${t}_{s}^{-1}$ ≫ QB) or by magnetic effects (${t}_{s}^{-1}$ ≪ QB). Hence, we run four tests up to t = 10t_s for a fiducial stopping time t_s  = 1 and have considered the following values for the pair (Q, B): test A (0.1, 0.1), test B (0.1, 1), test C (0.1, 10), and test D (1, 10). For all the cases (see Figure 13) the error in position is of the order of 10⁻³ or lower, and the error in velocity is at most of the order of 10⁻², which confirms that our integrator for the Lorentz term works well and preserves the second-order accuracy. We want to note that for a magnetic-dominated medium (QB ∼ 100) it is necessary to diminish even more the timestep in Equation (6) because with the current implementation, the integration diverges after a few timesteps. However, for studies of dust in the ISM such harsh conditions are not expected, so the implementation of the integrator presented here is adequate for our purposes.

Zoom In Zoom Out Reset image size

A.2. Coulomb Drag

This test is defined to test the Coulomb drag term in the full Equation (1). We run test D (Q = 1, B = 10) for ν = 0.1, 1, 10 to study how it effects the transition from the aerodynamic regime to the magnetic one; the results are shown in Figure 14. The main effect of the electric drag is the damping of the magnetic oscillation of the charged grain, that can be rapidly suppressed for large values of ν; however, it will be typically lower than Q since ν ∝ Q² (Equation (2)). This implies that the overall accuracy of the integrator is kept second order, since the oscillations in the precision arise from the effects of the Lorentz force.

Zoom In Zoom Out Reset image size