There has been interest in recent years to assess the ability of astrophysical hydrodynamics codes to correctly model the Kelvin-Helmholtz instability. Smoothed particle hydrodynamics (SPH), in particular, has received significant attention, though there has yet to be a clear demonstration that SPH yields converged solutions that are in agreement with other methods. We have performed SPH simulations of the Kelvin-Helmholtz instability using the test problem put forward by Lecoanet et al (2016). We demonstrate that the SPH solutions converge to the reference solution in both the linear and non-linear regimes. Quantitative convergence in the strongly non-linear regime is achieved by using a physical Navier-Stokes viscosity and thermal conductivity. We conclude that standard SPH with an artificial viscosity can correctly capture the Kelvin-Helmholtz instability.

We perform simulations of the Kelvin-Helmholtz instability using smoothed particle hydrodynamics (SPH). The instability is studied both in the linear and strongly non-linear regimes. The smooth, well-posed initial conditions of Lecoanet et al. (2016) are used, along with an explicit Navier-Stokes viscosity and thermal conductivity to enforce the evolution in the non-linear regime. We demonstrate convergence to the reference solution using SPH. The evolution of the vortex structures and the degree of mixing, as measured by a passive scalar `colour' field, match the reference solution. Tests with an initial density contrast produce the correct qualitative behaviour. The L2 error of the SPH calculations decreases as the resolution is increased. The primary source of error is numerical dissipation arising from artificial viscosity, and tests with reduced artificial viscosity have reduced L2 error. A high-order smoothing kernel is needed in order to resolve the initial velocity amplitude of the seeded mode and inhibit excitation of spurious modes. We find that standard SPH with an artificial viscosity has no difficulty in correctly modelling the Kelvin-Helmholtz instability and yields convergent solutions.

We present Phantom, a fast, parallel, modular and low-memory smoothed particle hydrodynamics and magnetohydrodynamics code developed over the last decade for astrophysical applications in three dimensions. The code has been developed with a focus on stellar, galactic, planetary and high energy astrophysics and has already been used widely for studies of accretion discs and turbulence, from the birth of planets to how black holes accrete. Here we describe and test the core algorithms as well as modules for magnetohydrodynamics, self-gravity, sink particles, H_2 chemistry, dust-gas mixtures, physical viscosity, external forces including numerous galactic potentials as well as implementations of Lense-Thirring precession, Poynting-Robertson drag and stochastic turbulent driving. Phantom is hereby made publicly available.

In numerical simulations, artificial terms are applied to the evolution equations for stability. To prove their validity, these terms are thoroughly tested in test problems where the results are well known. However, they are seldom tested in production-quality simulations at high resolution where they interact with a plethora of physical and numerical algorithms. We test three artificial resistivities in both the Orszag-Tang vortex and in a star formation simulation. From the Orszag-Tang vortex, the Price et. al. (2017) artificial resistivity is the least dissipative thus captures the density and magnetic features; in the star formation algorithm, each artificial resistivity algorithm interacts differently with the sink particle to produce various results, including gas bubbles, dense discs, and migrating sink particles. The star formation simulations suggest that it is important to rely upon physical resistivity rather than artificial resistivity for convergence.

We perform numerical simulations of dusty, supersonic turbulence in molecular clouds. We model 0.1, 1 and 10 {\mu}m sized dust grains at an initial dust-to-gas mass ratio of 1:100, solving the equations of combined gas and dust dynamics where the dust is coupled to the gas through a drag term. We show that, for 0.1 and 1 {\mu}m grains, the dust-to-gas ratio deviates by typically 10-20% from the mean, since the stopping time of the dust due to gas drag is short compared to the dynamical time. Contrary to previous findings, we find no evidence for orders of magnitude fluctuation in the dust-to-gas ratio for 0.1 {\mu}m grains. Larger, 10 {\mu}m dust grains may have dust-to-gas ratios increased by up to an order of magnitude locally. Both small (0.1 {\mu}m) and large ($\gtrsim$ 1 {\mu}m) grains trace the large-scale morphology of the gas, however we find evidence for 'size-sorting' of grains, where turbulence preferentially concentrates larger grains into dense regions. Size-sorting may help to explain observations of 'coreshine' from dark clouds, and why extinction laws differ along lines of sight through molecular clouds in the Milky Way compared to the diffuse interstellar medium.

We study the evolution of galactic magnetic fields using 3D smoothed particle magnetohydrodynamics (SPMHD) simulations of galaxies with an imposed spiral potential. We consider the appearance of reversals of the field, and amplification of the field. We find magnetic field reversals occur when the velocity jump across the spiral shock is above $\approx$20km s$^{-1}$, occurring where the velocity change is highest, typically at the inner Lindblad resonance (ILR) in our models. Reversals also occur at corotation, where the direction of the velocity field reverses in the co-rotating frame of a spiral arm. They occur earlier with a stronger amplitude spiral potential, and later or not at all with weaker or no spiral arms. The presence of a reversal at a radii of around 4--6 kpc in our fiducial model is consistent with a reversal identified in the Milky Way, though we caution that alternative Galaxy models could give a similar reversal. We find that relatively high resolution, a few million particles in SPMHD, is required to produce consistent behaviour of the magnetic field. Amplification of the magnetic field occurs in the models, and while some may be genuinely attributable to differential rotation or spiral arms, some may be a numerical artefact. We check our results using Athena, finding reversals but less amplification of the field, suggesting that some of the amplification of the field with SPMHD is numerical.

We present an updated constrained hyperbolic/parabolic divergence cleaning algorithm for smoothed particle magnetohydrodynamics (SPMHD) that remains conservative with wave cleaning speeds which vary in space and time. This is accomplished by evolving the quantity $\psi / c_h$ instead of $\psi$. Doing so allows each particle to carry an individual wave cleaning speed, $c_h$, that can evolve in time without needing an explicit prescription for how it should evolve, preventing circumstances which we demonstrate could lead to runaway energy growth related to variable wave cleaning speeds. This modification requires only a minor adjustment to the cleaning equations and is trivial to adopt in existing codes. Finally, we demonstrate that our constrained hyperbolic/parabolic divergence cleaning algorithm, run for a large number of iterations, can reduce the divergence of the field to an arbitrarily small value, achieving $\nabla \cdot B=0$ to machine precision.

Obtaining a stable magnetohydrodynamical (MHD) formalism in SPH - i.e. smoothed particle magnetohydrodynamics (SPMHD) - has proven remarkably difficult. To implement MHD requires two steps: a modification to the momentum equation and an induction equation, and both present challenges. We first provide an overview of how SPMHD is implemented, and then discuss how this implementation fails and the limitation of various corrective methods - with particular reference to the effects of particle disorder. Although there are many problems for which, with careful choice of corrective measures, good results can be obtained, we then show that, at the very limits of the state of the art, the ability to perform stable MHD calculations in SPH is curtailed by numerical issues.

Supersonic turbulence is believed to be at the heart of star formation. We have performed smoothed particle magnetohydrodynamics (SPMHD) simulations of the small-scale dynamo amplification of magnetic fields in supersonic turbulence. The calculations use isothermal gas driven at rms velocity of Mach 10 so that conditions are representative of star-forming molecular clouds in the Milky Way. The growth of magnetic energy is followed for 10 orders in magnitude until it reaches saturation, a few percent of the kinetic energy. The results of our dynamo calculations are compared with results from grid-based methods, finding excellent agreement on their statistics and their qualitative behaviour. The simulations utilise the latest algorithmic developments we have developed, in particular, a new divergence cleaning approach to maintain the solenoidal constraint on the magnetic field and a method to reduce the numerical dissipation of the magnetic shock capturing scheme. We demonstrate that our divergence cleaning method may be used to achieve $\nabla \cdot {\bf B}=0$ to machine precision, albeit at significant computational expense.

We perform a comparison between the smoothed particle magnetohydrodynamics (SPMHD) code, Phantom, and the Eulerian grid-based code, Flash, on the small-scale turbulent dynamo in driven, Mach 10 turbulence. We show, for the first time, that the exponential growth and saturation of an initially weak magnetic field via the small-scale dynamo can be successfully reproduced with SPMHD. The two codes agree on the behaviour of the magnetic energy spectra, the saturation level of magnetic energy, and the distribution of magnetic field strengths during the growth and saturation phases. The main difference is that the dynamo growth rate, and its dependence on resolution, differs between the codes, caused by differences in the numerical dissipation and shock capturing schemes leading to differences in the effective Prandtl number in Phantom and Flash.

Numerical methods to improve the treatment of magnetic fields in smoothed field magnetohydrodynamics (SPMHD) are developed and tested. Chapter 2 is a review of SPMHD. In Chapter 3, a mixed hyperbolic/parabolic scheme is developed which cleans divergence error from the magnetic field. Average divergence error is an order of magnitude lower for all test cases considered, and allows for the stable simulation of the gravitational collapse of magnetised molecular cloud cores. The effectiveness of the cleaning may be improved by explicitly increasing the hyperbolic wave speed or by cycling the cleaning equations between timesteps. In the latter, it is possible to achieve DivB=0. Chapter 4 develops a switch to reduce dissipation of the magnetic field from artificial resistivity. Compared to the existing switch in the literature, this leads to sharper shock profiles in shocktube tests, lower overall dissipation of magnetic energy, and importantly, is able to capture magnetic shocks in the highly super-Alfvenic regime. Chapter 5 compares these numerical methods against grid-based MHD methods (using the Flash code) in simulations of the small-scale dynamo amplification of a magnetic field in driven, isothermal, supersonic turbulence. Both codes exponentially amplify the magnetic energy at a constant rate, though SPMHD shows a resolution dependence that arises from the scaling of the numerical dissipation terms. The time-averaged saturated magnetic spectra have similar shape, and both codes have PDFs of magnetic field strength that are log-normal, which become lopsided as the magnetic field saturates. We conclude that SPMHD is able to reliably simulate the small-scale dynamo amplification of magnetic fields. Chapter 6 concludes the thesis and presents some preliminary work demonstrating that SPMHD can activate the magneto-rotational instability in 2D shearing box tests.

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.

We describe a new switch to reduce dissipation from artificial resistivity in Smoothed Particle Magnetohydrodynamics simulations. The switch utilises the gradient of the magnetic field to detect shocks, setting alpha_B = h |gradB| / |B|. This measures the relative degree of discontinuity, and the switch is not dependent on the absolute field strength. We present results comparing the new resistivity switch to the switch of Price & Monaghan (2005), showing that it is more robust in capturing shocks (especially in weak fields), while leading to less overall dissipation. The design of this switch is generalised to create similar switches for artificial viscosity and thermal conduction, with proof of concept tests conducted on a Sod shock tube and Kelvin-Helmholtz instabilities.

We have performed smoothed particle radiation magnetohydrodynamics (SPRMHD) simulations of the collapse of rotating, magnetised molecular cloud cores to form protostars. The calculations follow the formation and evolution of the first hydrostatic core, the collapse to form a stellar core, the launching of outflows from both the first hydrostatic core and stellar cores, and the breakout of the stellar outflow from the remnant of the first core. We investigate the roles of magnetic fields and thermal feedback on the outflow launching process, finding that both magnetic and thermal forces contribute to the launching of the stellar outflow. We also follow the stellar cores until they grow to masses of up to 20 Jupiter-masses, and determine their properties. We find that at this early stage, before fusion begins, the stellar cores have radii of $\approx 3$ R$_\odot$ with radial entropy profiles that increase outward (i.e. are convectively stable) and minimum entropies per baryon of $s/k_{\rm B} \approx 14$ in their interiors. The structure of the stellar cores is found to be insensitive to variations in the initial magnetic field strength. With reasonably strong initial magnetic fields, accretion on to the stellar cores occurs through inspiralling magnetised pseudo-discs with negligible radiative losses, as opposed to first cores which effectively radiate away the energy liberated in the accretion shocks at their surfaces. We find that magnetic field strengths of >10 kG can be implanted in stellar cores at birth.

Artificial resistivity is included in Smoothed Particle Magnetohydrodynamics simulations to capture shocks and discontinuities in the magnetic field. Here we present a new method for adapting the strength of the applied resistivity so that shocks are captured but the dissipation of the magnetic field away from shocks is minimised. Our scheme utilises the gradient of the magnetic field as a shock indicator, setting {\alpha}_B = h|gradB|/|B|, such that resistivity is switched on only where strong discontinuities are present. The advantage to this approach is that the resistivity parameter does not depend on the absolute field strength. The new switch is benchmarked on a series of shock tube tests demonstrating its ability to capture shocks correctly. It is compared against a previous switch proposed by Price & Monaghan (2005), showing that it leads to lower dissipation of the field, and in particular, that it succeeds at capturing shocks in the regime where the Alfv\'en speed is much less than the sound speed (i.e., when the magnetic field is very weak). It is also simpler. We also demonstrate that our recent constrained divergence cleaning algorithm has no difficulty with shock tube tests, in contrast to other implementations.

We present SPH formulations of Dedner et al's hyperbolic/parabolic divergence cleaning scheme for magnetic and velocity fields. Our implementation preserves the conservation properties of SPH which is important for stability. This is achieved by deriving an energy term for the Psi field, and imposing energy conservation on the cleaning subsystem of equations. This necessitates use of conjugate operators for divB and gradPsi in the numerical equations. For both the magnetic and velocity fields, the average divergence error in the system is reduced by an order of magnitude with our cleaning algorithm. Divergence errors in SPMHD are maintained to < 1%, even for realistic 3D applications with a corresponding gain in numerical stability. Density errors for an oscillating elliptic water drop using weakly compressible SPH are reduced by a factor of two.

We present a constrained formulation of Dedner et al's hyperbolic/parabolic divergence cleaning scheme for enforcing the \nabla\dot B = 0 constraint in Smoothed Particle Magnetohydrodynamics (SPMHD) simulations. The constraint we impose is that energy removed must either be conserved or dissipated, such that the scheme is guaranteed to decrease the overall magnetic energy. This is shown to require use of conjugate numerical operators for evaluating \nabla\dot B and \nabla{\psi} in the SPMHD cleaning equations. The resulting scheme is shown to be stable at density jumps and free boundaries, in contrast to an earlier implementation by Price & Monaghan (2005). Optimal values of the damping parameter are found to be {\sigma} = 0.2-0.3 in 2D and {\sigma} = 0.8-1.2 in 3D. With these parameters, our constrained Hamiltonian formulation is found to provide an effective means of enforcing the divergence constraint in SPMHD, typically maintaining average values of h |\nabla\dot B| / |B| to 0.1-1%, up to an order of magnitude better than artificial resistivity without the associated dissipation in the physical field. Furthermore, when applied to realistic, 3D simulations we find an improvement of up to two orders of magnitude in momentum conservation with a corresponding improvement in numerical stability at essentially zero additional computational expense.

We have performed Smoothed Particle Magnetohydrodynamics (SPMHD) simulations demonstrating the production of collimated jets during collapse of 1 solar mass molecular cloud cores to form the `first hydrostatic core' in low mass star formation. Recently a number of candidate first core objects have been observed, including L1448 IRS2E, L1451-mm and Per Bolo 58, although it is not yet clear that these are first hydrostatic cores. Recent observations of Per Bolo 58 in particular appear to show collimated, bipolar outflows which are inconsistent with previous theoretical expectations. We show that low mass first cores can indeed produce tightly collimated jets (opening angles <~ 10 degrees) with speeds of ~2-7 km/s, consistent with some of the observed candidates. We have also demonstrated, for the first time, that such phenomena can be successfully captured in SPMHD simulations.