# Smoothed Particle Magnetohydrodynamics Simulations of Protostellar Jets and Turbulent Dynamos

## Abstract

We presents results from Smoothed Particle Magnetohydrodynamics simulations of collapsing molecular cloud cores, and dynamo amplification of the magnetic field in the presence of Mach 10 magnetised turbulence. Our star formation simulations have produced, for the first time ever, highly collimated magnetised protostellar jets from the first hydrostatic core phase. Up to 40% of the initial core mass may be ejected through this outflow. The primary difficulty in performing these simulations is maintaining the divergence free constraint of the magnetic field, and to address this issue, we have developed a new divergence cleaning method which has allowed us to stably follow the evolution of these protostellar jets for long periods. The simulations performed of supersonic MHD turbulence are able to exponentially amplify magnetic energy by up to 10 orders of magnitude via turbulent dynamo. To reduce numerical dissipation, a new shock detection algorithm is utilised which is able to track magnetic shocks throughout a large range of magnetic field strengths.

2013 \SetConfTitleMagnetic Fields in the Universe IV \listofauthorsT. S. Tricco, D. J. Price, C. Federrath, & M. R. Bate \indexauthorTricco, T. S. \indexauthorPrice, D. J. \indexauthorFederrath, C. \indexauthorBate, M. R. \resumen\addkeywordstars: formation \addkeywordISM: Jets and outflows \addkeywordISM: magnetic fields \addkeywordturbulence \addkeywordmethods: numerical

## 1 Introduction

Magnetic fields play an important role in all phases of star formation. On the scale of molecular clouds, the sites where stars are formed, magnetic fields can affect the structure of the supersonic turbulence. It is known that this turbulence drives the local compression needed to trigger gravitational collapse of molecular cloud cores, and simulations by Padoan & Nordlund (2011) and Federrath & Klessen (2012) show that the additional pressure from magnetic fields can help support against gravitational collapse and reduce star formation rates by factors of 2–3.

On the scale of individual protostars, magnetic fields are responsible for launching protostellar jets and outflows. These jets and outflows are both well observed (Richer et al., 2000; Wu et al., 2004) and simulated (e.g. Machida et al., 2008; Commerçon et al., 2010; Price et al., 2012; Machida & Hosokawa, 2013). They reduce the efficiency of star formation (Matzner & McKee, 2000; Hansen et al., 2012), and are drivers of turbulence in the interstellar medium (Nakamura & Li, 2007; Carroll et al., 2010).

In this work, we have performed Smoothed Particle Magnetohydrodynamic (SPMHD) (Price & Monaghan, 2004a, b, 2005) simulations of the gravitational collapse of a prestellar core to form the first hydrostatic core (before the protostar is formed), and also of magnetised Mach 10 turbulence which exponentially amplifies an initial seed magnetic field. These simulations use new numerical techniques to maintain the divergence free constraint on the magnetic field, detect shocks, and reduce numerical dissipation of the magnetic field.

This paper begins with a brief review of SPMHD in §2, discussing how magnetic fields are simulated with SPMHD, and the benefits of using SPMHD for star formation simulations. In §3, the constrained hyperbolic divergence cleaning method is introduced, which is used to maintain . In §4, we present results from simulations of a collapsing prestellar core. The jet produced in our simulations during protostellar collapse is discussed in §4.2. In §5, results from simulations of Mach 10 magnetised turbulence are presented. These simulations focus on the dynamo amplification of the magnetic field, with comparison to results from grid based methods. A summary and discussion is given in §6.

## 2 Smoothed Particle Magnetohydrodynamics

Smoothed Particle Hydrodynamics (SPH) is a numerical method for simulating fluid flow (see reviews by Monaghan 2005; Price 2012). The hydrodynamic equations are solved by discretising the fluid into a set of particles which contain a portion of the mass, energy, and momentum of the fluid. Fluid quantities, such as density, are calculated per particle by interpolating from neighbouring particles using a kernel weighted summation.

SPH is widely used in astrophysics. It can easily handle complex geometries, has excellent conservation properties, and couples easily with N-body methods for gravity. Regions of higher density contain more mass, and therefore more resolution elements. This is useful for star formation, because as gas collapses to form dense objects, the mass, and hence particles, trace the collapse providing continuous resolution of that process until the Jeans mass falls below the mass resolution (Bate & Burkert, 1997).

### 2.1 The equations of SPMHD

The ideal MHD equations solved with SPMHD are given, for a particle , by

(1) | ||||

(2) | ||||

(3) | ||||

(4) |

where and are the velocity and magnetic fields, denotes , and is an interpolation kernel. We use the cubic spline kernel (Monaghan & Lattanzio, 1985) in this work.

Variable resolution is obtained by self-consistently deriving the density, , and smoothing length, , through iteration of equations (1) and (2). The smoothing length is related to the local particle spacing by , with corresponding to the number of dimensions. Variable smoothing length gradients are handled by (see Springel & Hernquist 2002).

The momentum equation (3) is derived from the Lagrangian and exactly conserves energy and momentum. It is expressed in terms of the Maxwell stress tensor,

(5) |

with the thermal pressure, , obtained through a suitable equation of state. Terms containing are subtracted from the SPMHD momentum equation, introducing a small amount of non-conservation of energy and momentum, but greatly enhancing stability and performance. The induction equation (4) is a representation of .

### 2.2 Shock capturing

Hydrodynamic and magnetic shocks are captured by the addition of artificial viscosity and resistivity to the momentum and induction equations. The artificial viscosity used here was formulated by Monaghan (1997) by analogy to Riemann solvers, and is given by

(6) |

The signal velocity represents the characteristic speed of information propagation across the shock, given by

(7) |

For ideal MHD, the sound speed, , is replaced by the fast MHD wave speed,

(8) |

where corresponds to the Alfvén speed.

Artificial resistivity was formulated through a similar procedure by Price & Monaghan (2005). The corresponding term in the induction equation (4) is given by,

(9) |

The signal velocity for artificial resistivity is chosen as , which is the averaged fast MHD wave speeds.

The dimensionless parameters and are of order unity. To reduce the dissipation from artificial viscosity and resistivity away from shocks (where it is unnecessary), and may be set individual for each particle and used to regulate the strength of the applied dissipation. Morris & Monaghan (1997) proposed integrating according to

(10) |

with . This equation increases in regions of converging flow, with a post-shock decay timescale, , of approximately five smoothing lengths ().

Price & Monaghan (2005) created a similar switch for artificial resistivity, using

(11) |

with a range .

Recently, we have proposed a new switch for artificial resistivity that is more robust at detecting shocks and leads to less overall dissipation (Tricco & Price, submitted). It sets

(12) |

in the range . This increases artificial resistivity in regions of strong magnetic field gradients. Since , this leads to a quantity which is related to the Alfvénic Mach number. By normalising the gradient against the magnitude of the magnetic field, the switch responds to the relative degree of discontinuity and does so independently of the absolute magnetic field strength (this is important for the dynamo amplification simulations in §5). Setting the value of directly in this manner improves the responsiveness to shocks by removing the time delay present in equation (11).

## 3 Constrained hyperbolic divergence cleaning

The zero divergence constraint on the magnetic field is maintained using constrained hyperbolic divergence cleaning (Tricco & Price, 2012). The cleaning algorithm couples an additional scalar field, , to the magnetic field, and divergence error in the magnetic field is dispersed and diffused using a series of damped waves. By spreading the divergence error over a larger volume, it is able to be removed faster than using just a diffusion term alone, and furthermore, the impact of any single large source of error is reduced.

The method is a Hamiltonian version of hyperbolic divergence cleaning (Dedner et al., 2002) which has been derived by defining the energy content of the field, and including it as part of the Lagrangian. This leads to continuum equations,

(13) | ||||

(14) |

and SPMHD equations

(15) | ||||

(16) |

The divergence wave speed, , is chosen as the maximum allowable according to the Courant timestep criterion (typically this would be the fast MHD wave speed). The damping term, , is best chosen in the regime of critical damping, which from empirical tests (Tricco & Price, 2012), is for 2D and for 3D. Equation (14) differs from Dedner et al. (2002) by the addition of the term, which describes how changes as the fluid is expanded or compressed. The constraint from energy conservation also leads to a specific choice of derivative operators in the SPH formulation, as given by equations (15) and (16).

By building the method from the ground up in the context of the Lagrangian equations of motion, it inherently retains the stability properties of SPH and is guaranteed to always decrease the divergence of the magnetic field. This fixes problems in the previous implementation by Price & Monaghan (2005) particularly at density contrasts and free boundaries.

## 4 Molecular cloud core collapse

Larson (1969) showed collapsing prestellar cores undergo a two stage collapse for low mass star formation. For both stages the collapse is nearly isothermal, but has a brief adiabatic phase when molecular hydrogen reaches densities higher than g . At these densities, the gas becomes optically thick, trapping radiation. These first hydrostatic core objects have a lifetime of only several thousand years, ending when the gas reaches K and the molecular hydrogen disassociates. The radiation is then able to freely escape, and the core undergoes a second near-isothermal collapse to form a protostellar core.

Observations of first core objects have been difficult because of the low luminosity of these objects, and it is only within recent years that candidate detections have been made. Pineda et al. (2011) observed a low-mass dense core in the Perseus Molecular Cloud, with upper limits on bolometric luminosity and temperature of L and K, and found traces of a slow ( km ), poorly collimated outflow. Per Bolo 58 has been studied by Enoch et al. (2010), finding it to be a promising first hydrostatic core candidate with an internal luminosity of . Dunham et al. (2011) detected a well collimated () bipolar outflow in Per Bolo 58 with characteristic velocity of km .

We have performed simulations of a collapsing prestellar core during the first stage of collapse to form the first hydrostatic core (Price, Tricco, & Bate, 2012).

### 4.1 Initial conditions and numerical details

The simulations are performed for a spherical core with radius cm ( AU), giving an initial density of g . It is set in solid body rotation with angular velocity rad . The free fall time is yr. A barotropic equation of state is used (as described in Price et al. 2012), where the gas is isothermal below a critical density of g , and adiabatic above this density. The speed of sound is cm .

The initial magnetic field is uniform along the rotation axis with mass-to-flux ratio 5, or . Edge effects with the magnetic field are avoided by embedding the core in an ambient medium in a periodic box of length . The medium has a density contrast of 1:30 and is set in pressure equilibrium with the core.

The core is simulated using particles. Self-gravity is included by use of a hierarchical binary tree (Benz et al., 1990), with the SPH smoothing kernel used for gravitational force softening (Price & Monaghan, 2007). A sink particle (Bate et al., 1995) is inserted once the density reaches g and accretes material within AU.

Only a minimal amount of artificial resistivity is applied to the magnetic field, using the switch described in equation (11) in the range . The constrained hyperbolic divergence cleaning algorithm (§3) is used to remove errors arising from the divergence of the magnetic field. Using this algorithm, the average divergence error of the field is kept to within (Figure 1).

### 4.2 First core jet

The simulations are performed until . At , the winding of the magnetic field by the infalling material launches a well collimated jet along the axis of rotation (Figure 2). The jet has mean velocity km , with top end velocities of 5–7 km (see Figure 3), which is consistent with observed speeds and the escape velocity for an object of this mass and radius.

The jet is highly efficient at removing mass. By , 40% of the original material in the core has been ejected through this outflow, consistent with outflow studies by Matzner & McKee (2000) and Hansen et al. (2012). The material around the sink particle settles into a disc-like object due to the conservation of angular momentum, but has sub-Keplerian orbital velocites by a factor of 3–4.

The magnetic field near the sink particle is mG, wound in a toroidal geometry. This leads to the “wiggles” observed in the jet as a result of the expanding magnetic field in the -direction. The magnetic field lines at are shown in Figure 4, where the field from each particle is represented with an opacity proportional to field strength.

## 5 Mach 10 magnetised turbulence

Why are observed magnetic fields in the universe as strong as they are? It is becoming increasingly understood that small scale turbulent dynamos drive exponential amplification of the magnetic field, such that even if initial magnetic fields had tiny field strengths, they would have rapidly reached observed values (see reviews by Beck et al. 1996; Brandenburg & Subramanian 2005; Widrow et al. 2012).

We have been performing a code comparison between SPMHD and finite difference methods on the dynamo amplification of magnetic fields in Mach 10 turbulence, representative of turbulence found in molecular clouds in the Milky Way (for reviews, see Evans 1999; Elmegreen & Scalo 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007). These simulations begin with an initially weak magnetic field, which is exponentially amplified through the conversion of turbulent energy until the magnetic energy reaches equipartition with the kinetic energy (an increase of 10 orders of magnitude in this case). Our comparison extends the purely hydrodynamic code comparison by Price & Federrath (2010) on the statistics of driven Mach 10 turbulence between SPH and grid-based methods. Their conclusion was that excellent agreement was found in the properties of the turbulence, with SPH performing better at resolving dense structures, and grid-based methods better for volumetric quantities.

### 5.1 Initial conditions

The simulations are run at resolutions of and particles. The initial conditions are similar to the parameter study performed by Federrath et al. (2011). The density is uniform, with , and an isothermal equation of state is used with . The initial magnetic field is such that the initial plasma beta is .

The turbulence is driven and sustained using an acceleration based upon the Ornstein-Uhlenbeck process (Eswaran & Pope, 1988; Federrath et al., 2010), which is a stochastic process with a finite autocorrelation timescale that drives motion at low wave numbers. The driving force is constructed in Fourier space, allowing it to be decomposed into solenoidal and compressive components and for this case we only use the solenoidal component.

### 5.2 Detecting magnetic shocks

It is important to correctly capture shocks in the magnetic field for these simulations. However, the dissipation from the added artificial resistivity reduces the magnetic Reynolds number. Since molecular clouds are known to have high kinetic and magnetic Reynolds numbers (–), it is important to reduce sources of numerical dissipation. Thus, a switch is used to “turn off” artificial resistivity in regions away from shocks. We found the switch proposed by Price & Monaghan (2005) was not able to detect shocks while the magnetic field was weak, leading to significant noise in the magnetic field and spurious growth rates.

We have developed a new artificial resistivity switch (Tricco & Price, submitted) for this work, as described in §2.2. It sets , which measures the relative degree of discontinuity in the magnetic field. In this way, it is able to detect and capture shocks throughout the several orders of magnitude change in magnetic field strength.

An important note about dissipation terms in SPMHD is that they are resolution dependent. The dissipation can be effectively halved by doubling the resolution of the simulation.

### 5.3 Turbulence results

The turbulence is simulated for 60 turbulent turnover times using the Flash code (Fryxell et al., 2000; Dubey et al., 2008) at grid cells, and with the Phantom SPMHD code at and particles. The SPMHD simulations have been run for both artificial resistivity applied with a fixed parameter, and using the new artificial resistivity switch to reduce dissipation.

Renderings of the evolution of the -integrated column density and magnetic field are shown in Figure 5 for the particle simulation using the new artificial resistivity switch. The structure in the magnetic field closely resembles the shock structures in the density field, which is to be expected since the magnetic field is weak. Even once the field is reaching saturation, shocks in the magnetic field are still driven primarily by the forcing though subtle differences start to become apparent.

The growth of magnetic energy as a function of time is shown in Figure 6. Comparable growth rates for grid cells can be achieved using either particles when applying fixed artificial resistivity everywhere, or at the same resolution of particles when using the new resistivity switch. Thus it can be concluded that using the new resistivity switch produces an effect similar to doubling the resolution.

The saturation level is in agreement between the Flash results and the SPMHD results when using a fixed artificial resistivity parameter. When using the new artificial resistivity switch to reduce magnetic dissipation, however, the saturation level for the SPMHD results are 2–3 higher for both the and simulations.

## 6 Summary and discussion

In this work, we have performed Smoothed Particle Magnetohydrodynamic (SPMHD) simulations of the gravitational collapse of a prestellar core to form the first hydrostatic core, and also of magnetised Mach 10 turbulence, representative of conditions in molecular clouds. These simulations used the constrained hyperbolic divergence cleaning method (Tricco & Price, 2012) to maintain the constraint on the magnetic field. It is a Hamiltonian version of hyperbolic divergence cleaning (Dedner et al., 2002) that we have formulated for SPMHD through derivation from the discretised Lagrangian. Thus it possesses the conservation and stability properties inherent to SPH, and was found to reduce errors in the magnetic field by 10.

The collapsing core simulations (Price et al., 2012) produced a slow, well collimated jet from the central object that has properties consistent with candidate observations of first hydrostatic cores. It is efficient at removing mass out of the core, with up to 40% of the material being removed by the time the remaining mass in the core has been accreted onto the central sink particle.

Our simulations of Mach 10 magnetised turbulence model dynamo amplification of the magnetic field through the conversion of turbulent energy into magnetic energy. An initially weak magnetic field is present which was exponentially increased 10 orders of magnitude in energy until it reaches saturation. Both the saturation level and growth rates were consistent with results from the grid based code Flash when using a fixed artificial resistivity parameter. Simulations were also run using a switch for artificial resistivity to reduce numerical dissipation (Tricco & Price, submitted), which lead to similar growth rates as if double the resolution had been used, but with a saturation level 2–3 higher.

T. Tricco thanks the conference organisers for the opportunity to speak. T. Tricco is supported by Endeavour IPRS and APA postgraduate research scholarships. We are grateful for funding via Australian Research Council Discovery Projects grants DP1094585 and DP110102191. This research was undertaken with the assistance of resources provided at the Multi-modal Australian ScienceS Imaging and Visualisation Environment (MASSIVE) through the National Computational Merit Allocation Scheme supported by the Australian Government. The Flash simulations were run at LRZ (grant pr32lo) and JSC (grant hhd20).

### Footnotes

- affiliation: Monash Centre for Astrophysics, School of Mathematical Sciences, Monash University, Vic, 3800, Australia (terrence.tricco@monash.edu, daniel.price@monash.edu, christoph.federrath@monash.edu).
- affiliation: Monash Centre for Astrophysics, School of Mathematical Sciences, Monash University, Vic, 3800, Australia (terrence.tricco@monash.edu, daniel.price@monash.edu, christoph.federrath@monash.edu).
- affiliation: Monash Centre for Astrophysics, School of Mathematical Sciences, Monash University, Vic, 3800, Australia (terrence.tricco@monash.edu, daniel.price@monash.edu, christoph.federrath@monash.edu).
- affiliation: School of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, EX4 4QL, United Kingdom (mbate@astro.ex.ac.uk)

### References

- Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362
- Bate, M. R. & Burkert, A. 1997, MNRAS, 288, 1060
- Beck, R., Brandenburg, A., Moss, D., Shukurov, A., & Sokoloff, D. 1996, ARA&A, 34, 155
- Benz, W., Cameron, A. G. W., Press, W. H., & Bowers, R. L. 1990, ApJ, 348, 647
- Brandenburg, A. & Subramanian, K. 2005, Phys. Rep., 417, 1
- Carroll, J. J., Frank, A., & Blackman, E. G. 2010, ApJ, 722, 145
- Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3
- Dedner, A., Kemm, F., Kröner, D., Munz, C.-D., Schnitzer, T., & Wesenberg, M. 2002, J. Comput. Phys., 175, 645
- Dubey, A., Fisher, R., Graziani, C., Jordan, IV, G. C., Lamb, D. Q., Reid, L. B., Rich, P., Sheeler, D., Townsley, D., & Weide, K. Astronomical Society of the Pacific Conference Series, Vol. 385, , Numerical Modeling of Space Plasma Flows, ed. N. V. PogorelovE. Audit & G. P. Zank, 145
- Dunham, M. M., Chen, X., Arce, H. G., Bourke, T. L., Schnee, S., & Enoch, M. L. 2011, ApJ, 742, 1
- Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 211
- Enoch, M. L., Lee, J.-E., Harvey, P., Dunham, M. M., & Schnee, S. 2010, ApJL, 722, L33
- Eswaran, V. & Pope, S. B. 1988, Computers and Fluids, 16, 257
- Evans, II, N. J. 1999, ARA&A, 37, 311
- Federrath, C., Chabrier, G., Schober, J., Banerjee, R., Klessen, R. S., & Schleicher, D. R. G. 2011, Physical Review Letters, 107, 114504
- Federrath, C. & Klessen, R. S. 2012, ApJ, 761, 156
- Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.-M. 2010, A&A, 512, A81
- Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W., & Tufo, H. 2000, ApJS, 131, 273
- Hansen, C. E., Klein, R. I., McKee, C. F., & Fisher, R. T. 2012, ApJ, 747, 22
- Larson, R. B. 1969, MNRAS, 145, 271
- Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
- Machida, M. N. & Hosokawa, T. 2013, MNRAS, 431, 1719
- Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2008, ApJ, 676, 1088
- Matzner, C. D. & McKee, C. F. 2000, ApJ, 545, 364
- McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565
- Monaghan, J. J. 1997, J. Comput. Phys., 136, 298
- —. 2005, Reports on Progress in Physics, 68, 1703
- Monaghan, J. J. & Lattanzio, J. C. 1985, A&A, 149, 135
- Morris, J. P. & Monaghan, J. J. 1997, J. Comput. Phys., 136, 41
- Nakamura, F. & Li, Z.-Y. 2007, ApJ, 662, 395
- Padoan, P. & Nordlund, Å. 2011, ApJ, 730, 40
- Pineda, J. E., Arce, H. G., Schnee, S., Goodman, A. A., Bourke, T., Foster, J. B., Robitaille, T., Tanner, J., Kauffmann, J., Tafalla, M., Caselli, P., & Anglada, G. 2011, ApJ, 743, 201
- Price, D. J. 2012, J. Comput. Phys., 231, 759
- Price, D. J. & Federrath, C. 2010, MNRAS, 406, 1659
- Price, D. J. & Monaghan, J. J. 2004a, MNRAS, 348, 123
- —. 2004b, MNRAS, 348, 139
- —. 2005, MNRAS, 364, 384
- —. 2007, MNRAS, 374, 1347
- Price, D. J., Tricco, T. S., & Bate, M. R. 2012, MNRAS, 423, L45
- Richer, J. S., Shepherd, D. S., Cabrit, S., Bachiller, R., & Churchwell, E. 2000, Protostars and Planets IV, 867
- Springel, V. & Hernquist, L. 2002, MNRAS, 333, 649
- Tricco, T. S. & Price, D. J. 2012, J. Comput. Phys., 231, 7214
- —. submitted, MNRAS
- Widrow, L. M., Ryu, D., Schleicher, D. R. G., Subramanian, K., Tsagas, C. G., & Treumann, R. A. 2012, Space Sci. Rev., 166, 37
- Wu, Y., Wei, Y., Zhao, M., Shi, Y., Yu, W., Qin, S., & Huang, M. 2004, A&A, 426, 503