\thechapter Introduction


March 7, 2018 \degreeDoctor of Philosophy \schoolSchool of Mathematical Sciences \facultyFaculty of Science \prevdegreesB.Sc.(Hons), M.Sc.


\prefaceCopyright Notice

Under the Copyright Act 1968, this thesis must be used only under the normal conditions of scholarly fair dealing. In particular no results or conclusions should be extracted from it, nor should it be copied or closely paraphrased in whole or in part without the written consent of the author. Proper written acknowledgement should be made for any assistance obtained from this thesis.

I certify that I have made all reasonable efforts to secure copyright permissions for third-party content included in this thesis and have not knowingly added copyright content to my work without the owner’s permission.


Numerical methods to improve the treatment of magnetic fields in smoothed field magnetohydrodynamics (SPMHD) are developed and tested. A mixed hyperbolic/parabolic scheme is developed which “cleans” divergence error from the magnetic field. The method introduces a scalar field which is coupled to the magnetic field. A conservative form for the hyperbolic equations is obtained by first defining the energy content of the new field, then using it in the discretised Lagrangian to obtain equations which manifestly conserve energy. This is shown to require conjugate first derivative operators in the SPMHD cleaning equations. Average divergence error is shown to be 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 in SPMHD. The method is adapted to work with a velocity field, demonstrating that it can reduce density variations in weakly compressible SPH simulations by a factor of 2.

A switch to reduce dissipation of the magnetic field from artificial resistivity is developed. Discontinuities in the magnetic field are located by monitoring jumps in the gradient of the magnetic field at the resolution scale relative to the magnitude of the magnetic field. This yields a simple yet robust method to reduce dissipation away from shocked regions. 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-Alfvénic regime.

These numerical methods are compared against grid-based MHD methods by comparison of the small-scale dynamo amplification of a magnetic field in driven, isothermal, supersonic turbulence. We use the SPMHD code, Phantom, and the grid-based code, Flash. We find that the growth rate of Flash is largely insensitive to the numerical resolution, whereas Phantom shows a resolution dependence that arises from the scaling of the numerical dissipation terms. The saturation level of the magnetic energy in both codes is about of the mean kinetic energy, increasing with higher magnetic Reynolds numbers. Phantom requires lower resolution to saturate at the same energy level as Flash. The time-averaged saturated magnetic spectra have a similar shape between the two methods, though Phantom contains twice as much energy on large scales. Both codes have PDFs of magnetic field strength that are log-normal, which become lopsided as the magnetic field saturates. We find encouraging agreement between grid- and particle methods for ideal MHD, concluding that SPMHD is able to reliably simulate the small-scale dynamo amplification of magnetic fields. We note that quantitative agreement on growth rates can only be achieved by including explicit, physical terms for viscosity and resistivity, because those are the terms that primarily control the growth rate and saturation level of the turbulent dynamo.

\prefaceDeclaration of Published Material

In accordance with Monash University Doctorate Regulation 17.2 Doctor of Philosophy and Research MasterÕs regulations the following declarations are made:

I hereby declare that this thesis contains no material which has been accepted for the award of any other degree or diploma at any university or equivalent institution and that, to the best of my knowledge and belief, this thesis contains no material previously published or written by another person, except where due reference is made in the text of the thesis.

This thesis includes 3 original papers published in or submitted to peer reviewed journals and 2 unpublished conference proceedings. The ideas, development and writing up of all the papers in the thesis were the principal responsibility of myself, the candidate, working within the School of Mathematical Sciences under the supervision of Dr. Daniel Price.

The inclusion of co-authors reflects the fact that the work came from active collaboration between researchers and acknowledges input into team-based research.

I have renumbered sections of submitted or published papers in order to generate a consistent presentation within the thesis.

The following chapters and/or sections have appeared in conference proceedings, have been published in peer reviewed journals, or have been submitted for publication:

  • Sections 46, 9, and Appendix \thechapter appear in
    Tricco, T. S. and Price, D. J.: 2012, ‘Constrained hyperbolic divergence cleaning for smoothed particle magnetohydrodynamics’. J. Comput. Phys. 231, 7214–7236.

  • Section 8 appears in
    Tricco, T. S. and Price, D. J.: 2012, ‘Hyperbolic divergence cleaning for SPH’. Proceedings of the 7th International SPHERIC Workshop.

  • Sections 10, 11, and 13 appear in
    Tricco, T. S. and Price, D. J.: 2013, ‘A switch to reduce resistivity in smoothed particle magnetohydrodynamics’. MNRAS 436, 2810–2817.

  • Section 12 appears in
    Tricco, T. S. and Price, D. J.: 2013, ‘A switch for artificial resistivity and other dissipation terms’, Proceedings of the 8th International SPHERIC Workshop.

  • Chapter \thechapter and Appendices \thechapter and \thechapter to appear in
    Tricco, T. S., Price, D. J., and Federrath, C.: 2014, ‘A comparison between grid and particle methods on small-scale dynamo amplification of magnetic fields in supersonic turbulence’, Submitted to ApJ.


I remember the moment I realised I should go to Australia. I enjoyed working with SPMHD and thought there was a lot of room for innovation in that area. At first I hadn’t considered doing a Ph.D. outside of Canada, but then the decision became obvious — I had been reading primarily the work of Daniel Price, and he seemed like a good guy: I should do my Ph.D. with him.

Moving to Melbourne in general has really been an positive influence on my life. I’ve thoroughly enjoyed my time at Monash, in Melbourne, and in Australia. I have had more opportunities for my Ph.D. than I could have hoped for, and I have made a lot of great friends and had great experiences. I’m sad to see this chapter of my life coming to a close.

I first and foremost sincerely wish to thank my supervisor, Daniel Price. He has been an incredible teacher, and there are many anecdotes I would like to share to demonstrate that, but then this would end up becoming very long. Suffice to say, he has really taught me a lot — directly and indirectly — on what takes to be a good researcher. The support, encouragement, and understanding that he has given me exceeded what I could have expected. The opportunities that he has given me have made a significant impact on the quality and enjoyment of my Ph.D. life. Many times I wished to nominate him for supervisor of the year, but was unable to because the regulations for these awards required he had to have supervised a number of students already. He definitely deserves it.

On a practical note, I would like to acknowledge the financial support of the following: The Australian government for supporting my candidature through an Australian Postgraduate Award and an Endeavour International Postgraduate Research Scholarship. The Astronomical Society of Australia for their travel grant supporting my attendance at Protostars and Planets VI in Heidelberg, Germany. Monash University, the Faculty of Science, the School of Mathematical Sciences, and the Monash Centre for Astrophysics for their respective contributions for travel funding throughout my candidature. Daniel Price for supplying a top-up to my stipend scholarship and funding conference travel through his ARC Discovery Grant DP1094585.

Computational resources have been utilised at the Multi-modal Australian ScienceS Imaging and Visualisation Environment (MASSIVE) through the National Computational Merit Allocation Scheme supported by the Australian Government (project NCIdz3), gStar through Astronomy Australia Ltd’s Astronomy Supercomputer Time Allocation Committee, and the Leibniz-Rechenzentrum (grant pr32lo).

Thank you to Joe Monaghan, Matthew Bate, Christoph Federrath, and Guillaume Laibe for their help throughout various stages of my Ph.D. I would also like to acknowledge my fellow Ph.D. students at Monash who have overlapped their time with mine, both ones finished and ones recently started. I would like to make particular mention of Tim Dolley, Nicolas Bonne, and David Palamara, the original trio who started at the same time as me. As well, Joelene Buntain and Hauke Worpel for being great officemates.

Thank you to my friends and loved ones back home: Kate Murphy; Michael Healey, Stephen Hinchey; John Hawkin, Hugh Newman, Marek Bromberek, and the rest of the MUN physics crowd; Geoff Barnes, Deanna Norman, and the rest of the ‘newfs’; Sharon Griffin, Donna Acorn, and all of the Griffins; Jack & Doreen Tricco, Doug Tricco, Sheila Wadland, and all of the Triccos. I’ve missed you guys a lot, and you’ve all helped in your own ways whether you’ve known it or not. Thank you to my brother, Jon, and a sincere thank you to my mother, Brenda. Reaching this point is in no small part due to her continual support. To my father, Paul, wish you were still around so you could see how far I’ve got.

List of Figures:

Chapter \thechapter Introduction

… magnetic fields may be included without difficulty …

Gingold and Monaghan (1977)

Magnetic fields are ubiquitous throughout the Universe. It is believed that even if the Universe began unmagnetised, battery effects would lead to an initial magnetisation of baryonic matter (e.g., the Biermann battery, Biermann 1950, see also the review by Widrow et al. 2012). Since magnetic monopoles do not exist in nature, there are no ‘sinks’ of magnetic field and it is difficult to destroy them. Therefore, once an initial magnetisation is present, dynamo processes can lead to ever stronger magnetic fields.

Nearly all current theoretical problems in astrophysics involve magnetic fields to some degree. Neutron stars have some of the strongest magnetic fields in the Universe. The magnetic fields of galaxies are thought to be dynamically relevant for their evolution, and are responsible for determining the propagation of cosmic rays. The magnetic field of the Sun is responsible for sunspots and solar flares.

Magnetic fields also play an important role in all stages of the star formation process. Stars form in cold (K) clouds of molecular gas (primarily ), which contain between of material. Supersonic turbulence in these clouds plays a key role in regulating star formation (see review by McKee and Ostriker, 2007). As the supersonic shock waves collide in the cloud, they create dense filaments which act as the nucleation sites along which stars can form. These provide the dense cores that begin the star formation process (Larson, 1981). The extra pressure from magnetic fields help guard against gravitational collapse, and numerical studies has shown that this can reduce star formation rates (e.g., Nakamura and Li, 2008; Price and Bate, 2008, 2009; Padoan and Nordlund, 2011; Federrath and Klessen, 2012).

On the scale of individual protostars, magnetic fields are responsible for driving jets and outflows — a signature of star formation. As a molecular cloud core collapses under its gravitational weight, conservation of angular momentum leads to an increase in angular velocity, winding up magnetic field lines. There are two ways in which the magnetic field may drive an outflow. One occurs when the tension in the field lines becomes too strong, driving material outwards as the magnetic field pops out of the plane of the disc (the ‘magnetic tower’, Lynden-Bell, 1996, 2003). The other is when material is centrifugally accelerated along poloidal magnetic field lines, essentially being ‘sling shotted’ away from the protostar (Blandford and Payne, 1982). Outflows are important sources of removing angular momentum from the star-disc system and in reducing the efficiency of gas conversion into stars.

Magnetic fields also play an important role in the accretion discs around young stars. Magnetised, differentially rotating flows are well known to be susceptible to the magneto-rotational instability (Balbus and Hawley, 1991). Consider two pieces of material on nearby orbits that are joined by a magnetic field line. As they drift apart, the tension in the magnetic field line resists the motion. This pulls the two pieces towards each other, causing the material in the inner orbit to slow down, and the material in the outer orbit to speed up. However, this only causes the inner material to drop to a lower orbit, and the outer material to drift outwards, exacerbating the problem. By this process, angular momentum is transported outwards through the disc. This instability leads to turbulence, and is thought to play a key role in driving accretion onto young stars.

Observations of magnetic fields may be obtained directly through Zeeman splitting measurements of spectral lines, or indirectly by the linear polarisation of thermal emission from dust grains. However, Zeeman measurements only yield information about the magnetic field along the line of sight, and the polarisation of dust grains only about its orientation in the plane of the sky. Therefore, the full information about the magnetic field is difficult to obtain. Furthermore, performing these observations may require a considerable amount of time, for example, Troland and Crutcher (2008)’s survey of Zeeman measurements of magnetic fields in molecular cloud cores involved 500 hours of observing time.

An important approach to test astrophysical theories is through the use of numerical simulation, and it is crucial that these numerical experiments reflect reality as closely as possible in order to yield meaningful results. This is accomplished through careful design and calibration of numerical methods.

This thesis is focused on improving the treatment of magnetic fields in smoothed particle magnetohydrodynamics (SPMHD), a Lagrangian particle based numerical method built on smoothed particle hydrodynamics (SPH). The general picture of SPH is to solve the equations of hydrodynamics by discretising a fluid into a collection of particles that mimic fluid behaviour. SPH has many advantages for astrophysics. One, the resolution is tied to the mass. Regions of higher mass have more particles, thus more resolution, which is advantageous as the densest areas are typically the most interesting (e.g., stars forming in a molecular cloud). Two, it is trivial to incorporate gravitational N-body methods since SPH is a particle based scheme. Three, advection is done perfectly, that is, without any dissipation, since it is a Lagrangian method. Four, it can easily handle complex geometries. Five, the Courant timestep does not depend upon the fluid velocity, thus allowing larger timesteps. And six, perhaps its strongest attribute, it has exact simultaneous conservation of mass, momentum, angular momentum, energy, and entropy to the precision of the time-stepping algorithm. This makes SPH significantly robust and stable since it reflects the conservation properties of nature.

SPH is widely used in astrophysics for the preceding reasons. For example, the cosmological code Gadget 2 (Springel, 2005) has over 1900 citations at the time of writing. The impetus to include other physics in SPH, such as magnetic fields, is clear. The foundation of SPMHD has been laid substantially through the Ph.D. research of Daniel Price (see Price and Monaghan 2004a, b, 2005 and also the recent review by Price 2012), building on the earlier work of Phillips and Monaghan (1985) and Morris (1996). This thesis follows as its spiritual successor, shoring up the remaining deficiencies to build a method that is able to accurately simulate a wide range of astrophysical problems.

The thesis is structured as follows: In Chapter \thechapter, the current state of SPMHD is reviewed. First, the continuum equations of ideal MHD are derived, which is more than mere exercise as this will elucidate some of the numerical issues to be discussed. The numerical method is built up step-by-step. We provide an overview of how to estimate the density for the set of SPH particles, how to adaptively set the resolution based on density, how to perform basic interpolation of quantities, and how to obtain the discretised form of the induction, energy, and momentum equations. The numerical instability present in the equations of motion will be discussed, along with strategies for its treatment. Methods for the capturing of shocks and discontinuities are outlined.

In Chapter \thechapter, the constrained hyperbolic divergence cleaning method to uphold the divergence-free constraint of the magnetic field in SPMHD is developed and tested. Past approaches to treat the divergence-free constraint in SPMHD are summarised first. In Chapter \thechapter, a method to reduce numerical dissipation of the magnetic field is presented and tested. Chapter \thechapter presents a comparison of the SPMHD methods developed against grid-based methods on the simulation of small-scale dynamo amplification of a magnetic field. The thesis is summarised in Chapter \thechapter.

Chapter \thechapter Smoothed particle magnetohydrodynamics

Smoothed particle magnetohydrodynamics (SPMHD) is a numerical method for solving the equations of magnetohydrodynamics (MHD) based on the smoothed particle hydrodynamics (SPH) method (Gingold and Monaghan, 1977; Lucy, 1977). The basic premise is to discretise the fluid by mass into a set of particles. To recover continuum behaviour from the collection of point particles, a weighting kernel is used to smooth their quantities over a local volume. Fluid properties can be reconstructed at any point in space by summing the weighted contributions of nearby particles.

The first attempts to include magnetic fields in SPH were performed by Gingold and Monaghan (1977) who considered magnetic polytropes, though in a form which did not conserve momentum or angular momentum. The basic SPMHD method has its roots in the work by Phillips and Monaghan (1985), who formulated equations of motion that conserve momentum, and applied the method to three-dimensional simulations of gravitationally collapsing gas clouds (Phillips, 1986a, b). The modern SPMHD method was developed by Price and Monaghan (2004a, b, 2005), who constructed fully conservative equations incorporating varying resolution, formulated magnetic shock capturing terms, and investigated approaches to treat the divergence-free constraint on the magnetic field. Since then, SPMHD has been applied to studies of protostar formation (Price and Bate, 2007; Bürzle et al., 2011a, b; Price et al., 2012; Bate et al., 2014b), star cluster formation (Price and Bate, 2008, 2009), neutron star mergers (Price and Rosswog, 2006), and magnetic fields in galaxies and galaxy clusters (Price and Bate, 2008; Donnert et al., 2009; Kotarba et al., 2009, 2010, 2011; Bonafede et al., 2011; Beck et al., 2012, 2013).

Technical difficulties in SPMHD arise from the divergence-free constraint on the magnetic field – an issue faced by any numerical MHD method. Magnetic monopoles are introduced if this constraint is not upheld, which is not only physically inaccurate, but leads to spurious monopole accelerations. This causes numerical instability when it exceeds the isotropic pressure. The details of the issues surrounding in SPMHD will be discussed as the method is presented throughout this chapter and in Chapter \thechapter.

We begin by deriving the continuum equations of ideal MHD along with the MHD wave solutions. This process is instructive and leads to further understanding of some of the finer points of the numerical scheme. The SPMHD discretised version of the MHD equations will then be constructed. It begins with a method to estimate density in SPH, and a review of basic SPH interpolation theory. With this, the discretised induction equation used to evolve the magnetic field can be obtained. Using the density estimate and discretised induction equation, the conservative equations of motion are built through a Lagrangian approach. The instability present in these equations will be discussed, along with approaches for removing it, making particular note of how it is related to the divergence-free constraint of the magnetic field. Dissipation terms for capturing shocks are presented, followed by methods to reduce dissipation.

1 Continuum magnetohydrodynamics

Ideal MHD is the merger of fluid dynamics with electromagnetic theory. Several useful textbooks for magnetohydrodynamic theory are Choudhuri (1998); Griffiths (1999); Batchelor (2000); Bellan (2006). The relevant equations are given by Euler fluid flow, describing the motion of an inviscid fluid,


Maxwell’s equations of electromagnetism,


and the Lorentz force law,


Here, is the fluid velocity, is the density, is the thermal pressure, is the electric field, is the magnetic field, is the charge density, is the current density, is the permeability of free space, and is the permittivity.

One key assumption is made: The fluid is highly ionised. This means that while the fluid can carry a magnetic field, on macroscopic scales (relevant for astrophysical systems), the positive and negative charges will average out and the fluid will be electrically neutral. Furthermore, since there is a significant number of free electrons, the fluid can be treated as an ideal conductor. Both of these conditions imply that the stationary electric field inside the fluid can be treated as negligible.

1.1 Momentum equation

Forces from the magnetic field are due to the Lorentz force law (Equation 7). Assuming that and does not vary with time, we can use Equation 5 to define the current density in terms of the magnetic field, such that the Lorentz force becomes


This may be rewritten as


from which it becomes clear that the magnetic field exerts two forces on the fluid. One is an isotropic magnetic pressure, which pushes fluid down gradients of magnetic field strength. The second is an attractive force directed along magnetic field lines, which functions like a tension in the magnetic field lines.

The total force on the fluid is the combination of pressure and magnetic forces. The momentum equation is therefore the addition of Equation 2 and 9, yielding


Here we have introduced the material derivative, , which follows the frame of reference of a parcel of fluid along its streamline. As SPH is a particle based method, it is natural to write equations using the material derivative.

The momentum equation can be written in terms of a stress tensor. Assuming that the magnetic field is divergence-free, the stress tensor can be defined as


which leads to


Expanding this, the momentum equation becomes


This is similar to Equation 10 except in the magnetic tension term. It contains an extra tensional force, , which appears due to the assumption that . The conservative form of the SPMHD momentum equation is obtained by using the stress tensor, though since it may be unsafe to assume the magnetic field is divergence-free when solving the equations numerically, this extra force term requires careful consideration. This is discussed in Sections 2.7 and 2.8.

1.2 Induction equation

The current density may be defined using Ohm’s law,


which expresses the current density, , in terms of the electric field in the co-moving frame of an observer and the electrical conductivity, , of the material, which is treated as constant. In a fixed frame of reference, the electric field is given by


giving Ohm’s law as


which is the combination of current induced by an electric field and by moving through a magnetic field. This permits the electric field in to be expressed as


Taking the curl of Equation 17, the electric field may be replaced using Equation 3 to obtain an evolution equation for the magnetic field as


In the limit of ideal MHD (infinite conductivity, ), the term involving the current density drops out. For non-ideal MHD, may be replaced using Equation 5 (neglecting the displacement current, ), obtaining


where is the magnetic resistivity.

In ideal MHD, the conductivity of the fluid is taken to be infinite, therefore . Expanding the first term of Equation 19, we can write


or by using and the material derivative,


The first term affects the magnetic field through shearing motion, while the second will increase the magnetic field when undergoing compression.

1.3 Summary of ideal MHD equations

The concise set of ideal MHD equations to be solved are


1.4 Wave solutions

The ideal MHD wave equations permit three wave modes, not just sound waves as found in simple fluids (i.e., Euler fluids). Understanding the ideal MHD wave solutions will be useful when introducing shock capturing schemes to our numerical method.

The ideal MHD wave modes can be obtained as follows. Assume a uniform density fluid at rest with a constant magnetic field. The equation of state is taken to be isothermal, , where is the speed of sound. Small perturbations are introduced to the density, velocity, and magnetic fields such that


where and are the background density and magnetic fields, with , , and perturbations to each field. The perturbations are taken to be sufficiently small so as to not disturb the equilibrium values of the fluid, thus the background fields remain constant in time. Inserting Equations 2628 into the ideal MHD equations, the set of linearised equations are then


Second order effects are assumed negligible so those terms involving multiple perturbations are discarded.

We assume that the perturbations have wave-like solutions of the form , where is the wave vector. Equations 29 and 30 become


where in deriving Equation 33, we have made use of Equation 32 to substitute . Taking the time derivative of Equation 33, and using Equation 31, we obtain


where we have defined , known as the Alfvén speed.

Taking the magnetic field to be in the -direction and vector in the - plane, that is and , Equation 34 yields the following set of equations,


The first component is


Since this is directed along , orthogonal to both the direction of wave propagation and magnetic field, this is a transverse oscillation known as an Alfvén wave. These have phase velocity , directed along the orientation of the magnetic field. Alfvén waves can be understood as occurring from the tension present in magnetic field lines, operating in a similar fashion to vibrating strings in a string instrument.

Waves in the - plane are found from the determinant,


Using the quadratic formula, we can solve for , finding phase velocity


The wave velocity is in the same plane as the wave vector, therefore these are longitudinal waves. If no magnetic field is present, that is , these reduce to ordinary sound waves. These two wave types occur from the combination of magnetic and thermal pressure. One wave mode has its speed boosted by the addition of magnetic pressure (called fast MHD wave). The second wave mode is preferentially guided by the magnetic field, such that propagation across magnetic field lines is halted (slow MHD waves).

2 Discretised magnetohydrodynamics

The beauty of SPH is in its simplicity and intuitiveness. The basic method can be derived from first principles such that the discretised equations which are solved are the physical equations governing the system of discrete particles. Therefore, SPH inherently has the conservation properties of real physics, giving the method numerical stability and robustness.

In this section, an overview of SPMHD is presented. Focus will be on how the MHD equations (2225) are solved numerically, highlighting how the discretised equations are obtained and the numerical challenges specific to SPMHD. For a complete background on interpolation theory, properties of smoothing kernels, stability analysis, and other deeper technical issues, the reader is referred to the PhD thesis of Morris (1996) and the reviews by Monaghan (2005) and Price (2012).

2.1 Estimating the density

The first step is to obtain an estimate of the density. This is accomplished by taking a weighted summation of the mass of neighbouring particles within a characteristic radius , known as the smoothing length. The density in SPH is estimated as


where is the weighting function known as the smoothing kernel. We assume that each particle is allowed its own smoothing length and that it is spatially varying.

The density sum of Equation 39 can be used to calculate the density of a particle whenever required. It is unnecessary under general circumstances to evolve the density of a particle using its time derivative. However, for constructing the SPMHD equations, the discretised version of the continuity equation is useful. It can be obtained by taking the time derivative of Equation 39, yielding


where we have made use of the chain rule and introduced the shorthand notation . Using and the anti-symmetry of the kernel gradient, that is , this can be simplified to yield


where is the gradient taken respect to the coordinates of particle and


The term is important to correctly account for spatially varying smoothing lengths (see Springel and Hernquist, 2002; Monaghan, 2002). In general, they should be used when computing any derivative estimate. In particular, inclusion of these terms into the SPMHD momentum, induction, and energy equations has been shown to improve the representation of wave propagation and shocks (Price and Monaghan, 2004b). Obtaining may be done through Equation 43 below, as given by Equation 44.

2.2 Setting the smoothing length

The smoothing length is individually set per particle by mutually solving


with the density summation (Equation 39). The derivative is given by


and Equation 43 leads to an expression for the density as


Here, , , is the dimension of the system and is a dimensionless quantity specifying the ratio of smoothing length to particle spacing. For the spline family of kernels (Schönberg, 1946; Monaghan, 1985; Monaghan and Lattanzio, 1985), this is typically chosen to be . Since the density itself is a function of smoothing length, this requires iteration until both quantities converge. This is an expensive process, as the density summation needs to be re-evaluated for each iteration, which may further necessitate performing a neighbour search.

A root finding technique can be used to find the smoothing length and density for each particle. The function to find the root of is


That is, the root, , occurs when the expected density (from Equation 45) agrees with the density as calculated through summation (Equation 39).

The root can be found with the Newton-Raphson technique (e.g. Price and Monaghan, 2004b, 2007). The smoothing length may be iterated according to


where is the first derivative of ,


For convenience of implementation, this may be rewritten using as


where is the expression in Equation 45. By using the tangent of to iterate towards the root, the method has second order convergence and as such is an efficient means to solve for and . However, it may fail if the tangent of is nearly parallel, leading to the next iteration to significantly overshoot the root and diverge. This risk can be curbed by including a check to restrict modification of between iterations by no more than, say, . For a typical SPMHD calculation, the risk of the method diverging is minimal (and usually indicative of a more serious problem elsewhere).

An approach guaranteed to converge, though slower with only linear convergence, is to use a bisection method. The method is simple. If the converged value of is known to lie within a specified interval, then the interval can be halved until it is found. Which half of the interval to reject can be determined through , where if , then should be decreased from its current value, otherwise should be increased on the next iteration.

Convergence can be determined by monitoring the relative difference in (or ) between iterations. This is determined according to


where . Note that , the value of the smoothing length before the first iteration, is used so that the denominator remains fixed and convergence occurs only when agrees with .

A simple approach to reduce the number of overall iterations is to time integrate the smoothing length, predicting a value close to the root before beginning the root finding technique. Taking the time derivative of Equation 43, and using the continuity equation (Equation 22), we obtain


Given that is usually calculated in an SPH code, this adds almost no additional computational cost, yet can significantly increase the overall efficiency of the code.

At the beginning of each timestep, the step-by-step approach to setting the smoothing length per particle is: {enumerate*} Compute the density using summation. Iterate to using Equation 47 (or other root-solving technique). Do and agree within the specified tolerance (Equation 50)? {enumerate*} Yes, accept and proceed to step 4. No, begin again from step 1 Predict for the next timestep using Equation 51.

2.3 Interpolation basics

In order to define the discretised induction equation, it is necessary to understand the basics of SPH interpolation theory. This is not a comprehensive review, but will introduce the basics necessary to formulate the SPMHD equations. Throughout this section, we assume that the smoothing length is uniform and constant so that the presentation may be clearer. Inserting terms to account for variable smoothing lengths may be appropriately inserted where gradients of the smoothing kernel are taken.

In the continuum limit, the value of a quantity can be obtained by using a delta function to pluck that value at a specified location, that is


In SPH, the smoothing kernel plays the role of the delta function. It has property such that


The kernel is assumed to be spherically symmetric, and normalised such that . In a discretised system, the integral is replaced by a summation over elements. The quantity can be obtained through


where acts like the volume element of the integral. This reduces to the density summation if , and is the traditional way to introduce SPH. In this thesis however, we use the mass weighted summation


where acts like the normalisation on the summation. If is a constant, this returns exactly. Though for clarity of presentation, we continue using Equation 54.

The derivative of Equation 54 is


However, this yields a poor estimate of the gradient. For example, constant functions will yield a non-zero result. Higher accuracy gradient estimates may be obtained by taking a Taylor series expansion to obtain error terms, then subtracting those errors from the gradient (Morris, 1996; Price, 2004). The Taylor series expansion of in Equation 56 about is


Thus, the gradient estimate can be made first order by subtracting the first error term in the Taylor series, yielding


which is exact for constant functions. This is the most common form for calculating gradients in SPH. Obtaining divergence and curl estimates for vector quantities may be obtained through a similar procedure as the preceding, with the divergence and curl of the magnetic field presented in Section 2.4.

A second order estimate can be obtained by subtracting the second error term of the Taylor series, yielding




is a matrix that acts as a correction to the kernel gradient. This approach has been used by Bonet and Lok (1999). This second order derivative estimate adds more computational expensive since it requires a matrix inversion and storage of the matrix elements.

While second derivatives may be estimated by taking the derivative of Equation 58, this is quite sensitive to particle disorder and leads to a poor estimate. A less noisy estimate is to take two first derivatives, applying Equation 58 (or Equation 59) twice in succession. This is more expensive as it requires two loops over particle neighbours. An alternative approach is that of Brookshaw (1985), whereby a second derivative estimate is obtained from


This may also be written as


making use of the definition , where is the scalar portion of the kernel derivative. We note that there is inconsistent usage of this definition in the literature. Notably the SPH review by Price (2012) uses the aforementioned definition, whereas the SPH review by Monaghan (2005) instead uses , differing by a factor . A reader should be careful of this difference in the literature. In this thesis, we adopt usage consistent with Price (2012) ().

This second derivative estimate may be obtained through Taylor series expansion of about in the integral approximation,


This functionally approximates a second derivative by dividing the first derivative of the smoothing kernel by the particle spacing, .

2.4 and

The divergence and curl of the magnetic field may be obtained in a similar manner to the first derivative estimates of scalar quantities. The divergence and curl equivalents of Equation 56 are given by




It is possible to obtain first order accurate estimates of the first derivatives through other means than the Taylor series expansion presented in the preceding section. For the divergence operator, consider the identity


Inserting the simple first derivative operators from Equations 56, 64 and 65 will yield first order accurate estimates. In this thesis, we use the mass weighted summations, with the first order accurate divergence and curl of the magnetic field given by (including variable smoothing length terms)




Other identities lead to other first derivative estimates. The identity


may be used to obtain an entirely different form for the first derivative. The divergence and curl with this operator (including variable smoothing length terms) are




The error for these estimates is large (), and yield non-zero results for constant functions.

We refer to the first derivative operators as the ‘difference’ measure (Equation 67 and 68) and the ‘symmetric’ measure (Equation 70 and 71). It is noteworthy that the symmetric measure of is what will appear in the equations of motion, and the implications of this derivative estimate will be discussed in Section 2.8.

2.5 Energy equation

The equations of motion need to be coupled to an equation of state to determine the thermal pressure. If the pressure is a function of internal energy, , a suitable equation must be used to evolve in time. Consider the thermodynamic relation,


where is the temperature, is the entropy, is the volume, and quantities are expressed per unit mass. Since SPMHD is inherently dissipationless, may be taken to be zero. (Alternatively, if the pressure is a function of entropy, the entropy per particle may be passively advected and increased only from added sources of dissipation.) Converting to be per unit mass, the time derivative of is


Using the SPH continuity equation (Equation 41), the discretised internal energy equation is thus


2.6 Induction equation

Using basic interpolation theory, the induction equation (Equation 24) may be discretised as


Alternatively, the quantity could be evolved. Rewriting the induction equation as


the SPMHD form is


Both approaches are utilised throughout this thesis, depending on the code used. Neither approach confers any significant advantage over the other (Price, 2012).

2.7 Conservative equations of motion

The equations of motion for SPMHD will be obtained by using the Lagrangian for the discretised system (Price, 2004; Price and Monaghan, 2004b). This will yield the equations of motion that physically govern the system, providing a method that has exact conservation of momentum, energy, and entropy. Consider the SPMHD Lagrangian,


The action integral, , is stationary. Therefore, small perturbations must not change the solution, that is


If small deviations are introduced into the Lagrangian about , then


Using the density summation (Equation 39) and induction equation (Equation 75) as constraints, the variations and can be written in terms of according to


Inserting these into Equation 80 and using the thermodynamic relation (Equation 73), we obtain


The perturbations in and can be removed by multiplying the latter terms by , introducing delta functions into the equation and yielding


Simplifying out the delta functions and using the anti-symmetry of the kernel gradient (), we obtain


Inserting Equation 85 into 79 and integrating the velocity term by parts, the equations of motion are found to be