© ESO, 2017

1. Introduction

Observations of star forming regions from telescopes such as Herschel (André et al. 2010, 2014; Arzoumanian et al. 2011; Motte et al. 2010), and IRAM (Hacar et al. 2013; Henshaw et al. 2014; Peretto et al. 2014; Tafalla & Hacar 2015; Beuther et al. 2015a,b), among others (e.g., Kirk et al. 2013; Li et al. 2015), have shown the ubiquity of filaments in star forming regions, however the conditions for the formation of these structures are not yet well known. Observations and theory (Fiege & Pudritz 2000; Padoan et al. 2001; Alves et al. 2008; Nakamura & Li 2008; Heitsch 2013; Arzoumanian et al. 2013; Hacar et al. 2013; Henshaw et al. 2014; Tomisaka 2014; Kirk et al. 2015; Seifried & Walch 2015; Moeckel & Burkert 2015; Li et al. 2015; Franco & Alves 2015, among others) indicate that turbulence and magnetic fields both have a role in the formation of these structures, however to what degree and in what form (e.g., driven or decaying turbulence) is still up for debate. In addition, the degree of ionisation within the molecular cloud needs to be taken into account. In the past couple of decades, observations of star forming regions have revealed a plethora of neutral and ionic molecules (e.g., CO, NH₃, N₂H⁺, etc.) that exist at and help probe various density thresholds. In a magnetised medium, the ions are bound to the field lines, however, depending on the radiation field and the ionisation fraction, neutral particles may be able to slip past the field lines to create density enhancements which may collapse into clumps and filaments. The ability for neutral particles to slip past the field lines is dependent on the degree of ambipolar diffusion within a region, however, the majority of simulations only consider a flux frozen medium (i.e., a medium where the neutral particles are coupled to the magnetic field lines due to frequent collisions with ions). Ambipolar diffusion of the neutral fluid relative to the ions is a necessary ingredient to understand both core formation and turbulence in partially ionised molecular clouds. For example, the preferred core fragmentation scale is a sensitive function of ionisation fraction and mass-to-flux ratio (Ciolek & Basu 2006; Bailey & Basu 2012). Bailey & Basu (2013) showed how this can lead to a much broader core mass function than pure hydrodynamic fragmentation. Furthermore, turbulent flows in a strongly magnetised medium can induce rapid ambipolar diffusion in compressed regions, leading to localised star formation with low overall efficiency (e.g., Nakamura & Li 2005; Basu et al. 2009a). Here we extend such models using planar models of molecular clouds threaded by a large scale magnetic field, including turbulence, ambipolar diffusion, and a step-like ionisation profile that reflects ionisation by background ultraviolet starlight in the low column density regions and cosmic ray induced ionisation in higher column density regions. These models have a background magnetic field that is perpendicular to the plane of the sheet, consistent with the presence of a dynamically important magnetic field, self-gravity, and turbulence that is Alfvénic or sub-Alfénic. These conditions allow settling of matter along the field lines and the formation of a flattened layer before further evolution occurs. Here we do not explore the opposite limit of an initially weak magnetic field that can be swept up in a sheet by turbulence or expanding supernova shells and be aligned parallel to a sheet. For example Nagai et al. (1998) performed a linear stability analysis of such a sheet and found that it can form elongated structures either perpendicular or parallel to the magnetic field direction in the cases of low or high ratio of external pressure to self-gravitational pressure respectively. Previous three dimensional simulations with dynamically important magnetic fields (Kudoh & Basu 2008, 2011) have shown similar results as the corresponding thin sheet results of Nakamura & Li (2005) and Basu et al. (2009a). The thin sheet models are actually more useful in studying initially near-critical mass-to-flux ratio clouds, since the critical mass-to-flux ratio is known analytically in a thin sheet geometry (Ciolek & Basu 2006) but not identifiable analytically for three-dimensional models with finite extent along the magnetic field direction. Some of the outstanding questions in this field are the values of the initial mass-to-flux ratio in molecular clouds and the level of initial turbulence and whether it is freely decaying or continues to be driven.

With the above issues in mind, our aim is to self consistently explore all of these effects. In a previous set of papers, we looked at the effects of the ionisation profile on the structure and kinematics of core collapse within non-ideal magnetic molecular clouds (Bailey & Basu 2014; Bailey et al. 2015, hereafter BB2014 and BBC2015, respectively). To further this analysis, we include the effects of turbulence in our simulations in order to see how/if the density, mass-to-flux ratio structures, and kinematics can change. We assess if the ionisation structure is playing a role in the dynamical evolution and structure of a cloud just before star formation. To do so, we explore the effects of changing the turbulent Mach number or initial mass-to-flux ratio in two separate sets of models (Set 1 and Set 2, respectively). In this paper we focus on the density and mass-to-flux ratio structures found within the two sets of simulations while a second paper will look at how the addition of turbulence affects the kinematics when different ionisation profiles are considered. The rest of this paper will be structured as follows. Section 2 describes the numerical code used to produce the simulations. Section 3 describes the model parameters investigated as well as the results of the simulations. Section 4 looks at the physical properties of filaments and cores found within each model and compares this analysis to other studies. Finally, in Sects. 5 and 6 we present our discussion and summary, respectively.

2. Numerical code

For our simulations, we use the Basu & Ciolek (2004) IDL code to model a weakly ionised, magnetic, self-gravitating, isothermal molecular cloud. Specifically, we are using the modified version (Bailey & Basu 2014) that includes density dependent ionisation profiles which regulate the degree of ambipolar diffusion within the molecular cloud. We explore the resulting kinematics, density, mass-to-flux ratio and velocity structures of simulated clump-core complexes within partially ionised, isothermal magnetic turbulent interstellar molecular clouds. The model assumes planar clouds with periodic boundary conditions in the x- and y-directions and a local vertical half thickness Z. This so-called “thin-sheet” approximation allows us to explore a significant parameter space while minimising computational resources and has been shown to provide results that are similar to three-dimensional non-ideal MHD simulations (Kudoh & Basu 2011; Nakamura & Li 2005; Basu et al. 2009a). Three-dimensional effects of e.g., turbulent support along the magnetic field direction (Kudoh & Basu 2003) and ionisation fraction variation along the vertical direction are not accounted for in this thin sheet model. A full description of the assumptions, non-axisymmetric equations and formulations and possible short-comings of this model can be found in Basu & Ciolek (2004), Ciolek & Basu (2006), Basu et al. (2009a,b), however we highlight those essential for the analysis within this paper.

2.1. Magnetic fields

The model assumes a magnetic field that threads the cloud perpendicular to the xy plane and includes the effects of ambipolar diffusion. The timescale for collisions between neutral particles and ions is ${\mathit{\tau }}_{\mathrm{ni}}\mathrm{=}\mathrm{1.4}\left(\frac{{\mathit{m}}_{\mathrm{i}}\mathrm{+}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{m}}_{\mathrm{i}}}\right)\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathrm{⟨}\mathit{\sigma w}{\mathrm{⟩}}_{\mathrm{i}{\mathrm{H}}_{\mathrm{2}}}}\mathit{,}$(1)where m_i is the ion mass, m_H2 is the mass of molecular hydrogen, n_i is the number density of ions, and ⟨ σw ⟩ _iH2 is the neutral-ion collision rate. Typical ions within a molecular cloud include singly ionised Na, Mg, and HCO which have an average mass of 25 amu. Assuming collisions between H₂ and HCO⁺, the neutral-ion collision rate is 1.69 × 10^-9 cm^-3 s^-1 (McDaniel & Mason 1973). The factor of 1.4 in Eq. (1) accounts for the fact that the inertia of helium is neglected in the calculation of the slowing-down time of neutrals by collisions with ions (Mouschovias & Ciolek 1999; Ciolek & Basu 2006).

Collapse within a molecular cloud is regulated by the normalised mass-to-flux ratio of the background reference state, ${\mathit{\mu }}_{\mathrm{0}}\mathrm{\equiv }\mathrm{2}\mathit{\pi }{\mathit{G}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{\sigma }}_{\mathrm{n}\mathit{,}\mathrm{0}}}{{\mathit{B}}_{\mathrm{ref}}}\mathit{,}$(2)where (2πG^1/2)^-1 is the critical mass-to-flux ratio for gravitational collapse in the adopted model (Ciolek & Basu 2006), σ_n,0 is the initial column density and B_ref is the magnetic field strength of the background reference state. Regions with mass-to-flux ratios greater than, approximately equal to, and less than one are known to be supercritical, transcritical and subcritical, respectively. In the limit where τ_ni → 0, the medium is defined to be flux-frozen, that is, frequent collisions between the neutral particles and ions couple the neutrals to the magnetic field. Under these conditions, subcritical regions are supported by the magnetic field and only supercritical regions may collapse within a finite time frame. Non-zero values of τ_ni are inversely dependent on the ion number density and therefore on the degree of ionisation for a fixed neutral density.

2.2. Turbulence

Turbulence within our simulations is incorporated via an initial “turbulent” velocity field added to the background reference state of uniform column density and vertical magnetic field strength B_ref. The velocity field is generated in Fourier space using the same method described in Basu et al. (2009a). For a physical grid containing N zones in both the x- and y-directions, there exists a corresponding wavenumber in Fourier space, k_x and k_y, respectively. For each (k_x, k_y) pair, a Fourier velocity amplitude is assigned from a Gaussian distribution and scaled by the power spectrum ${\mathit{v}}_{\mathit{k}}^{\mathrm{2}}\mathrm{\propto }{\mathit{k}}^{\mathit{n}}$, where $\mathit{k}\mathrm{=}\mathrm{(}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. The resulting Fourier velocity field is then transformed back into physical space. The distributions of v_x and v_y are chosen independently in this manner and each rescaled such that the rms amplitude is equal to the chosen turbulent velocity amplitude v_a. The turbulent velocity field is only added at the beginning of each simulation and allowed to decay. For our simulations, we test the effect of varying both the turbulent velocity amplitude and the initial mass-to-flux ratio within the context of varying ionisations profiles.

To enable comparison between models, we construct the Gaussian distributions for each model using the same realisation (i.e., the same random number generator seed). However, within this realisation, we ensure that the x and y velocity fields are generated using two distinct seed values. By maintaining the same realisation, we are able to distinguish which physical changes are due to changes in physical parameters as opposed to the random nature of the initial velocity distribution. The exact choice of random seed will affect the run time of the simulation, therefore quoted simulation times within this paper should be taken as representative but not exact values.

2.3. Ionisation profiles

The modified version of the IDL code includes density dependent ionisation profiles. For our analysis, we utilise two different ionisation profiles which describe the ionisation fraction (χ_i) as a function of density. These profiles set the neutral-ion collision time within each pixel since τ_ni ∝ 1 /χ_i. The first profile used is a step-like ionisation profile that includes both the ultraviolet (UV) and cosmic ray (CR) regimes. The form of this step-like ionisation function is given by $\log {\mathit{\chi }}_{\mathrm{i}}\mathrm{=}\{\begin{array}{c}\\ \log {\mathit{\chi }}_{\mathit{i,}\mathrm{0}}\mathrm{+}\mathrm{0.5}\mathrm{(}\log {\mathit{\chi }}_{\mathit{i,c}}\mathrm{-}\log {\mathit{\chi }}_{\mathit{i,}\mathrm{0}}\mathrm{)}& \\ \hspace{0.17em}\mathrm{\times }\left(\mathrm{1}\mathrm{+}\tanh \frac{{\mathit{A}}_{\mathrm{V}}\mathrm{-}{\mathit{A}}_{\mathit{V,}\mathrm{crit}}}{{\mathit{A}}_{\mathit{V,d}}}\right)& {\mathit{A}}_{\mathrm{V}}\mathrm{\le }{\mathit{A}}_{\mathrm{V}\mathit{,}\mathrm{CR}}\\ \log [\mathrm{1.148}\mathrm{\times }{\mathrm{10}}^{-7}\mathrm{(}\mathrm{1}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}& \\ \mathrm{\times }{\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}{\frac{\mathrm{2.75}\mathrm{mag}}{{\mathit{A}}_{\mathrm{V}}}}^{\mathrm{)}}]& {\mathit{A}}_{\mathrm{V}}\mathit{>}{\mathit{A}}_{\mathrm{V}\mathit{,}\mathrm{CR}}\mathit{,}\end{array}$(3)where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}$ is the dimensionless external pressure and A_V,CR = 6.365 mag is the location of the transition from the UV regime to the CR regime. The step function parameters are set to log χ_i,0 = −4.0, log χ_i,c = −7.362, A_V,crit = 4.0 mag, and A_V,d = 1.05 mag. This ionisation profile is based upon those presented by Ruffle et al. (1998), where we have chosen our high and low ionisation levels (log χ_i,0 and log χ_i,c, respectively) to depict a rough average of the different profiles presented. The location of the step was chosen to correspond with the typical A_V where dust shielding sets in, thus eliminating the effects of the UV radiation. Naturally, other step-like profiles can be constructed in this way, however we have chosen to focus on this one to facilitate comparison with our previous analysis and simulations (see Bailey & Basu 2012, 2014; Bailey et al. 2015). With the inclusion of this ionisation profile, the ionisation fraction regime (and neutral-ion collision time) within the cloud depend on the column density/visual extinction of each region and therefore changes as the cloud evolves. The derived neutral-ion collision time for each pixel then depends on which ionisation regime it is in. For the UV regime, the neutral-ion collision time is computed via $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{ni}}\mathit{/}{\mathit{t}}_{\mathrm{0}}\mathrm{=}\mathrm{0.262}{\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\left(\frac{\mathrm{0.01}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}{{\mathit{\sigma }}_{\mathrm{n}\mathit{,}\mathrm{0}}}\right)\left(\frac{{\mathrm{10}}^{-7}}{{\mathit{\chi }}_{\mathrm{i}}}\right)\mathrm{(}\mathrm{1}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}{\mathrm{)}}^{-1}\mathit{,}& & \end{array}$(4)where t₀ is the dimensionless time scale (t₀ = c_s/ 2πGσ_n,0). For the CR regime, we set the dimensionless neutral-ion collision time to τ_ni,0/t₀ = 0.2.

The second ionisation profile corresponds to a cosmic ray (CR) only profile which is achieved by assuming a dimensionless neutral-ion collision time of τ_ni,0/t₀ = 0.2 for all densities.

2.4. Dimensionless and dimensional values

Our model is characterised by several dimensionless free parameters including a dimensionless form of the initial neutral-ion collision time (τ_ni,0/t₀ ≡ 2πGσ_n,0τ_ni,0/c_s) and a dimensionless external pressure ($\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}\mathrm{\equiv }\mathrm{2}{\mathit{P}}_{\mathrm{ext}}\mathit{/}\mathit{\pi G}{\mathit{\sigma }}_{\mathrm{n}\mathit{,}\mathrm{0}}^{\mathrm{2}}$). Here, c_s = (k_BT/m_n)^1/2 is the isothermal sound speed; k_B is the Boltzmann constant, T is the temperature in Kelvin, and m_n is the mean mass of a neutral particle (m_n = 2.33 amu). We normalise column densities by σ_n,0, length scales by ${\mathit{L}}_{\mathrm{0}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{2}\mathit{\pi G}{\mathit{\sigma }}_{\mathrm{n}\mathit{,}\mathrm{0}}$ and time scales by t₀ = c_s/ 2πGσ_n,0. Based on these parameters, typical values of the units used and other derived quantities are $\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathrm{n}\mathit{,}\mathrm{0}}\mathrm{=}\frac{\mathrm{3.63}\mathrm{\times }{\mathrm{10}}^{-3}}{\mathrm{(}\mathrm{1}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\left(\frac{{\mathit{n}}_{\mathrm{n}\mathit{,}\mathrm{0}}}{{\mathrm{10}}^{\mathrm{3}}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}{\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}\mathit{,}\\ & & {\mathit{c}}_{\mathrm{s}}\mathrm{=}\mathrm{0.188}{\left(\frac{\mathit{T}}{\mathrm{10}\mathrm{K}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{km}{\mathrm{s}}^{-1}\mathit{,}\\ & & \\ & & \end{array}$where n_n,0 is the initial neutral number density. For our analysis, we assume a dimensionless external pressure $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{P}}_{\mathrm{ext}}\end{array}\mathrm{=}\mathrm{0.1}$ (P_ext ≈ 10³ cm^-3 K), a neutral number density of 1.1 × 10³ cm^-3, and a temperature T = 10 K.

3. Simulations

3.1. Model parameters

The models and analyses presented in this paper form an extension to the analyses presented in BB2014 and BBC2015 to include turbulence. We assume an initially diffuse cloud with an initial background column density which corresponds to a visual extinction A_V,0 = 1 mag. From the prescription of Pineda et al. (2010) (see also Bailey & Basu 2012) and assuming a mean molecular weight of 2.33 amu, the resulting conversion between visual extinction and mass column density is ${\mathit{\sigma }}_{\mathrm{n}}\mathrm{=}\mathrm{3.638}\mathrm{\times }{\mathrm{10}}^{-3}\mathrm{(}{\mathit{A}}_{\mathrm{V}}\mathit{/}\mathrm{mag}\mathrm{)}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}\mathit{.}$(9)All simulations begin with an initial turbulent velocity perturbation assuming a power spectrum exponent n = −4, such that ${\mathit{v}}_{\mathit{k}}^{\mathrm{2}}\mathrm{\propto }{\mathit{k}}^{-4}$ and are performed on a 512 × 512 periodic box. The box size is 64πL₀, which translates to a size of 15.16 pc, or a pixel size of 0.0296 pc, for T = 10 K and σ_n,0 = 3.638 × 10^-3 g cm^-2. This results in a total mass of 4000 M_⊙ within the simulation region. All simulations terminate when any pixel within the region attains a density enhancement that is 10 times the initial column density (σ_n,0). This threshold was chosen as it corresponds to the runaway gravitational collapse of a cloud into a dense core. Basu et al. (2009a) verified, with high resolution runs, that collapse does indeed continue past this value. In addition, we chose this stopping conditions because we are interested in how the physical conditions of the environment affect the ability of a molecular cloud to achieve runaway collapse and not the subsequent evolution of the high density regions formed. For our simulations, we are looking at both the effects of changing the turbulent Mach number and the degree of magnetic support in the context of different ionisation profiles. To achieve this, we present the results of two sets of simulations. The first set (Models I−X) assumes the same initial mass-to-flux value while varying the initial Mach number. For the second set (Models XI−XVI), we have fixed the initial Mach number and vary the initial mass-to-flux ratio. The initial parameters for both model sets are listed in Table 1. Here, v_a/c_s is the turbulent velocity amplitude scaled by the sound speed (see Sect. 2.2), which is equivalent to the Mach number, t_r refers to the total run time of the simulation, and σ_n,max, μ_max and μ(σ_n,max) refer to the maximum density, maximum mass-to-flux ratio and mass-to-flux ratio at the location of the maximum density within the simulation at the final time, respectively. The abbreviations quoted for the ionisation profile indicate whether the simulations assumes a step-like (SL) or cosmic-ray-only (CR) ionisation profile. These two profiles correspond to an initial neutral-ion collision time τ_ni,0/t₀ = 0.001 and τ_ni,0/t₀ = 0.2, respectively. The characteristic free-fall time for our system is 1.01 Myr. From the linear analysis of Bailey & Basu (2012), the ambipolar diffusion time is equivalent to the fragmentation time scale for μ₀ = 0.5 and is dependent on the initial ionisation fraction of the medium. This corresponds to values of 1.56 Gyr and 7.88 Myr for the SL and CR-only models, respectively. For the initial mass-to-flux ratio value used in Set 1 (μ₀ = 1.1), the characteristic times found through linear analysis are 7.2 Myr and 4.1 Myr for the SL and CR-only models, respectively.

The simulations track the evolution of various physical properties of the cloud including the density enhancement, velocity, magnetic field strength, mass-to-flux ratio and the ionisation. Given the wealth of information available from these simulations, we are splitting the discussion of these parameters into two separate papers. In this paper, we will focus on the structures formed within the simulations by looking at the density enhancements and mass-to-flux ratio maps. The second paper will focus on the kinematic information revealed through the velocity maps and synthetic spectra of the simulations. The following subsections present the effect of varying the set-specific initial conditions on the density enhancement and mass-to-flux ratio maps, paying particular attention to the difference between the two different ionisation profiles.

3.2. Effect of varying the initial mach number

3.2.1. Density enhancement structures

In this section, we look at the effect of changing the turbulent Mach number on the density enhancement structures formed within a collapsing molecular cloud (Models I–X). We are specifically interested in the differences due to the underlying ionisation profile. Figure 1 shows a sample of the density enhancement maps (with respect to the initial density) for three values of the turbulent Mach number (0.03 (left), 1.0 (centre), and 2.0 (right)) for both the step-like (top row) and CR-only (bottom row) ionisation profiles. Models VII–X do not show significant differences compared to models V and VI and therefore are not shown, however, these models do exhibit shorter collapse timescales than the regions with lower turbulence (see Table 1). The colour scheme depicts the logarithm of the density enhancement. Looking across the three columns, we see that the assumed ionisation profile results in the largest difference in the simulations with the least amount of turbulence (Models I and II) and the smallest, virtually indistinguishable difference in the simulations with the strongest turbulence (Models V and VI). This indistinguishability is also observed in models for Mach 3 and 4 (not shown). Focussing on the SL models (top row), we see that the degree of fragmentation increases from a single monolithic filament/clump in Model I (v_a/c_s = 0.03) to multiple thin filaments in Model V (v_a/c_s = 2.0). Conversely, the CR-only models (bottom row) show a decrease in the degree of fragmentation from single isolated cores in Model II to thin filaments in Model VI. Finally, we note that our “non-turbulent” simulations (Models I and II) are very similar to the corresponding micro-turbulent simulations in BB2014 and BBC2015.

The formation of either a filament or a core is a direct consequence of the state of the medium. As shown in Bailey & Basu (2012, 2014), the fragmentation length scale of a cloud depends on the mass-to-flux ratio and the ionisation profile. The transcritical nature of these simulations puts them in the regime where a change in the ionisation fraction can drastically change the fragmentation length scale. As shown in Fig. 2 (right panel) of Bailey & Basu (2012), for a nearly flux frozen medium (τ_n,i = 0.001, i.e., models with SL ionisation profile) the fragmentation length scale is much larger than in a medium with a smaller ionisation fraction (τ_n,i = 0.2, i.e., models with CR-only ionisation profile). The thinning out of the structures found in Models I and II to those observed in Models V and VI is due to the turbulence present in the medium.

3.2.2. Mass-to-flux ratio structure

In this section, we focus on the effect that varying the initial turbulent Mach number and ionisation profile have on the distribution of mass-to-flux ratio within the first set of simulations. The effect of varying the initial mass-to-flux ratio will be explored in the following section. Figure 2 shows representative mass-to-flux ratio maps for CR-only models with Mach number equal to 0.03, 2 and 4 (Models II, VI and X, respectively). From these maps, we see that the structure of the mass-to-flux ratio within the regions changes on smaller and smaller scales as the degree of turbulence increases. This is also evident in the corresponding SL simulations (not shown), however due to the nearly flux frozen nature of these simulations, the range of mass-to-flux ratios present within the region are not significantly different from the initial value (μ₀ = 1.1) and therefore do not show up within a linearly scaled map. As indicated in Table 1, there is a significant difference between the maximum mass-to-flux ratio values found in the SL and CR-only models; the average maximum value for SL models is 1.186 while the average maximum value for the CR-only models is 1.352. Looking at the low density regions, we see that the structure within the CR-only models exhibits ~0.1 differences between the high and low mass-to-flux ratio regions (green vs. blue) while in the SL simulations, the differences are on the order of 0.003 (not shown). The CR-only models are able to achieve a higher peak mass-to-flux ratio due to the lower ionisation fraction within these simulations as compared to the SL models. This allows for more neutrals to slip past the magnetic field lines to increase the density while locally keeping the magnetic field strength the same and thus increasing the mass-to-flux ratio.

The top panel of Fig. 3 shows a representative map of the mass-to-flux ratio structure within the dense region of a SL model (Model V). Looking at the clump/core regions in detail we see that the models exhibit mass-to-flux ratio structures similar to those described in Bailey & Basu (2014) (Model A), i.e., a high mass-to-flux ratio corresponding to the highest density region surrounded by an extremely low mass-to-flux ratio on either side. Specifically, for the core depicted in Fig. 3, we see that the high density central region has a mass-to-flux ratio of ~1.18, while the low density outer regions have a mass-to-flux ratio of only ~1.04. The bottom panel of Fig. 3 shows profiles of the mass-to-flux ratio across the core as a means to further illustrate this interesting structure. The black, solid line shows the profile for the core in the top panel, while the red, dashed line shows the profile for the corresponding core in Model VI. As shown by these two profiles, this characteristic structure/profile is evident within all the models, i.e., a spike at the location of the core and depressions on the outskirts as described above. This is a clear signature of flux redistribution by ambipolar diffusion. The profile is more pronounced in the SL models however due to the smaller variations in the background mass-to-flux ratio. In contrast, the feature can be much broader within the CR-only simulations. In this case, the flux redistribution occurs over a larger region than in the SL models. Initial analysis of the kinematics show that Model V (SL) and Model VI (CR-only) have similar velocity structures near the core, while Model III (SL) has a steeper velocity gradient near the core than in Model IV (CR-only), thus suggesting that the SL and CR simulations present distinguishable velocity structures for apparently similar density enhancement structures. Further analysis of the kinematics will be presented in our next paper.

Figure 4 shows the maximum mass-to-flux ratio within the simulation region as a function of time for the ten models in Set 1 (μ₀ = 1.1). Solid lines indicate simulations with SL ionisation profiles while dashed lines indicate CR-only ionisation profiles. Colours indicate different Mach values (Black = 0.03, Red = 1, Green = 2, Blue = 3, and Purple = 4). The inset shows a zoom-in of the first 10t₀. Note that the location of the maximum mass-to-flux ratio (μ_max) does not necessarily correspond to the location of the maximum density enhancement within the simulation. This is evident by the fact that the μ_max values quoted in Table 1 do not necessarily equal the mass-to-flux ratio at the location of the maximum density enhancement μ(σ_n,max). This mismatch between mass-to-flux ratio values can occur throughout the evolution of the region.

Looking at the set as a whole, it is clear that the increase in the Mach number results in a marked decrease in the evolutionary time. In addition we note that for low value Mach numbers (i.e., those in Models I–IV), the SL simulation takes up to two times longer to evolve compared to the corresponding CR simulation, however for Mach numbers greater or equal to 2, the evolutionary times for the two different profiles are more or less identical. These temporal comparisons can also be seen in Table 1. Focussing on the evolution of the mass-to-flux ratio itself, we see that all SL profile models exhibit a constant mass-to-flux ratio for the majority of the run with a sharp increase at near the end. This mimics the evolution of the maximum density and is a result of the nearly flux frozen nature of the medium at lower densities. In contrast, the CR models show a steady, mostly linear increase over the entire duration of the simulation. Again this is a direct consequence of the increased ambipolar diffusion throughout the simulation region in these models. Finally, we can see that the CR simulations result in greater final maximum mass-to-flux ratios than the SL models. These different mass-to-flux ratios should imply different dynamical histories. Therefore, as mentioned before, we expect the velocity structure to be different between the two ionisation profiles which can cause variations in simulated observations of molecular cloud tracers (Bailey et al., in prep.) and enable us to quantify the effects of the ionisation structure and magnetic fields via molecular line observations.

3.3. Effect of mass-to-flux ratio on turbulent simulations

In addition to the models discussed in the previous sections (Set 1: Models I–X) we also ran a second set of models (Set 2: Models XI–XVI) which look at the evolution of the medium for different initial mass-to-flux ratios while keeping the initial turbulent Mach number fixed (v_a/c_s = 2.0). For this set, we have explored both sub- and supercritical initial mass-to-flux ratios which we can compare to the transcritical runs from Set 1 (Models V and VI for SL and CR, respectively). The initial parameters and resulting values are noted in Table 1. In the following subsections, we look at the effect of changing the initial mass-to-flux ratio on the density enhancement structures and mass-to-flux ratio structures. Note that Models XI and XIII did not collapse into high density structures, therefore the values in Table 1 reflect those at the time when these models were terminated manually.

3.3.1. Density enhancement structures

Linear analysis from Bailey & Basu (2012) shows that although fragmentation of a subcritical region may be possible, it can take up to infinitely longer than media with trans- or supercritical initial mass-to-flux ratios, however this analysis does not take into account the possible effects of turbulence on the fragmentation of a subcritical region. Looking at the results of the SL-simulations, we see that the addition of a Mach 2 turbulent component does not add enough energy to create regions where ambipolar diffusion can overcome the magnetic field and allow a subcritical medium (Models XI and XIII) to evolve to a point where runaway gravitational collapse can occur. For Model XIII, an initial run which assumed our fiducial stopping condition seemed to suggest that runaway collapse could occur even though the medium is sub-critical. However, a subsequent run with an increased stopping condition (σ/σ_n,0 = 20) revealed that although a maximum density enhancement of 12.82 was achieved at t/t₀ = 11, runaway collapse was not achieved. Instead, these dense regions were still subcritical and transient in nature. Figure 5 shows density enhancement maps at the final time for each simulation in Set 2, as quoted in Table 1. The top row shows the SL models while the bottom row shows the CR-only models. The initial mass-to-flux ratio values increase from left to right (μ = 0.5,0.8 and 2.0, respectively). The colour scale shows the logarithm of the density enhancement. Looking at the subcritical SL models which did not result in runaway gravitational collapse (Model XI and Model XIII, top left and middle respectively), we see that the medium has become quite smooth with some evidence of low density filaments. These filaments tend to be transient and do not persist long enough for cores to form within them. For the supercritical SL case (Model XV), the medium fragments into dense filaments as expected. In this case, the higher mass-to-flux ratio allows this fragmentation to occur on a shorter timeframe than in the corresponding model from Set 1 (Model V). Looking at the bottom row, we see that with the CR-only ionisation profile, all three models are able to fragment into dense cores/filaments. Interestingly, change in the resulting structures via the increase in mass-to-flux ratio mimics the evolution seen with increasing turbulence, i.e. the smallest mass-to-flux ratio yields small isolated compact cores while larger mass-to-flux ratios result in a stretching of these single cores into filamentary structures with evidence of multiple cores within (see also Basu et al. 2009b,a).

3.3.2. Mass-to-flux ratio structures

Figure 6 shows mass-to-flux ratio maps at the final time for the Set 2 CR-only simulations as quoted in Table 1. As with the SL models from Set 1, the variations in the mass-to-flux ratios for SL models in this set are not significantly different from the initial value and therefore do not show up well on a linear scale. The initial mass-to-flux ratio value increases from left to right. The colour scale shows the mass-to-flux ratio on a linear scale. Note that the minimum and maximum scale for each panel is not the same. The maximum mass-to-flux ratio values present at the displayed times can be found in Table 1. Looking at differences between the SL and CR simulations, we see that the CR-only simulations result in larger changes from the initial mass-to-flux ratio values than the SL only simulations. This is due to the SL simulations being held at near flux-frozen conditions for the majority of the run. Note that the large difference between the final mass-to-flux ratio values for the μ = 0.5 (Models XI and XII) and μ = 0.8 (Models XIII and XIV) is due to the fact that the SL simulations were not able to collapse and therefore were always held at near flux frozen conditions. Looking at the bottom row of Fig. 6, the CR-only simulations show that as the mass-to-flux ratio increases, the amount of structure within the mass-to-flux ratio map increases. However, the relative amplitude of variations is greater for the lower mass-to-flux ratio models. Looking back at Model VI (μ₀ = 1.1), we see that the degree of structure present in that simulation fits in with the pattern in the three CR-only simulations in Set 2. Analysis of the profiles across the cores within the models of this set shows that, as with the other set, both the SL and CR-only simulations show the same characteristic shapes with the SL models again exhibiting spikier profiles due to the smaller variation in mass-to-flux ratio between the low and high density regions. We note that the fluctuations in mass-to-flux ratio are introduced in our initial conditions, and we are tracking the difference in their subsequent evolution.

4. Physical properties of filaments and cores

The structures visible within a star forming region are highly dependent on the sensitivity of the instrument used to observe them. Data from the Herschel space telescope, the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope, the James Clerk Maxwell Telescope (JCMT) and, now even more from the Atacama Large Millimetre/submillimetre Array (ALMA) and the Jansky Very Large Array (JVLA), show that complex structures are typically observed in the direction of star forming regions. These structures can take on many different forms including subfragmentation of large clump regions into several cores (i.e., Sadavoy et al. 2012), ubiquitous filaments in star forming regions (André et al. 2010; Hacar et al. 2013; Henshaw et al. 2014), bundles of fibres within filaments (Hacar et al. 2013; Tafalla & Hacar 2015) and even filaments within dense cores (Pineda et al. 2011). Our simulations show formation of both filaments and dense cores.

As shown in Figs. 1 and 5, the majority of the models exhibit the formation of dense structures in four main regions. Figure 7 shows a representative map of these four visually recognisable regions which contain high density filaments within our simulations¹. This map is constructed using Model V, however please note that this map has been shifted by 50 pixels (20 map units) in the x direction in order to ensure that the filaments contained within region 1 are not split at the periodic boundary of the box. To avoid confusion however, all quoted locations of regions will be based upon the original unshifted figures (i.e., those shown in Fig. 1). The following sections will examine the physical properties of the filaments and embedded cores within the four regions defined above. Specifically we will determine the masses of the cores and filaments and analyse the profiles across the filaments to determine the widths via two different fitting functions.

4.1. Filament and core masses

To determine the mass within the filaments and dense structures in the simulations we utilise CLUMPFIND2D (Williams et al. 1994). This is a set of IDL routines which find the extent of structures within observations or simulated data. Structures are determined by assuming linearly spaced contours based upon user definitions and connecting pixels at each contour level that are within one resolution element of each other (Williams et al. 1994). For our simulations we have defined the outer edges of the filaments to be at A_V = 3 mag and dense structures to be regions above A_V = 7 mag. Recall that for our simulations, A_V refers to the visual extinction perpendicular to the 2D plane. Output from these routines include the “intensity” of the pixels within the identified structure. In the case of our simulations, this intensity is the sum of the column density over all pixels within the defined structure. From these data we then determine the enclosed mass of the filaments and dense structures. Finding filaments and cores within observational data usually requires more sophisticated techniques such as morphological component analysis (MCA) or use of programs such as getfilaments (Men’shchikov 2013) or DisPerSe (Sousbie 2013). Comparison of these techniques by Könyves et al. (2015) showed they all give approximately the same output. In our case, however, we do not require such sophisticated methods as our models are not affected by the noise and background fluctuations that are present in astronomical observations.

Table 2 shows the results of this analysis for the filaments within each of the models. For some models, the regions show evidence of multiple filamentary structures. This occurs when the density falls below our defined lower threshold such that the lowest defined contours enclose more than one filamentary structure. Note that we have not performed this analysis on Models II, XI, XII, XIII and XIV since they either do not exhibit filamentary structure within the region² or the maximum density within the simulation at the final time lies below our density definition of a filament. As shown in Table 2, the mass contained within a filament tends to decrease dramatically as the degree of turbulence increases. This continues until a Mach value of 2. For higher turbulence values, the mass contained seems to remain relatively constant. Comparing the two ionisation profiles, we see that for simulations with Mach <2, the filaments within the SL models tend to have larger masses than the CR-only simulations. Above Mach 2, the masses contained within the filaments tend to be equal across the two different ionisation profiles. Comparing the two model sets (Models V and VI vs. Models XV and XVI) we see that the regions with higher mass-to-flux ratios have smaller filament masses. This is due to both the shorter collapse time and smaller collapse length scales associated with regions containing supercritical mass-to-flux ratios.

Table 3 shows the results of this analysis for the cores (i.e., coherent regions with A_V> 7 mag). Looking at the locations of the resulting cores, we see that the majority of the cores tend to form in regions 1 and 4. This indicates that although the other two regions have formed filaments, the density is not high enough to have started to form cores yet. Note that the cores in the Set 2 models do not coincide with any of the defined filament regions. From the resulting masses, we again see that, as the degree of turbulence increases, the mass contained within the cores decreases. Note that these masses are for the full core regions found within each simulation at the end of the run, however this does not mean that each core will collapse into the same mass star. Each core still requires further evolution before a star can be formed which could include further fragmentation into several smaller cores. Comparing the two ionisation profiles, we see that on average, the cores evolving within the SL models accumulate more mass than those in the CR-only models. This is in contrast to the findings of Chen & Ostriker (2014) who found that their core masses did not depend on the degree of ambipolar diffusion. Such a discrepancy could be due to the differences in model assumptions.

Comparing the masses from Set 2 to Models V and VI, we see that in all cases, the models with sub- and supercritical mass-to-flux ratios exhibit smaller filament and core masses. The smaller masses found within the Set 2 models is entirely due to the sub- or supercritical initial mass-to-flux ratio. As shown in the linear analysis of Bailey & Basu (2012), structures formed within regions with magnetic fluxes that are significantly larger or smaller than the mass result in smaller fragmentation length scales than in those regions where the mass and magnetic flux are comparable.

4.2. Filament width as a function of Mach number

The discovery of filaments pretty much everywhere within star forming regions has promoted significant efforts in determining their widths and overall structure within different environments. Initial studies of the widths of filaments discovered through sub-millimetre dust continuum observations of the Gould Belt with the Herschel space telescope tend to show a characteristic width of 0.1 pc (Arzoumanian et al. 2011). Studies of other regions such as the Taurus molecular cloud (TMC; Malinen et al. 2012) and B211 (Palmeirim et al. 2013) also revealed this seemingly canonical value. Studies of Galactic cold core regions as chosen from the Cold Clump Catalogue of Planck Objects (C3PO) by Juvela et al. (2012) found that some regions (e.g., G94.15+6.50 and G298.31-13.05) conformed to a characteristic width of 0.1 pc, however other regions revealed much larger widths of up to 1 pc. Other studies also find regions with filaments that do not exhibit a width of 0.1 pc. A search for filamentary structure in ¹³CO data of Taurus reveal a broad distribution of profile widths peaking at 0.4 pc (Panopoulou et al. 2014), observations of filaments within the Cygnus X region find a mean profile width between 0.26 and 0.34 pc (Hennemann et al. 2012) and Herschel observations of the Galactic plane via the Hi-GAL project find widths ranging from 0.1 to 2.5 pc (Schisano et al. 2014). Numerical studies have also found a wide range of filament widths. For example, Smith et al. (2014) finds a range of FWHM between 0.2 pc and 0.35 pc while Kirk et al. (2015) find filament FWHMs that range from 0.05 pc to 0.12 pc. Ntormousi et al. (2016) found a difference in the peak of the filament width distribution depending on whether the region assumed ideal or non-ideal MHD, however, without the inclusion of self-gravity in their models, they were unable to conclude whether this was due to the different filament selection or the physics inherent to filament formation. Theoretical considerations by Hocuk et al. (2016) of the Jeans length within a non-magnetic medium, where the gas and dust temperature is followed in detail, suggests that the length scale is on the order of 0.1 pc, while studies of both non-self-gravitating and self-gravitating filaments tends to suggest that dissipation via ion-neutral friction may be a promising mechanism to explain the structure of self-gravitating filaments (Hennebelle 2013; Hennebelle & André 2013). In addition to this, recent studies of the L1495/B213 regions in Taurus tend to suggest that large scale filaments may be composed of thread-like bundles (Hacar et al. 2013). Given the discrepancy between all these different studies, it is only natural for us to analyse our own filaments to determine how the width varies with the Mach number and if our simulations agree/disagree with previous studies.

To determine the width of the filamentary structures within our simulations we utilise the fitting techniques of previous studies (e.g., Arzoumanian et al. 2011; Smith et al. 2014; Kirk et al. 2015). Within these techniques, one constructs an average profile across the filament and then fits this profile with a proper function from which a width can be extracted. The average profile is constructed by taking several density profile cuts along the length of the filament. The maximum of each profile is then shifted to coincide with each other and the average profile computed. For our simulations, we will examine the average profiles for the filaments found via CLUMPFIND2D listed in Table 2. Our profiles, are constructed by taking cuts across the filament for each pixel along the length of the filament. We assume that the maximum of each profile within the average occurs at r_c = 0, (where r_c is the radius from the centre of the filament) and the profile is followed out to 1 pc from this maximum on both sides. We chose this range to both ensure that the entire width of the filament is represented within the average profile and to enable comparisons with previous studies that also adopt this profile extent.

As described by Arzoumanian et al. (2011), Smith et al. (2014) and Kirk et al. (2015), among others, the average profiles are fit by either a Plummer profile or a Gaussian profile. For our analysis, we will also fit these two profiles. We utilise the 2D surface density Plummer profile from Smith et al. (2014), $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{A}}_{\mathit{p}}{\mathit{n}}_{\mathrm{c}}{\mathit{R}}_{\mathrm{flat}}}{\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{flat}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{]}}^{\mathrm{(}\mathit{p}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}}\mathrm{+}\mathit{B}\mathrm{\left[}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{\right]}\mathit{,}$(10)where n_c is the central density, R_flat is the radius within which the profile is flat, p is the slope of the power law decrease outside this radius and B is the background density. Note that p = 2 corresponds to Bonnor-Ebert model clouds (e.g. Dapp & Basu 2009). The A_p value in the relation above is a finite constant normalisation factor based on an assumed distribution of filament inclinations with respect to the observer as given by ${\mathit{A}}_{\mathit{p}}\mathrm{=}\frac{\mathrm{1}}{\cos \mathit{i}}{\mathrm{\int }}_{\mathrm{\infty }}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{u}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{u}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{p}\mathit{/}\mathrm{2}}}\mathit{,}$(11)where i is the inclination of the filament on the plane of the sky, generally assumed to be 0 (Arzoumanian et al. 2011). Arzoumanian et al. (2011) show that non-zero inclinations can affect the derived density by up to 60% but have no effect on the shape of the profile, however given the two dimensional nature of our simulations, it is natural to assume an inclination angle of zero. For the cases where p = 2 or 4, with i = 0, A_p = π/ 2, which is the same value used by Smith et al. (2014). For our analysis, we will also adopt this value for A_p but still allow the value of p to vary. Errors in the average profile are assumed to be equal to the maximum and minimum density values of the individual profiles at each distance from the maximum density of the filament. When these errors are taken into account in computing the profiles, the wings of the profile can sometimes prevent finding a good fit. To combat this problem, in these situations we truncate the average profile to only include the main density peak centred around r_c. This ensures that we are only fitting the primary filament within each region and that the widths are not contaminated by other nearby filaments.

Our Gaussian function takes the form $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{A}}_{\mathrm{g}}\exp \left(\frac{\mathrm{-}\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\mathit{,}$(12)where A_g is the Gaussian amplitude, r₀ is the mean of the Gaussian (which by definition should be zero for all of the profiles) and σ is the dispersion of the Gaussian.

4.2.1. Results

For this analysis, we attempt to fit the above two profiles to all of the filaments identified in Table 2. Figure 8 shows representative examples of three of the fits. The left and middle panels show the examples of fits for the lowest and highest Mach values, Model I (Mach 0.03) and Model IX (Mach 4) respectively. The right hand panel shows an example of a profile with a complicated wing structure (Model VII, Mach 3, Filament 3), for which we show an example of a truncated Plummer profile. In all panels, the red curve shows the average profile, the grey bars indicate the error range based upon the maximum and minimum density values at that particular radius from the centre of the filament. The blue curve depicts the Gaussian fit while the black curve depicts the Plummer profile. The results of all the fits can be found in Table 4. Note that filaments IV-1b and V-1b are missing from the table because it was not possible to find a satisfactory fit for either the Gaussian function or Plummer profile. Figure 9 shows the results from Table 4 in a visual form. Closed symbols depict the SL models while open symbols depict CR-only models. The red squares denote the filament which contains the region with the maximum density within the simulation. X-axis labels indicate the model number (roman numerals) and associated Mach number. The top left plot shows the FWHM of the Gaussian fit as a function of model. The other three panels show the Plummer profile parameters: p (bottom left), R_flat (top right) and n_c (bottom right). Focussing on the upper left panel first, we find that there is a scatter in the FWHM between the filaments within each model, however the overall FWHM tends to decrease as a function of Mach number. Above Mach 2, the scatter within a single model decreases significantly. Looking at the filaments which contain the highest density region (red squares) we see that these filaments tend to have the lowest FWHM amongst all the filaments within a particular model for Mach values less than 2. Above Mach 2, these regions no longer exhibit the smallest FWHM.

Looking at the three panels for the Plummer profile parameters, we see that there does not seem to be a discernible trend in the fitted p value, however we see that for increasing Mach number, the flat region of the Plummer profile (R_flat) decreases in size while the central density of the filament tends to increase. Looking at the right panels, we see that the filaments which contain the densest regions within the simulations tend to have the smallest central flat Plummer profile regions (upper right panel) and for the most part, the largest central densities (lower right panel). Finally, comparing the two sets of models (Models V and VI vs. Models XV and XVI) we see that the larger mass-to-flux ratio results in overall smaller FWHM and R_flat values. Looking at the p values, we see that their range is still scattered, however the filaments that contain the densest region of the simulation have the smallest p values. Finally, the models with the higher mass-to-flux ratios have larger central filament densities, particularly for the densest regions in the simulations (red squares).

4.2.2. Comparison to previous studies

Computational

In the past couple of years since the discovery of the seemingly characteristic width of 0.1 pc for filaments within star forming regions (Arzoumanian et al. 2011), several studies investigating the origin of this value (e.g., Smith et al. 2014; Kirk et al. 2015; Federrath 2016, among others) have been published. Here we will compare the results of our filament analysis to these three studies.

Before comparing, however, we first note the similarities and differences between all six simulations. Table 5 shows the important defining features of each of the six set ups in an easy to compare manner. For papers which include various physical set ups, (e.g., driven and decaying turbulence, HD and MHD etc.) we will only compare, if possible, to the simulations and results that closest match our set up (i.e., simulations with decaying turbulence and/or magnetic fields). As shown in Table 5, the most notable differences between the four different studies are the absence of magnetic fields in Smith et al. (2014), the absence of ambipolar diffusion in all three (Smith et al. 2014; Kirk et al. 2015; Federrath 2016), the absence of a treatment for chemistry in Federrath (2016) and the two-dimensional nature of our simulations compared to the others. In addition, note that both our simulations and those of Smith et al. (2014) assume an initial cloud that is 7.5 times larger than the clouds in the other two.

The two parameters related to the filament width are the FWHM and R_flat for the Gaussian and Plummer profiles, respectively. As defined in Arzoumanian et al. (2011) the FWHM from the Gaussian profile approximately corresponds to 3 times R_flat for p = 2. By fitting the Gaussian function to the profiles, Smith et al. (2014) found average FWHM values ranging from ~0.2 pc to 0.35 pc depending on the radius from the centre of the filament considered (0.35 pc vs. 1.0 pc) while Kirk et al. (2015) find FWHM values between 0.05 pc and 0.12 pc. For our simulations, we find that the FWHM decreases as function of Mach number with a range between 0.11 pc and 1 pc. Federrath (2016) does not quote FWHM values, however, they do find that they are consistent with the 3R_flat definition. Looking at this criteria for our simulations, as well as those by Smith et al. (2014) and Kirk et al. (2015), we see that the majority of the models do not satisfy this relation. This is due to the fact that in these analyses, the value of p was not fixed to 2 and thus the relation does not hold.

For the extent of the flattened inner region of the filament (R_flat), Smith et al. (2014) find an average value of ~0.07 pc, Kirk et al. (2015) find a range of values from 0.006 pc to 0.03 pc depending on whether the simulation included magnetic fields or not. Federrath (2016) find R_flat values of 0.05 to 0.1 for non-magnetic and magnetic simulations respectively, however, note that their fits are restricted to the inner 0.2 pc of the filament profile. Looking back at Fig. 9 we see that our simulations with the highest Mach values tend to agree with the values found by Smith et al. (2014) and Federrath (2016), but are larger than those reported by Kirk et al. (2015).

Comparing the other two Plummer profile parameters, we see that Smith et al. (2014) find an average power law index (p) of 2.2 while Kirk et al. (2015) quote a range of 1.3−2, both of which are within the range of 1.5–2.5 quoted by Arzoumanian et al. (2011). Federrath (2016) set their power law index value to 2 for all fits. For our simulations, we find a range between ~2 to 4 with an overall average of 3.00. If we split the fitted filaments into those evolving in the SL and CR-only environments, we find average p values of 3.10 and 2.95 respectively. Finally, looking at the central density of the filament, the fits of our simulations yield a maximum value of only 2.55 × 10⁴ cm^-3 which is up to 3 orders of magnitude smaller than the values quoted by Smith et al. (2014) and Kirk et al. (2015) (1.07 × 10⁵ cm^-3 and 10⁴−10⁷ cm^-3, respectively). Federrath (2016) does not give results for this parameter.

Looking back at Table 5, we also note that out of all four studies, our simulations exhibit the poorest resolution while those of Kirk et al. (2015) exhibit the best resolution with the simulations from Smith et al. (2014) and Federrath (2016) falling in between. Therefore, it is not surprising that Kirk et al. (2015) simulations exhibit the smallest FWHM and R_flat while ours exhibit some of the largest. The reason for discrepancies in the fitting parameters can be attributed to the resolution of the simulations, the evolutionary time frame and maximum achieved densities within our simulations. That is, our simulations represent the early phases of cloud/filament formation compared to Smith et al. (2014) and Kirk et al. (2015).

Observational

As discussed at the beginning of Sect. 4.2 there have been several observational studies of various regions with the intent of determining the width of the filaments found within. Figure 10 shows the deconvolved filament width as a function of the central column density for our simulations combined with culled data from surveys of IC5146 (Arzoumanian et al. 2011, diamonds), Galactic Cold Cores (Juvela et al. 2012, pluses) and the DR21 filament in the Cygnus X region (Hennemann et al. 2012, X’s). The coloured symbols depict the particular models as indicated by the legends. For our model data, filled symbols depict SL models while shaded symbols depict CR-only models. Circles indicate models with an initially transcritical mass-to-flux ratio (μ₀ = 1.1) while squares indicate models with an initially super-critical mass-to-flux ratio (μ₀ = 2.0). The solid black line depicts the computed Jeans length as a function of column density [${\mathit{\lambda }}_{\mathit{J}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathit{G}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{\right)}$] for a gas temperature of 10 K. The shaded area within the dashed lines indicate the range of filaments widths in the Aquila and Polaris regions as depicted in Arzoumanian et al. (2011, Fig. 7) for which no data were provided.

From this figure we note a few interesting trends within our model data itself. First, the scatter in the measured widths tends to decrease as the central column density increases although both Models I and III show larger widths even at higher column densities. Comparing the blue coloured symbols, we see again that the filaments forming within supercritical mass-to-flux ratio regions (squares) exhibit smaller widths than those that are evolving in transcritical regions (circles). Finally, we see that the overall width tends to decrease as the degree of turbulence increases. Comparing our data to the expected Jeans length at 10 K (solid line), we see that our simulations tend to exhibit widths that are much larger than the Jeans length, however a few of the super-critical filaments do come close. Comparing our filament widths to the three observed regions, we see that our simulations tend to agree the most with those of the Galactic cold core survey (Juvela et al. 2012), however given the earlier stage of evolution that our simulations represent, no direct conclusion can be drawn.

Looking at the figure as a whole however, we see that there exists a wide range of filament widths and corresponding central densities among the observed regions, particularly within the cold core observations of Juvela et al. (2012). This is to be expected given the fact that these observations were performed over a wide range of points along the Galactic plane and therefore examine varying evolutionary environments between the filaments. The DR21 filament within Cygnus X is thought to be the convergence point of several filaments where high mass star formation is likely occurring (Hennemann et al. 2012). This easily accounts for the significant difference between the densities observed for this region as compared to the others regions and our simulations, all of which are assumed to be undergoing low mass star formation. Finally, the parameter space covered by the shaded area covers filaments within Aquila, IC5146 and Polaris (Arzoumanian et al. 2011) and depicts regions of high, moderate and little to no star formation, respectively. Looking back at the original figure from Arzoumanian et al. (2011, Fig. 7), we see that there seems to be an anti-correlation between the degree of star formation and the filament width, i.e., very active star forming regions show larger filament widths than regions with little to no active star formation.

5. Discussion

The goal of this paper has been to determine the effect of different ionisation profiles on the formation of structures and cores within a turbulent magnetised medium. In addition, we explored the effect of varying the strength of the turbulence and the initial mass-to-flux ratio on the timescale and extent of the structures formed. To this end, we performed two sets of simulations which explored these various effects. The following subsections discuss the main results found through the analysis.

5.1. Effect of ionisation profile

The main effect explored by the simulations in this paper are those of the two ionisation profiles on the formation of structures in a turbulent medium. From our analysis, we found that the ionisation profile does not seem to matter structurally for Mach values greater than or equal to 2. As depicted in Fig. 1, the structures present in Models V and VI are virtually indistinguishable. Upon close inspection, the filaments in Model V (SL ionisation profile) are slightly puffed up compared to those in Model VI (CR-only ionisation profile) however the difference is too small to see without overlaying the two models. We conclude that for highly turbulent driven flows, approximate flux-freezing holds, and there is not enough time for ambipolar diffusion to affect the evolution, even in the CR ionisation cases.

Looking at the density enhancement maps for the final times as shown in Fig. 1, it is easy to see a pattern with increasing Mach number that seems to depend on the choice of ionisation profile. For the SL simulations, it appears that as the degree of turbulence increases, the size of the final structures formed decreases from a monolithic single filament in Model I (v_a/c_s = 0.03) to thinner, more fragmented filaments in Model V (v_a/c_s = 2). Conversely, for the CR-only simulations, the opposite trend is observed i.e., the final size of the structures increases from small core like structures in Model II to long filamentary structures in Model VI. However, when we look at the evolution of each individual model over time, such a pattern does not always hold up. Instead, all models tend to have an initial monolithic filament form at lower density values within which the final density structure forms, whether it be a continuation of a monolithic structure as in Model I, or thinner, more fragmented filamentary structures as observed in Models III–VI. In addition to the differences in the density enhancement structure, we also looked at the mass-to-flux ratio structure. Comparing the two ionisation profiles, we found that SL simulations show very little deviation from the initial mass-to-flux ratio while CR simulations allow for a more significant increase within the medium. Specifically, the nearly flux frozen state of the low density gas in SL simulations results in an amplification in the magnetic field enhancement on the outskirts of the core in comparison to the CR simulation. Although the simulations show that the different ionisation profiles do not imply structural differences for models with Mach number greater than two or a supercritical mass-to-flux ratio, this may not hold true within the kinematics of the region. These possible differences will be explored in the second paper of this series.

5.2. Effect of mach number

The models in Set 1 (see Table 1) explore the effect of changing the strength of the injected turbulence via the Mach number. Here we found that as the turbulent Mach number increases, the structure within the mass-to-flux ratio maps becomes more complicated. In addition, we find that the value of the Mach number has an effect on the physical properties of the structures themselves. Looking at the mass analysis, we find that as the Mach number increases, the mass contained within a particular structure decreases. Comparing corresponding models for the two ionisation profiles, we see that for Mach values greater than 2, there is not a significant difference between the derived masses of the structures.

The fitting analysis tends to agree with all of the effects mentioned above. Specifically, the trends in the FWHM and R_flat parameters shown in Fig. 9 agree with the assertion that there is little structural difference between the SL and CR simulations once the Mach number rises above 2. The fitting also shows that as the Mach number increases, the central density of the filaments also increases, however one must remember that the value of this parameter is solely a result of the fitting of Eq. (10). Based on our fitting analysis, we conclude that the values of the fitted parameters could give some insight into the degree of turbulence within a region assuming that the magnetic field strength (or mass-to-flux ratio) is comparable to the value used in our simulations. Kinematics studies will be needed to better quantify the consistency of our model predictions with observations.

5.3. Effect of initial mass-to-flux ratio

The models in Set 2 in combination with Models V and VI from Set 1 (see Table 1) explore the effect of changing the initial mass-to-flux ratio on the formation of structures within a turbulent medium. Through these eight models, we were able to look at both a highly and moderately subcritical medium, a transcritical medium and supercritical medium evolving within either a SL or CR-only ionisation profile. The main effect of the different initial mass-to-flux ratio observed within these simulations is on the structures formed. For the subcritical simulations (Models XI to XIV), the SL ionisation profile simulations are not able to form a self-gravitating core, even with the added turbulence. Conversely, the CR-only models are able to overcome the strong magnetic field associated with the subcritical mass-to flux ratio. This is a direct result of the low density regions of the SL simulations being held in a nearly flux frozen state within a strong magnetic field, while the ambipolar diffusion within the CR-only simulation still allows for the neutral particles to slip past the field lines even with a strong magnetic field. From these simulations, we found that the main effect of a higher mass-to-flux ratio is to form thinner filaments. This is evident in the filament profile fits performed in Sect. 4.2.

Based on all of the results presented in this subsection and the previous two, we can say that the underlying ionisation and turbulent content can affect the density structure when turbulence is below Mach 2. In a future paper, the velocity field in the different models will be analysed by using the simulations to create synthetic spectra of optically thin low- and high-density tracers. In this way, we can quantify whether different ionisation profiles (i.e., different chemistries) can result in observable differences in clouds with apparently similar structures which can only be observed in dust continuum emission.

5.4. Filaments and fibres

Observations of the Taurus L1495/B213 complex have shown the existence of filaments both with and without core formation (Hacar et al. 2013). Comparing to these observations, we see that the density structure of the filaments is best reproduced by our models with high mass-to-flux ratio (μ₀ = 2.0) and/or at least mildly supersonic turbulence. These observations have also shown the existence of fibres within the macro filamentary structure (Hacar et al. 2013). They also found that fibres could either be fertile or sterile, that is they either show evidence of core formation or they don’t. Based on these observations, two possible scenarios for the formation of such structures have been proposed. First, Tafalla & Hacar (2015) propose a “Fray and Fragment” scenario where a macro filament forms then frays into smaller fibres which then fragment into cores. The second scenario, “Fray and Gather”, is based upon turbulent hydrodynamic simulations from Smith et al. (2016). In their simulations, they found that the medium tends to fray into fibres rather than forming a macro filament. These fibres are then gathered into bundles by the large scale motions of the medium. Looking at the temporal evolution of our simulations (not shown), we find that they tend to support the “Fray and Gather” scenario, however there are a few exceptions, specifically Models I and II. In Model I, as discussed before, the medium undergoes monolithic collapse into a single large filament, which would be the first step of the “Fray and Fragment” scenario, however we do not see any evidence of fibres or multiple cores within the filament. Looking at the evolution of Model II, we see that it forms tenuous, large scale filaments at low densities which fragment into core-like structures with no evidence of fibres, thus skipping the “Fray” stage altogether.

An interesting and somewhat unexpected result from the second set of simulations come from the two subcritical SL models (Models XI and XIII). Although these simulations did not collapse into high density filaments or cores, the evolution of the medium over the course of the simulation shows the formation of thin transient filamentary structures that are similar to those observed by Hacar et al. (2013). Looking at the evolution of these simulations over time (not shown), we see that they form thin filaments that are gathered together by the bulk motions, as in the “Fray and Gather” scenario of Smith et al. (2016), however these gatherings are transient, thus not allowing for the formation of cores or chains of cores to develop along these fibres. This could be the reason behind the fertile and sterile fibres as defined by Hacar et al. (2013). Regions with fertile fibres likely indicate a trans- or supercritical mass-to-flux ratio within the region while the sterile fibres are likely subcritical in nature. In addition, the evolution of our subcritical simulations indicate that these sterile fibre bundles are likely transient in nature; forming and dispersing over the course of a couple Myr. In the following kinematics paper, we will study these spatially resolved fibres to see if they will show up as different kinematic features when observed along the plane of the sheet.

6. Summary

Based on the preceding analysis and discussion, we summarise several of the results below:

• The addition of turbulence within the simulations allows for the formation of filaments within the medium. The degree of turbulence required to form filamentary structure within the simulations depends on the ionisation profile of the medium. The step-like profile requires very minimal turbulence (i.e., Mach 0.03; Model I) while the CR-only simulations require larger degrees of turbulence (e.g., Mach ~1; Model IV).

• Above Mach 2, the effect of the ionisation profile on the density structures formed becomes minimal. This suggests that above a certain turbulent threshold, ambipolar diffusion is not occurring sufficiently to make a dynamical difference in the magnetic forces.

• Analysis of the mass-to-flux ratio maps reveals a patchy distribution of high and low values within the CR-only simulations. The patchiness was found to depend on the degree of turbulence with higher turbulent velocities resulting in finer structures.

• Analysis of the mass-to-flux ratio structures in the vicinity of the cores reveal a unique structure with high values that coincide with the high density regions and significantly lower values just outside. This structure is elongated along the direction of the filament. Profiles across this feature are sharply peaked for the SL ionisation profile which would provide a unique signature within star forming regions. The profiles found in CR-only simulations show the same peak at the high density core however due to the lower ionisation fraction, the profiles are not as sharp as in the SL simulations.

• Variation of the initial mass-to-flux ratio shows that the CR-only simulations are able to form self-gravitating cores within a subcritical (μ < 1) region while SL simulations are unable to, even with the addition of turbulence. The supercritical (μ > 1) simulations yield thinner filaments than the equivalent simulation with a transcritical mass-to-flux ratio (μ ≳ 1).

• Observed filaments of Hacar et al. (2013) and Tafalla & Hacar (2015) are best reproduced by models with supercritical mass-to-flux ratios and/or at least mildly supersonic turbulence.

• Analysis of the temporal evolution within the subcritical SL simulations (Models XI and XIII) showed the formation of thin filamentary fibres within the quiescent medium that are similar to those observed by Hacar et al. (2013). Analysis of the formation of these fibres within our simulations tend to favour the “Fray and Gather” scenario of Smith et al. (2016) as opposed to the “Fray and Fragment” scenario proposed by Tafalla & Hacar (2015).

• Our simulations suggest that fertile filaments observed by Hacar et al. (2013) may exist in regions with super- or transcritical mass-to-flux ratios. Conversely, our simulations suggest that the sterile filaments observed by Hacar et al. (2013) may exist in a highly magnetised medium and that fibre bundles may be transient.

• Filament and core masses are highly dependent on the physical properties of the medium. In general, the CR-only simulations yield filaments with smaller masses than the SL counterparts. This is particularly evident below Mach 2. The overall mass found within the filaments decreases as the degree of turbulence and/or mass-to-flux ratio increases. This is the same for the cores embedded within the filaments.

• Filament fitting showed trends as a function of turbulent velocity and mass-to-flux ratio. Specifically, the Gaussian FWHM and Plummer R_flat parameters show a decreasing trend with increasing Mach number. Conversely, the central density of the filament as determined through the Plummer fits showed an increasing trend with increasing Mach number. In addition, a larger mass-to-flux ratio seems to have the same effect as an increase in turbulence for the three parameters that showed a trend. Finally, there was no discernible trend between the power law index from the Plummer fits and the degree of turbulence.

• Our filaments tend to overlap best with the Cold Core observations (Juvela et al. 2012) while exhibiting lower densities and larger widths when compared to Hennemann et al. (2012) and Arzoumanian et al. (2011) respectively, however simulations with larger degrees of turbulence and/or supercritical mass-to-flux ratios do overlap with the parameter space occupied by the parameter space defined in Arzoumanian et al. (2011).

Although the simulations tend to show little difference between the resulting structures formed in highly turbulent media based on the density enhancement maps, there are differences between the two ionisation profiles when other physical parameters are analysed (i.e., mass-to-flux ratio structures). As such, it is not unrealistic to expect differences to appear in other physical parameters such as the velocity structure. Analysis of such differences/similarities between the evolution of the velocity and resulting synthetic spectra for varying initial turbulent Mach numbers, initial mass-to-flux ratio values and ionisation profiles will be explored in the second paper of this series.

---

Acknowledgments

S.B. was supported by a Natural Science and Engineering Research Council (NSERC) Discovery Grant. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ ERC Grant Agreement No. 320620-PALs.

References

Appendix A: Resolution effects

As indicated in Table 5 and studies such as those presented in Ntormousi et al. (2016), the resolution of the simulations can have an effect on the sizes and scale of structures produced. As such, in this Appendix, we explore the effect of resolution on our fiducial pair of simulations, namely Models V and VI. We have chosen these models as they represent the threshold degree of turbulence and mass-to-flux ratio required to produce filaments. We look at both Models V and VI to discern if either of the underlying ionisation profiles (SL or CR-only, respectively) react differently at higher resolution. For the high resolution simulations, we have doubled the number of pixels in each dimension (1024 × 1024) while maintaining the same physical size of the box. This results in a pixel size that is half that of the low resolution models (i.e., 0.0148 pc).

The first difference evident between the high and low resolution simulations is the run time. The original simulations for Models V and VI ran for 3.8 and 3.7 dimensionless time units, respectively, while the high resolution simulations ran until dimensionless times of 2.4 and 2.2, respectively. Figure A.1 shows a comparison of the column density enhancement maps for the low (top row) and high (bottom row) resolution simulations of Models V (left column) and VI (right column). At first look, the high and low resolution maps look different. This, however, is not surprising as there is no such thing as an identical random number field when looking at models with different resolutions, even if we use the same seed i.e., increasing the resolution increases the number of pixels and changes the arrangement of the applied random number field.

Focussing on the high resolution simulations, we see that, just as in the low resolution versions, both the SL and CR-only ionisation profiles produce similar structures. In terms of filament formation, we see that the high resolution simulations only produce one filament which conforms to our definition, namely structures with A_V> 2 however we can see that there are other regions that are in the process of forming filaments at the run termination time. High resolution runs of these models with increased stopping conditions (i.e., σ/σ_n,0 = 20) show the continued formation of these filaments, however for comparison purposes, we will focus on the simulations with the same stopping conditions (i.e., σ/σ_n,0 = 10). As such there is only one high density region within the high resolution simulations which we can use to compare to the low resolution runs. A further set of high/low resolution runs for Model V using different seeds for the turbulent velocity distribution (not shown) yielded simulations which ran to dimensionless times of 5.8 and 6.4, respectively and resulted in the formation of multiple filaments within the high resolution simulation.

Figure A.2 shows the mass-to-flux ratio profiles for the densest cores in the high (solid lines) and low (dashed lines) resolution simulations. The black curves depict the SL models (V and V_h) while the red curves depict the CR models (VI and VI_h). Here we are comparing the only high density region within the high resolution simulation to the highest density regions within the low resolution simulations, i.e., filaments V-4 and VI-4 for Models V and VI, respectively. The curves have all been shifted such that the peaks of the profiles coincide at y = 0. The x-axis depicts the distance away from the central peak in parsecs. Here, we see that the high resolution simulations again exhibit the same profile described in Sect. 3.2.2, however, the increased resolution allows for an increased degree of detail. Specifically, the peaked profile becomes more distinct in the high resolution simulations with the CR-only model again showing a wider feature with more variability in the wings compared to the SL model.

Finally, to further compare the two resolutions, we also extracted the average profile of the high resolution filaments. As stated earlier, when assuming the same stopping conditions as the low resolution simulations, the high resolution runs only produced one filament which adheres to such criteria. Table A.1 shows a comparison between the FWHM of the single high resolution filaments and the low resolution filaments originally presented in Table 4 for both ionisation profiles. For both models, we see that the FWHM of the high resolution simulation falls within the range of FWHMs set by the low resolution filaments. For Model VI, we also see that the high resolution filament yields a very similar value to half of the low resolution filaments. When comparing the high resolution filament directly to the low resolution filament which contains the highest density core, we see that the high resolution simulations results in larger FWHM by 21% and 32% for Models V and VI, respectively. Conversely, when compared to the highest FWHM within the low resolution set, the high resolution simulations differ by 59% and 44% for Models V and VI, respectively. Therefore, we conclude that resolution does have an effect on the measured widths of the filaments, however it does not result in smaller filaments and on a whole the FWHMs of the high resolution filaments agree with the range of FWHMs found in the low resolution simulations.

All Tables

All Figures