We present an algorithm for simulating the equations of ideal magnetohydrodynamics and other systems of differential equations on an unstructured set of points represented by sample particles. Local, third-order, least-squares, polynomial interpolations (Moving Least Squares interpolations) are calculated from the field values of neighboring particles to obtain field values and spatial derivatives at the particle position. Field values and particle positions are advanced in time with a second order predictor-corrector scheme. The particles move with the fluid, so the time step is not limited by the Eulerian Courant-Friedrichs-Lewy condition. Full spatial adaptivity is implemented to ensure the particles fill the computational volume, which gives the algorithm substantial flexibility and power. A target resolution is specified for each point in space, with particles being added and deleted as needed to meet this target. Particle addition and deletion is based on a local void and clump detection algorithm. Dynamic artificial viscosity fields provide stability to the integration. The resulting algorithm provides a robust solution for modeling flows that require Lagrangian or adaptive discretizations to resolve. This paper derives and documents the Phurbas algorithm as implemented in Phurbas version 1.1. A following paper presents the implementation and test problem results.

We present an algorithm for simulating the equations of ideal magnetohydrodynamics and other systems of differential equations on an unstructured set of points represented by sample particles. The particles move with the fluid, so the time step is not limited by the Eulerian Courant-Friedrichs-Lewy condition. Full spatial adaptivity is required to ensure the particles fill the computational volume, and gives the algorithm substantial flexibility and power. A target resolution is specified for each point in space, with particles being added and deleted as needed to meet this target. We have parallelized the code by adapting the framework provided by GADGET-2. A set of standard test problems, including 1e-6 amplitude linear MHD waves, magnetized shock tubes, and Kelvin-Helmholtz instabilities is presented. Finally we demonstrate good agreement with analytic predictions of linear growth rates for magnetorotational instability in a cylindrical geometry. This paper documents the Phurbas algorithm as implemented in Phurbas version 1.1.

Finite-difference simulations of fluid dynamics and magnetohydrodynamics generally require an explicit diffusion operator, either to maintain stability by attenuating grid-scale structure, or to implement physical diffusivities such as viscosity or resistivity. If the goal is stability only, the diffusion must act at the grid scale, but should affect structure at larger scales as little as possible. For physical diffusivities the diffusion scale depends on the problem, and diffusion may act at larger scales as well. Diffusivity undesirably limits the computational timestep in both cases. We construct tuned finite-difference diffusion operators that minimally limit the timestep while acting as desired near the diffusion scale. Such operators reach peak values at the diffusion scale rather than at the grid scale, but behave as standard operators at larger scales. We focus on the specific applications of hyperdiffusivity for numerical stabilization, and high Schmidt and high Prandtl number simulations where the diffusion scale greatly exceeds the grid scale.

Numerical simulations including magnetic fields have become important in many fields of astrophysics. Evolution of magnetic fields by the constrained transport algorithm preserves magnetic divergence to machine precision, and thus represents one preferred method for the inclusion of magnetic fields in simulations. We show that constrained transport can be implemented with volume-centered fields and hyperresistivity on a high-order finite difference stencil. Additionally, the finite-difference coefficients can be tuned to enhance high-wavenumber resolution. Similar techniques can be used for the interpolations required for dealiasing corrections at high wavenumber. Together, these measures yield an algorithm with a wavenumber resolution that approaches the theoretical maximum achieved by spectral algorithms. Because this algorithm uses finite differences instead of fast Fourier transforms, it runs faster and isn't restricted to periodic boundary conditions. Also, since the finite differences are spatially local, this algorithm is easily scalable to thousands of processors. We demonstrate that, for low-Mach-number turbulence, the results agree well with a high-order, non-constrained-transport scheme with Poisson divergence cleaning.

We study numerically the dependence of the critical magnetic Reynolds number Rmc for the turbulent small-scale dynamo on the hydrodynamic Reynolds number Re. The turbulence is statistically homogeneous, isotropic, and mirror--symmetric. We are interested in the regime of low magnetic Prandtl number Pm=Rm/Re<1, which is relevant for stellar convective zones, protostellar disks, and laboratory liquid-metal experiments. The two asymptotic possibilities are Rmc->const as Re->infinity (a small-scale dynamo exists at low Pm) or Rmc/Re=Pmc->const as Re->infinity (no small-scale dynamo exists at low Pm). Results obtained in two independent sets of simulations of MHD turbulence using grid and spectral codes are brought together and found to be in quantitative agreement. We find that at currently accessible resolutions, Rmc grows with Re with no sign of approaching a constant limit. We reach the maximum values of Rmc~500 for Re~3000. By comparing simulations with Laplacian viscosity, fourth-, sixth-, and eighth-order hyperviscosity and Smagorinsky large-eddy viscosity, we find that Rmc is not sensitive to the particular form of the viscous cutoff. This work represents a significant extension of the studies previously published by Schekochihin et al. 2004, PRL 92, 054502 and Haugen et al. 2004, PRE, 70, 016308 and the first detailed scan of the numerically accessible part of the stability curve Rmc(Re).

We report an extensive numerical study of the small-scale turbulent dynamo at large magnetic Prandtl numbers Pm. A Pm scan is given for the model case of low-Reynolds-number turbulence. We concentrate on three topics: magnetic-energy spectra and saturation levels, the structure of the field lines, and the field-strength distribution. The main results are (1) the folded structure (direction reversals at the resistive scale, field lines curved at the scale of the flow) persists from the kinematic to the nonlinear regime; (2) the field distribution is self-similar and appears to be lognormal during the kinematic regime and exponential in the saturated state; and (3) the bulk of the magnetic energy is at the resistive scale in the kinematic regime and remains there after saturation, although the spectrum becomes much shallower. We propose an analytical model of saturation based on the idea of partial two-dimensionalization of the velocity gradients with respect to the local direction of the magnetic folds. The model-predicted spectra are in excellent agreement with numerical results. Comparisons with large-Re, moderate-Pm runs are carried out to confirm the relevance of these results. New features at large Re are elongation of the folds in the nonlinear regime from the viscous scale to the box scale and the presence of an intermediate nonlinear stage of slower-than-exponential magnetic-energy growth accompanied by an increase of the resistive scale and partial suppression of the kinetic-energy spectrum in the inertial range. Numerical results for the saturated state do not support scale-by-scale equipartition between magnetic and kinetic energies, with a definite excess of magnetic energy at small scales. A physical picture of the saturated state is proposed.

We consider the problem of incompressible, forced, nonhelical, homogeneous, isotropic MHD turbulence with no mean magnetic field. This problem is essentially different from the case with externally imposed uniform mean field. There is no scale-by-scale equipartition between magnetic and kinetic energies as would be the case for the Alfven-wave turbulence. The isotropic MHD turbulence is the end state of the turbulent dynamo which generates folded fields with small-scale direction reversals. We propose that the statistics seen in numerical simulations of isotropic MHD turbulence could be explained as a superposition of these folded fields and Alfven-like waves that propagate along the folds.

The amplification of magnetic fields in a highly conducting fluid is studied numerically. During growth, the magnetic field is spatially intermittent: it does not uniformly fill the volume, but is concentrated in long thin folded structures. Contrary to a commonly held view, intermittency of the folded field does not increase indefinitely throughout the growth stage if diffusion is present. Instead, as we show, the probability-density function (PDF) of the field strength becomes self-similar. The normalized moments increase with magnetic Prandtl number in a powerlike fashion. We argue that the self-similarity is to be expected with a finite flow scale and system size. In the nonlinear saturated state, intermittency is reduced and the PDF is exponential. Parallels are noted with self-similar behavior recently observed for passive-scalar mixing and for map dynamos.

We consider the problem of incompressible, forced, nonhelical, homogeneous and isotropic MHD turbulence with no mean magnetic field and large magnetic Prandtl number. This type of MHD turbulence is the end state of the turbulent dynamo, which generates folded fields with small-scale direction reversals. We propose a model in which saturation is achieved as a result of the velocity statistics becoming anisotropic with respect to the local direction of the magnetic folds. The model combines the effects of weakened stretching and quasi-two-dimensional mixing and produces magnetic-energy spectra in remarkable agreement with numerical results at least in the case of a one-scale flow. We conjecture that the statistics seen in numerical simulations could be explained as a superposition of these folded fields and Alfv\'en-like waves that propagate along the folds.

We introduce a fluid dynamics algorithm that performs with nearly spectral accuracy, but uses finite-differences instead of FFTs to compute gradients and thus executes 10 times faster. The finite differencing is not based on a high-order polynomial fit. The polynomial scheme has supurb accuracy for low-wavenumber gradients but fails at high wavenumbers. We instead use a scheme tuned to enhance high-wavenumber accuracy at the expense of low wavenumbers, although the loss of low-wavenumber accuracy is negligibly slight. A tuned gradient is capable of capturing all wavenumbers up to 80 percent of the Nyquist limit with an error of no worse than 1 percent. The fact that gradients are based on finite differences enables diverse geometries to be considered and eliminates the parallel communications bottleneck.

We report a series of numerical simulations showing that the critical magnetic Reynolds number Rm_c for the nonhelical small-scale dynamo depends on the Reynolds number Re. Namely, the dynamo is shut down if the magnetic Prandtl number Pr=Rm/Re is less than some critical value Pr_c<1 even for Rm for which dynamo exists at Pr>=1. We argue that, in the limit of Re->infinity, a finite Pr_c may exist. The second possibility is that Pr_c->0 as Re->infinity, while Rm_c tends to a very large constant value inaccessible at current resolutions. If there is a finite Pr_c, the dynamo is sustainable only if magnetic fields can exist at scales smaller than the flow scale, i.e., it is always effectively a large-Pr dynamo. If there is a finite Rm_c, our results provide a lower bound: Rm_c<220 for Pr<=1/8. This is larger than Rm in many planets and in all liquid-metal experiments.

Fast particles diffusing along magnetic field lines in a turbulent plasma can diffuse through and then return to the same eddy many times before the eddy is randomized in the turbulent flow. This leads to an enhancement of particle acceleration by large-scale compressible turbulence relative to previous estimates in which isotropic particle diffusion is assumed.

We investigate field-line separation in strong MHD turbulence using direct numerical simulations. We find that in the static-magnetic-field approximation the thermal conductivity in galaxy clusters is reduced by a factor of about 50 relative to the Spitzer thermal conductivity of a non-magnetized plasma. This value is too small for heat conduction to balance radiative cooling in clusters.

We investigate field-line separation in strong MHD turbulence analytically and with direct numerical simulations. We find that in the static-magnetic-field approximation the thermal conductivity in galaxy clusters is reduced by a factor of about 5-10 relative to the Spitzer thermal conductivity of a non-magnetized plasma. We also estimate how the thermal conductivity would be affected by efficient turbulent resistivity.

The growth and saturation of magnetic field in conducting turbulent media with large magnetic Prandtl numbers are investigated. This regime is very common in low-density hot astrophysical plasmas. During the early (kinematic) stage, weak magnetic fluctuations grow exponentially and concentrate at the resistive scale, which lies far below the hydrodynamic viscous scale. The evolution becomes nonlinear when the magnetic energy is comparable to the kinetic energy of the viscous-scale eddies. A physical picture of the ensuing nonlinear evolution of the MHD dynamo is proposed. Phenomenological considerations are supplemented with a simple Fokker--Planck model of the nonlinear evolution of the magnetic-energy spectrum. It is found that, while the shift of the bulk of the magnetic energy from the subviscous scales to the velocity scales may be possible, it occurs very slowly -- at the resistive, rather than dynamical, time scale (for galaxies, this means that generation of large-scale magnetic fields cannot be explained by this mechanism). The role of Alfvenic motions and the implications for the fully developed isotropic MHD turbulence are discussed.

We study the intermittency and field-line structure of the MHD turbulence in plasmas with very large magnetic Prandtl numbers. In this regime, which is realized in the interstellar medium, some accretion disks, protogalaxies, galaxy-cluster gas, early Universe, etc., magnetic fluctuations can be excited at scales below the viscous cutoff. The salient feature of the resulting small-scale magnetic turbulence is the folded structure of the fields. It is characterized by very rapid transverse spatial oscillation of the field direction, while the field lines remain largely unbent up to the scale of the flow. Quantitatively, the fluctuation level and the field-line geometry can be studied in terms of the statistics of the field strength and of the field-line curvature. In the kinematic limit, the distribution of the field strength is an expanding lognormal, while that of the field-line curvature K is stationary and has a power tail K^{-13/7}. The field strength and curvature are anticorrelated, i.e. the growing fields are mostly flat, while the sharply curved fields remain relatively weak. The field, therefore, settles into a reduced-tension state. Numerical simulations demonstrate three essential features of the nonlinear regime. First, the total magnetic energy is equal to the total kinetic energy. Second, the intermittency is partially suppressed compared to the kinematic case, as the fields become more volume-filling and their distribution develops an exponential tail. Third, the folding structure of the field is unchanged from the kinematic case: the anticorrelation between the field strength and the curvature persists and the distribution of the latter retains the same power tail. We propose a model of back reaction based on the folding picture that reproduces all of the above numerical results.

This is a brief review of the main results of our recent studies of the nonlinear evolution of the small-scale MHD dynamo in the high-Prandtl-number regime and of the structure of the resulting saturated state of the isotropic homogeneous MHD turbulence. It is emphasized that the MHD regime without a uniform mean field (as is the case in our studies) is fundamentally different from the one in which such a field is externally imposed. The ability of the turbulence to bend and fold the magnetic-field lines leads to the emergence of a distinctive small-scale structure. The fields are organized in folds of characteristic length comparable to the size of the largest turbulent eddies with spatial-direction reversals at the resistive scale. These folds are very hard to destroy. In the nonlinear regime, the folding structure coexists with Alfven waves propagating along the folds. The turbulent energy injected by the forcing is dissipated in part resistively via the small-scale magnetic fields, and in part viscously via the Alfven waves.

Magnetic field amplification in astrophysics ultimately requires an understanding of magnetohydrodynamic turbulence. Kinetic helicity has long been known to be important for large scale field growth in forced MHD turbulence, and has been recently demonstrated numerically to be asymptotically consistent with slow mean field dynamo action in a periodic box. Here we show numerically that the magnetic spectrum at and below the forcing scale is also strongly influenced by kinetic helicity. We identify a critical value, $f_{h,crit}$ above which the magnetic spectrum develops maxima at wavenumber $= 1$ scale {\it and} at the forcing scale, For $f< f_{h,crit}$ the field peaks only at the resistive scale. Kinetic helicity may thus be important not only for generating a large scale field, but also for establishing observed peaks in magnetic spectra at the forcing scale. The turbulent Galactic disk provides an example where both large scale ($>$ supernova forcing scale) fields and small scale ($\le$ forcing scale, with peak at forcing scale) fields are observed. We discuss this, and the potential application to the protogalaxy, but also emphasize the limitations in applying our results to these systems.

A weak fluctuating magnetic field embedded into a turbulent conducting medium grows exponentially while its characteristic scale decays. In the ISM and protogalactic plasmas, the magnetic Pr is very large, so a broad spectrum of growing magnetic fluctuations is excited at subviscous scales. We study the statistical correlations that are set up in the field pattern and show that the magnetic-field lines possess a folding structure, where most of the scale decrease is due to rapid transverse field direction reversals, while the scale of the field variation along itself stays approximately constant. Specifically, we find that the field strength and the field-line curvature are anticorrelated, and the curvature possesses a stationary limiting distribution with the bulk located at the values of curvature comparable to the characteristic wave number of the velocity field and a power tail extending to large values of curvature. The regions of large curvature, therefore, occupy only a small fraction of the total volume of the system. Our theoretical results are corroborated by direct numerical simulations. The implication of the folding effect is that the advent of the Lorentz back reaction occurs when the magnetic energy approaches that of the smallest turbulent eddies. Our results also directly apply to the problem of statistical geometry of the material lines in a random flow.

We simulate the evolution of an initially weak magnetic field in forced turbulence for a range of Prandtl numbers. The field grows exponentially with the Kulsrud-Anderson $k^{3/2}$ spectrum until the magnetic energy approaches the viscous-scale kinetic energy, where the magnetic forces then backreact on the velocity. Further growth proceeds more slowly until a saturated state is reached where the magnetic and kinetic energies are equal, and where the magnetic energy exists primarily at the resistive scale. We discuss the structure of this turbulence and the extrapolation of the results to astrophysically-large Prandtl numbers.

We introduce Gradient Particle Magnetohydrodynamics (GPM), a new Lagrangian method for magnetohydrodynamics based on gradients corrected for the locally disordered particle distribution. The development of a numerical code for MHD simulation using the GPM algorithm is outlined. Validation tests simulating linear and nonlinear sound waves, linear MHD waves, advection of magnetic fields in a magnetized vortex, hydrodynamical shocks, and three-dimensional collapse are presented, demonstrating the viability of an MHD code using GPM. The characteristics of a GPM code are discussed and possible avenues for further development and refinement are mentioned. We conclude with a view of how GPM may complement other methods currently in development for the next generation of computational astrophysics.

We simulate incompressible MHD turbulence in the presence of a strong background magnetic field. Our major conclusions are: 1) MHD turbulence is most conveniently described in terms of counter propagating shear Alfven and slow waves. Shear Alfven waves control the cascade dynamics. Slow waves play a passive role and adopt the spectrum set by the shear Alfven waves, as does a passive scalar. 2) MHD turbulence is anisotropic with energy cascading more rapidly along k_perp than along k_parallel, where k_perp and k_parallel refer to wavevector components perpendicular and parallel to the local magnetic field. Anisotropy increases with increasing k_perp. 3) MHD turbulence is generically strong in the sense that the waves which comprise it suffer order unity distortions on timescales comparable to their periods. Nevertheless, turbulent fluctuations are small deep inside the inertial range compared to the background field. 4) Decaying MHD turbulence is unstable to an increase of the imbalance between the flux of waves propagating in opposite directions along the magnetic field. 5) Items 1-4 lend support to the model of strong MHD turbulence by Goldreich & Sridhar (GS). Results from our simulations are also consistent with the GS prediction gamma=2/3. The sole notable discrepancy is that 1D power law spectra, E(k_perp) ~ k_perp^{-alpha}, determined from our simulations exhibit alpha ~ 3/2, whereas the GS model predicts alpha = 5/3.