1. Introduction
Only relatively recently have Lagrangian statistics begun to be explored as a tool to understand properties of turbulence in the electrically conducting fluids that are described by magnetohydrodynamics (MHD) (Busse et al. Reference Busse, Müller, Homann and Grauer2007; Homann et al. Reference Homann, Grauer, Busse and Müller2007b; Busse & Müller Reference Busse and Müller2008; Busse, Müller & Gogoberidze Reference Busse, Müller and Gogoberidze2010; Eyink et al. Reference Eyink, Vishniac, Lalescu, Aluie, Kanov, Bürger, Burns, Meneveau and Szalay2013; Pratt et al. Reference Pratt, Busse, Müller, Watkins and Chapman2017; Pratt, Busse & Muller Reference Pratt, Busse and Müller2020a; Pratt, Busse & Müller Reference Pratt, Busse and Müller2020b). In contrast, the study of neutral-fluid turbulence from a Lagrangian viewpoint has a long history (e.g. Taylor Reference Taylor1922; Richardson Reference Richardson1926), with established scaling laws due to Batchelor and Richardson that have been explored both theoretically and experimentally. In neutral fluids, comprehensive studies have examined how Lagrangian statistics related to dispersion are dependent on the Reynolds number of the flow (Yeung & Borgas Reference Yeung and Borgas2004; Yeung, Pope & Sawford Reference Yeung, Pope and Sawford2006; Sawford, Yeung & Hackl Reference Sawford, Yeung and Hackl2008). In this work we explore these questions in the more physically complex system of fully nonlinear incompressible MHD turbulence.
Turbulent mixing is investigated through the diffusion and relative dispersion of fluid tracer particles. A full description of mixing in a plasma depends on the Reynolds number. Relative dispersion in neutral fluids quantifies the mixing of smoke and pollutants in the atmosphere, or of micro-plastics and debris in the oceans. In electrically conducting fluids, diffusion, dispersion, and the mixing processes they represent, determine how fusion products from the core of a star are mixed into the star's outer layers, thereby changing the course of stellar evolution. Relative dispersion also represents the spreading and mixing of plasma in the interstellar or interplanetary medium, and affects how energetic particles and cosmic rays are transported. In these contexts, Reynolds numbers are predicted that are many orders higher than can be achieved presently by direct numerical simulations (DNS) that fully resolve the turbulent motions. However, the features of diffusion and dispersion captured by fluid particles in DNS of MHD turbulence provide information relevant to turbulent mixing in these astrophysical applications (e.g. Zahn Reference Zahn1993; Heyer & Brunt Reference Heyer and Brunt2004; Utomo, Blitz & Falgarone Reference Utomo, Blitz and Falgarone2019). The details of how these statistics change as the Reynolds number increases allow us to make predictions for realistic mixing.
In astrophysical settings, the magnitude of the large-scale magnetic field is often moderate. For example, in the solar wind, ion foreshock and magnetosheath, ranges have been reported such that the magnetic field is between one and 2.5 times the root-mean-square (r.m.s.) fluctuations, i.e. 
$1 \leq B_0 / B_{ {rms}} \leq 2.5$ (see table 1 of Zimbardo et al. Reference Zimbardo, Greco, Sorriso-Valvo, Perri, Vörös, Aburjania, Chargazia and Alexandrova2010). In this work we examine a system with a weak anisotropy caused by such a large-scale magnetic field, and select the situation where the magnetic field is equal to the average r.m.s. magnetic field fluctuations, i.e. 
$B_{ {rms}} \approx B_0$.
This work is structured as follows. In § 2 we describe in detail the simulations performed. In § 3 we introduce a new resolution criterion for anisotropic MHD turbulence simulations, a necessity for producing a set of DNS that can be closely compared. In § 4 we present diffusion curves, a single-particle Lagrangian statistic. In § 5 we present several statistics derived from pairs of tracer particles, focused on understanding the relative dispersion. In § 6 we summarize our findings and draw broader conclusions from the statistics presented.
2. Simulations
We investigate the effect of the Reynolds number on statistically stationary, forced, homogeneous, incompressible MHD turbulence in the presence of a moderate static magnetic field. This magnetic field has a constant value and direction throughout the simulation. Locally the magnetic field also experiences time-dependent fluctuations which possess zero mean. The strength of these fluctuations can be measured with respect to the strength of the imposed magnetic field. A global average of the magnetic field would return the value due to the imposed field and, therefore, the term mean magnetic field best describes the type of field imposed. Such a magnetic field has alternatively been referred to as an external magnetic field, or a guide field, although these terms are less specific to the way that the field has been imposed.
In each DNS, we solve the non-dimensional equations
$$\begin{gather} \frac{\partial \boldsymbol{\omega} }{\partial t} - \boldsymbol{\nabla} \times (\boldsymbol{{v}} \times \boldsymbol{\omega} + \boldsymbol{j} \times \boldsymbol{B}) = \hat{\nu} \nabla^2 \boldsymbol{\omega} + \boldsymbol{f}^{\omega} , \end{gather}$$
$$\begin{gather}\frac{\partial \boldsymbol{B} }{\partial t} - \boldsymbol{\nabla} \times (\boldsymbol{{v}} \times \boldsymbol{B}) = \hat{\eta} \nabla^2 \boldsymbol{B}+ \boldsymbol{f}^{b}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} = 0,\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{v}} = 0, \end{gather}$$
using a pseudospectral method in a simulation volume with periodic boundary conditions. These equations include terms for the solenoidal velocity field 
$\boldsymbol {{v}}$, vorticity 
$\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {{v}}$, magnetic field 
$\boldsymbol {B}$ and current 
$\boldsymbol {j}=\boldsymbol {\nabla }\times \boldsymbol {B}$. Each of the quantities in (2.1)–(2.3) have been non-dimensionalized using relevant time and length scales, commonly referred to as Alfvénic units. Two dimensionless parameters, 
$\hat {\nu }$ and 
$\hat {\eta }$, appear in the equations. They derive from the kinematic viscosity 
$\nu$ and the magnetic diffusivity 
$\eta$. A fixed time step and a low-storage third-order Runge–Kutta method (Williamson Reference Williamson1980) are used for the time integration. The mean magnetic field is designated by 
$B_0$, and points purely in the positive 
$z$-direction. At any point in time it has a value close to unity with respect to the r.m.s. of the magnetic field fluctuations; a long-time average produces 
$B_0 \approx B_{ {rms}}$. The Alfvén ratio is approximately unity for all simulations discussed in this work; we calculate this ratio as 
$r_{{A}} = \langle \textit {E}_{ {v}} /\textit {E}_{ {b}} \rangle$ from the kinetic energy per unit mass 
$\textit {E}_{ {v}}= \boldsymbol {{v}}^2/ 2$ and the energy per unit mass contained in the magnetic fluctuations 
$\textit {E}_{ {b}}= B_{ {rms}}^2 / 2$. In the Alfvén ratio, the brackets indicate an average over time; this is the simulation time during which the dynamics of Lagrangian tracer particles are examined. The implications of an Alfvén ratio of one are that our simulations include the full nonlinear interaction of velocity and magnetic fields that contribute with equal weight.
To maintain the turbulence in a statistically stationary steady state, the vorticity and magnetic fields are forced on the largest scales of the simulation volume using a deterministic forcing, applied in Fourier space, that also allows the largest-scale motions of the system to evolve. Deterministic forcing has the advantage that no stochastic source of fluctuations is introduced on large scales. We call the deterministic forcing method that we use homogeneous forcing; it is distinct from forcing methods used in several earlier works on Lagrangian MHD turbulence (Busse et al. Reference Busse, Müller, Homann and Grauer2007; Homann et al. Reference Homann, Grauer, Busse and Müller2007b; Busse & Müller Reference Busse and Müller2008; Busse et al. Reference Busse, Müller and Gogoberidze2010) but identical to the forcing method used in Pratt et al. (Reference Pratt, Busse and Müller2020b). Homogeneous forcing establishes a constant injection of energy at large scales. In (2.1) and (2.2) forcing terms 
$\boldsymbol {f}^{\omega }$ and 
$\boldsymbol {f}^{b}$ are introduced which are non-zero only for the wavevector shells 
$1 \leq |\boldsymbol {k}| \leq 2.5$. The forcing terms are defined as
$$\begin{gather} \boldsymbol{\hat{f}}^{\omega}(\boldsymbol{k},t) = \gamma_{f,\omega} \frac{\hat{\boldsymbol{\omega}}(\boldsymbol{k},t)}{|\hat{\boldsymbol{\omega}}(\boldsymbol{k},t)|^2} , \end{gather}$$
$$\begin{gather}\boldsymbol{\hat{f}}^{b}(\boldsymbol{k},t) = \gamma_{f,b} \frac{\hat{B}(\boldsymbol{k},t)}{|\hat{B}(\boldsymbol{k},t)|^2}. \end{gather}$$
Variables with hats are used to describe the forcing because it is applied in Fourier space. The constants 
$\gamma _{f, \omega }$ and 
$\gamma _{f, b}$ regulate the energy injection. These two forcing constants are set equal in our simulations, and are identical across all of the simulations presented in this work. By forcing both fields equally, we achieve approximate equipartition of kinetic and magnetic energy at all scales, and study the regime of MHD turbulence defined by the full nonlinear interaction of the velocity and magnetic fields. This is in contrast to studies that force only the velocity field, oriented toward understanding the dynamo (e.g. McKay et al. Reference McKay, Linkmann, Clark, Chalupa and Berera2017; Brandenburg et al. Reference Brandenburg, Haugen, Li and Subramanian2018). In natural settings such as molecular clouds (Hennebelle & Falgarone Reference Hennebelle and Falgarone2012; Heiles & Troland Reference Heiles and Troland2005), and the solar wind (Müller & Grappin Reference Müller and Grappin2005; Boldyrev, Perez & Zhdankin Reference Boldyrev, Perez and Zhdankin2012) kinetic and magnetic energies are commonly observed to be close to equipartition, so this choice is a physically realistic one.
To inhibit the emergence of states dominated by Elsässer (Elsässer Reference Elsässer1950) positive (
$z^+$) or negative (
$z^-$) interactions, the cross helicity of the forced modes is set to zero. This is accomplished as part of the forcing scheme by enforcing orthogonality between the magnetic field vector and the velocity field vector (see, e.g. Müller, Malapaka & Busse Reference Müller, Malapaka and Busse2012). The magnitude of the total cross helicity, normalized by 
$\sqrt {\textit {E}_{ {v}}}\sqrt {\textit {E}_{ {b}}}$, never rises above 
$0.095$ in the simulations examined in this work. This prevents the MHD turbulent system from becoming imbalanced, which can lead to a breakdown of the nonlinear energy cascade (as discussed in Biskamp Reference Biskamp2003). In addition, the level of total magnetic helicity, normalized by 
$\textit {E}_{ {b}}^2/k_{ {forcing}}$, never rises above 
$5 \times 10^{-4}$. For a system in a quasi-stationary state, homogeneous forcing is expected to disturb the natural turbulent flow only mildly; this forcing method has been examined in Ghosal et al. (Reference Ghosal, Lund, Moin and Akselvoll1995) and Vorobev et al. (Reference Vorobev, Zikanov, Davidson and Knaepen2005). Based on the forcing, our simulations are consistent with strong turbulence as described by Perez & Boldyrev (Reference Perez and Boldyrev2007) and Verdini & Grappin (Reference Verdini and Grappin2012). Large-scale Alfvén waves are permitted and are observed when homogeneous forcing is used.
2.1. Lagrangian tracer particles
The positions of Lagrangian tracer particles are initialized in a homogeneous random distribution at a time when the turbulent flow has attained a statistically stationary state. The approximate number of tracer particles deployed in each simulation, in millions, is given in table 1. The particle numbers we use are comparable to earlier MHD turbulence studies (Busse et al. Reference Busse, Müller, Homann and Grauer2007; Homann et al. Reference Homann, Grauer, Busse and Müller2007b; Busse & Müller Reference Busse and Müller2008; Pratt et al. Reference Pratt, Busse and Müller2020b), and are larger than those used in Yeung et al. (Reference Yeung, Pope and Sawford2006) and Sawford et al. (Reference Sawford, Yeung and Hackl2008) to study homogeneous isotropic hydrodynamic turbulence. The consequence is that the Lagrangian statistics that we produce in this work are well resolved in space. At each time step the particle velocities are interpolated from the instantaneous Eulerian velocity field using a tricubic polynomial interpolation scheme (Lekien & Marsden Reference Lekien and Marsden2005; Homann, Dreher & Grauer Reference Homann, Dreher and Grauer2007a). Particle positions are calculated by numerical integration of the equation of motion using a low-storage third-order Runge–Kutta method that matches the one used for the Eulerian field integration. We record particle information every four time steps, at approximate time intervals of 
$0.1 \tau _{\eta }$, so that the Lagrangian statistics produced are also well resolved in time. The hydrodynamic simulation that we perform for comparison, simulation 3H, uses a slightly larger time interval of 
$0.2 \tau _{\eta }$. Particles are initially arranged in tetrads with a reference particle at the vertex and three other particles aligned along each Cartesian direction at fixed separations from the reference particle. Each dispersion result is defined by this initial separation distance. The smallest initial separation distances that we deploy are well resolved by the grids used in each simulation; this will be quantified in the section on the resolution. Each simulation is run for at least 400 of the Kolmogorov time scale 
$\tau _{\eta }$. At large times, the particles that make up a pair can reach a separation that is on the scale of the simulation volume. For pairs of particles that have separated that far, the probability is low for a statistically significant number to reapproach each other again. Although the forcing is applied on length scales comparable to the simulation volume, the velocity field experienced by pairs of Lagrangian tracer particles separated by the size of the simulation volume is generally decorrelated, and the two particles move approximately independently from each other.
Simulation parameters: the simulation number, the Eulerian grid size 
$N^3$ and the number of tracer particles 
$N_p$ is provided for each simulation. The Kolmogorov time scale 
$\tau _{\eta }$, large-eddy length scale 
$\textit {L}_{ {{E}}}$, the large-eddy turnover time 
$\textit {T}_{ {{E}}}$ and the time-averaged r.m.s. of magnetic fluctuations 
$B_{ {rms}}$ are given. The magnitude of the mean magnetic field is 
$B_0=1$ for case 
$1,2,3,4$. The resolution in the perpendicular direction 
$k_{ {{max}},\perp } \eta _{ {{kol}},\perp }$ and in the parallel direction 
$k_{ {{max}},\parallel } \eta _{ {{kol}},\parallel }$ are provided. The Reynolds number is calculated as described in (3.13) using the perpendicular Kolmogorov microscale 
$\eta _{ {{kol}},\perp }$. All simulations take place in a cubic simulation volume, and flow statistics are gathered for at least 
$400 \tau _{\eta }$. The simulation 3H is a purely hydrodynamic simulation performed for comparison with simulation 3.
3. Resolution of anisotropic MHD turbulence
In anisotropic MHD turbulence, attention must be given both to the shape of the simulation volume and the ability to resolve the smallest relevant scales of the turbulent dynamics. At both large scales and small scales, differences exist between the energy spectra in the directions parallel and perpendicular to the mean magnetic field (see figure 1). The asymptotic spectral scaling exponents of the one-dimensional energy spectra are approximately -1.5 (perpendicular) and 
$-1.6$ (parallel). These spectra are defined as 
$E(k)=\int \,{\rm d}^3\boldsymbol {k}' E(\boldsymbol {k}') \delta (|\boldsymbol {k}'\boldsymbol {\cdot} \hat {\boldsymbol {e}}_k|-k)$, where the wavenumber 
$k$ runs along the direction given by the appropriately chosen and fixed unit vector 
$\hat {\boldsymbol {e}}_k$. These scalings are consistent with the current understanding of inertial-range dynamics of incompressible MHD turbulence, which is based on the concepts of critical balance of characteristic time scales parallel and perpendicular to the local magnetic field (Goldreich & Sridhar Reference Goldreich and Sridhar1995b; Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2015) and dynamical alignment of fluctuations of velocity and magnetic field (Boldyrev Reference Boldyrev2006), in combination with specific log-Poisson intermittency corrections (Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015; Mallet & Schekochihin Reference Mallet and Schekochihin2017). We provide these scalings purely for the sake of completeness; they should be viewed with care. The main caveat is the limited width of the inertial scaling range, a consequence of the well-resolved dissipation region required for Lagrangian small-scale statistics. Another relevant issue is that anisotropic MHD energy transfer for systems with moderate mean magnetic fields should be analysed in a local frame of reference aligned with the local mean magnetic field, which includes contributions from magnetic fluctuations, instead of the global mean field 
$B_0$. Such an analysis of spectral transfer is outside the scope of this Lagrangian investigation. For the present work, the following simple observations of MHD spectral dynamics seem sufficient to characterize the system investigated. As the energy cascade proceeds to smaller scales, the eddies become more elongated in the direction parallel to the magnetic field; the amount of elongation is dependent on the length scale of motions in that direction, which in turn is dependent on the strength of the magnetic field (e.g. as discussed in Cho & Vishniac Reference Cho and Vishniac2000; Schekochihin, Cowley & Yousef Reference Schekochihin, Cowley and Yousef2008; Verdini et al. Reference Verdini, Grappin, Hellinger, Landi and Müller2015). We explored the corresponding time-scale dependence of this anisotropy from the Lagrangian point of view in Pratt et al. (Reference Pratt, Busse and Müller2020b). Due to the elongation of eddies in the direction aligned with the magnetic field, the smallest length scales relevant for the study of turbulence are different from those in the perpendicular direction. At the largest scales, correlation lengths are also longer in the direction aligned with the mean magnetic field than in the perpendicular direction.
Time-averaged kinetic energy spectra for the four MHD simulations described in table 1. The spectra are calculated for a one-dimensional wavevector 
$k$ taken in the (a) 
$x$-direction (perpendicular to the mean magnetic field), and (b) 
$z$-direction (parallel to the mean magnetic field). Grey lines indicate theoretical scaling laws relevant to the inertial range for MHD turbulence.
3.1. Resolution of the largest scales of anisotropic MHD turbulence
Corresponding to these correlation lengths, velocity structures have different characteristic scales in the directions parallel and perpendicular to the mean magnetic field. To assure that these velocity structures are both resolved, we consider how the periodic simulation volume should be shaped. We solve the MHD equations (2.1)–(2.3) in a simulation volume with sides of length 
$2 {\rm \pi}$ in the 
$x$ and 
$y$ directions, that are perpendicular to the mean magnetic field. Using a simulation volume that is elongated in the 
$z$-direction is common for simulations of anisotropic MHD that examine a strong mean magnetic field. To determine the necessary length in the 
$z$-direction, we consider the correlation length of the velocity field in each direction. In test simulations with a stronger mean magnetic field, we measure a correlation length of the velocity field 
$L_{ {c},\parallel }$ in the direction parallel to 
$B_0$ that is significantly larger than in the perpendicular direction. This measurement agrees with previous results (e.g. Cho, Lazarian & Vishniac Reference Cho, Lazarian and Vishniac2002; Boldyrev Reference Boldyrev2005; Chandran Reference Chandran2008). For the strong turbulence regime, a relative change in length scales can be predicted from the premise of critical balance (Goldreich & Sridhar Reference Goldreich and Sridhar1995a). Using this premise, the ratio of the largest-scale wavenumbers (
$k \sim 2 {\rm \pi}/ L$) in the perpendicular and parallel directions grows linearly with the magnitude of the mean magnetic field 
$B_0$,
\begin{equation} k_{{\parallel}} B_0 \sim k_{{\perp}} B_{{rms}}. \end{equation}
Thus, for the largest scales of the flow, critical balance predicts that the ratio of parallel and perpendicular length scales should grow linearly with 
$B_0$. The critical balance relation in (3.1) implies that the nonlinear eddy turnover time is of the order of the Alfvén time 
$\tau _{ {{NL}}} \sim \tau _{ {{A}}}$. For the 
$B_0=1$ simulations examined in this work, the box length in the parallel direction is much larger than the parallel correlation length, 
$L_z \gtrsim 10 L_{ {c},\parallel }$. We therefore are able to use a cubic simulation volume with no elongation for all simulations in table 1.
3.2. Resolution of the smallest scales of anisotropic MHD turbulence
We revisit the resolution criterion commonly used to measure whether isotropic turbulence is sufficiently resolved, and we extend this criterion to the anisotropic case. The smallest length and time scales that characterize Navier–Stokes turbulence are defined in terms of the kinetic energy dissipation rate 
$\epsilon _{{v}}$ and the kinematic viscosity. These are the Kolmogorov microscale 
$\eta _{ {{kol}}}=(\nu ^3/\epsilon _{ {v}})^{1/4}$ and the Kolmogorov time scale 
$\tau _{\eta }=(\nu /\epsilon _{ {v}})^{1/2}$. To test whether a DNS adequately resolves the smallest physically relevant spatial scales of turbulence, the Kolmogorov microscale 
$\eta _{ {{kol}}}$ is typically multiplied by the highest wavenumber resolved in a simulation 
$k_{{{max}}}$. This provides a non-dimensional number that indicates how well the grid spacing resolves this smallest physical length scale. The ‘standard’ criterion for adequate spatial resolution for a simulation in the case of homogeneous isotropic turbulence (Yeung & Pope Reference Yeung and Pope1989; Pope Reference Pope2000; Donzis, Yeung & Sreenivasan Reference Donzis, Yeung and Sreenivasan2008; Yeung, Sreenivasan & Pope Reference Yeung, Sreenivasan and Pope2018) is then
\begin{equation} k_{{{max}}} \eta_{{{kol}}} \gtrapprox 1.5. \end{equation}The presence of a mean magnetic field makes the gradients of the turbulent fields higher in the perpendicular directions compared with the gradients in the parallel direction. Therefore, we need to extend the criterion in (3.2) to ensure that the smallest eddies are well resolved both for the mean field parallel and perpendicular directions.
To extend this resolution criterion, we examine the definition of the kinetic energy dissipation rate, which is related to gradients in the velocity field
\begin{equation} \epsilon_{v}= \frac{1}{2}\sum_{i,j}\nu \left\langle \left(\frac{\partial {v}_i}{\partial x_j}\right)^2\right\rangle . \end{equation}In Fourier space this is expressed as
\begin{equation} \epsilon_{v}= \sum_{\boldsymbol{k}} \nu |\boldsymbol{k}|^2 \sum_{i=1}^{3}\langle \hat{{v}}_i^{{\ast}}(\boldsymbol{k},t) \hat{{v}}_i(\boldsymbol{k},t)\rangle. \end{equation}
By separating the modulus of the wavevector in (3.4) into components parallel 
$k_{\parallel }=\boldsymbol {k}\boldsymbol {\cdot} (B_0 \boldsymbol {\hat {z}})/B_0$ and perpendicular 
$\boldsymbol {k}_{\perp }=\boldsymbol {k}\times (B_0 \boldsymbol {\hat {z}})/B_0$ to the direction of anisotropy, the kinetic energy dissipation rate can be split into contributions that arise from gradients parallel and perpendicular to the mean magnetic field. We define these components as
$$\begin{gather} \epsilon_{v} =\frac{2}{3}\epsilon_{{v},\perp} + \frac{1}{3}\epsilon_{{v},\parallel}, \end{gather}$$
$$\begin{gather}\epsilon_{{v},\perp}=\frac{3}{2} \sum_{\boldsymbol{k}} \nu |\boldsymbol{k}_{{\perp}}|^2 \sum_{i=1}^{3}\langle \hat{{v}}_i^{{\ast}}(\boldsymbol{k},t) \hat{{v}}_i(\boldsymbol{k},t)\rangle, \end{gather}$$
$$\begin{gather}\epsilon_{{v},\parallel} = 3 \sum_{\boldsymbol{k}} \nu k_{{\parallel}}^2 \sum_{i=1}^{3}\langle \hat{{v}}_i^{{\ast}}(\boldsymbol{k},t) \hat{{v}}_i(\boldsymbol{k},t)\rangle. \end{gather}$$
In the case of homogeneous isotropic turbulence, these definitions recover 
$\epsilon _{ {v},\perp }= {\epsilon _{ {v},\parallel }= \epsilon _{ {v}}}$. In the case of a finite mean magnetic field, gradients perpendicular to the mean magnetic field will be larger than gradients parallel to the mean magnetic field, leading to 
$\epsilon _{ {v},\perp }> \epsilon _{ {v}}$ and 
$\epsilon _{ {v},\parallel }< \epsilon _{ {v}}$.
Two different length scales, for the direction perpendicular and parallel to the mean magnetic field, can be defined from 
$\epsilon _{ {v},\perp }$ and 
$\epsilon _{ {v},\parallel }$. These are
$$\begin{gather} \eta_{{{kol}},\perp}=\left(\hat{\nu}^3/\epsilon_{{v},\perp} \right)^{1/4}, \end{gather}$$
$$\begin{gather}\eta_{{{kol}},\parallel}=\left(\hat{\nu}^3/\epsilon_{{v},\parallel} \right)^{1/4} . \end{gather}$$
In the isotropic case, these length scales are equal to the Kolmogorov length scale. Using 
$\eta _{ {kol},\perp }$ and 
$\eta _{ {kol},\parallel }$ a generalized version of the small-scale resolution criterion in (3.2) can be defined for homogeneous anisotropic turbulence. This combined criterion for adequate resolution of the smallest scales both parallel and perpendicular to the mean magnetic field direction is
$$\begin{gather} k_{{{max}},\perp} \eta_{{{kol}},\perp} \gtrapprox 1.5 , \end{gather}$$
$$\begin{gather}k_{{{max}},\parallel} \eta_{{{kol}},\parallel} \gtrapprox 1.5 . \end{gather}$$
These two criteria are fulfilled in all simulations discussed in this work, and their values are listed in table 1. When we use these criteria to set-up our simulations, the smallest scales of anisotropic MHD turbulence appear to be well resolved in both directions. For example, Perez et al. (Reference Perez, Mason, Boldyrev and Cattaneo2014) notes that the numerical effects of insufficient small-scale resolution are a steepening of the spectrum at intermediate scales and a flattening closer to the grid scale. Those effects are not observed in the spectra of our simulations in figure 1. The smallest of initial separation distances for pairs of Lagrangian tracer particles are set to 
$2 \eta _{ {{kol}},\perp }$. Because the Kolmogorov length scale in this direction is well resolved, so are these particle separations.
3.3. Reynolds number for anisotropic MHD turbulence
For a homogeneous isotropic system, the Reynolds number is standardly defined from the kinetic energy, the viscosity and a characteristic length scale
\begin{equation} {{\it Re}}= \langle \textit{E}_{{v}}^{1/2} \textit{L}_{{{E}}}\rangle /\nu. \end{equation}
The characteristic length scale 
$\textit {L}_{ {{E}}}$ is defined as a dimensional estimate of the size of the largest eddies, 
$\textit {L}_{ {{E}}}=\textit {E}_{ {v}}^{3/2}/\epsilon _{ {v}}$. This length scale is summarized in table 1 for each of our simulations.
To calculate the Reynolds number for anisotropic flows, we use the more general definition of the Reynolds number (see Chapter 6.1.2 of Pope Reference Pope2000)
\begin{equation} {{\it Re}}= c (\eta_{{{kol}},\perp} / \textit{L}_{{{F}}})^{{-}4/3}, \end{equation}
where 
$\textit {L}_{ {{F}}}$ is a forcing length scale and 
$c$ is a constant that must be determined. Our method of homogeneous forcing affects a minimum length scale
\begin{equation} \textit{L}_{{{F}}}=2 {\rm \pi}/ k_{{{F}},{{max}}} = 0.8 {\rm \pi}. \end{equation}
We determine the constant 
$c$ by comparing the definitions of the Reynolds number in (3.12) with that in (3.14) for an isotropic 
$B_0=0$ simulation; this produces a value of the constant 
$c \approx 1.16$. The Reynolds number calculated using (3.13) is summarized in table 1 for each of our simulations. For an isotropic flow, the standard Reynolds number based on the Kolmogorov microscale, 
${{\it Re}}$, and the Taylor-scale Reynolds number, 
${{\it Re}}_{\lambda }$, provide equivalent information (see, for example, Pope Reference Pope2000; Sawford et al. Reference Sawford, Yeung and Hackl2008), and are related by
\begin{equation} 15 {{\it Re}} = {{\it Re}}_{\lambda}^2. \end{equation}
Thus, the Reynolds numbers in table 1 can be translated to Taylor-scale Reynolds numbers; for simulation 4, we find that 
${{\it Re}}_{\lambda } \approx 560$. The magnetic Reynolds number is defined from the Reynolds number and the magnetic Prandtl number, i.e. 
$\textit {Re}_{ {m}}=\textit {Pr}_{ {m}} \textit {Re}$. In all simulations in this work, the magnetic Prandtl number 
${Pr}_{ {m}}=1$ so that the magnetic Reynolds number is equal to the Reynolds number.
Table 1 provides an overview of the Reynolds number and other fundamental parameters. Among the fundamental parameters in this table are also two time scales: the Kolmogorov time scale 
$\tau _{\eta }$, which is associated with the smallest scales of motion, and a time scale associated with the largest scales of motion, called the large-eddy turnover time 
$\textit {T}_{ {{E}}} =\textit {E}_{ {v}}/\epsilon _{ {v}}$. Our measurements for the ratio of these time scales, 
$\tau _{\eta }/\textit {T}_{ {{E}}}$, decrease as 
${{\it Re}}^{-1/2}$, as expected for isotropic hydrodynamic turbulence (see Chapter 6.1.2 of Pope Reference Pope2000).
4. Single-particle Lagrangian Diffusion
The examination of single-particle diffusion curves, the average square distance that a particle has moved from its initial position, is common in Lagrangian studies of turbulence. We designate this quantity by 
$\langle \Delta ^2(t) \rangle$, where the brackets indicate an average over all Lagrangian tracer particles. In anisotropic MHD turbulence, such diffusion curves were first studied by Busse & Müller (Reference Busse and Müller2008), who examined the influence of mean magnetic fields of magnitude 
$B_0 = 2$ and 5 in units of 
$B_{ {rms}}$. Our present results are physically distinct from that work in three ways: (1) we examine a significantly weaker mean magnetic field, (2) we use a different forcing method, and (3) our simulations extend to a higher Reynolds number.
If diffusion curves are calculated from a single initial time, they can be influenced by idiosyncrasies of the flow at that time. We produce our diffusion curves by averaging the results for several independent initial times using the conventional averaging methods described in Dubbeldam et al. (Reference Dubbeldam, Ford, Ellis and Snurr2009), so that they are not dependent on any single initial state of the flow. Diffusion, as a single-particle statistic, is also notoriously sensitive to large-scale flow features because each particle's separation is measured relative to a fixed point in space. Large-scale flows are able to sweep along large numbers of Lagrangian particles, affecting the outcome of diffusion curves. Because of this large-scale sweeping, diffusion statistics cannot be formulated to provide information about a specific length scale and the related time scale.
Since diffusion is dominated by the largest-scale fluctuations of the system that carry most of the kinetic energy, the natural time scale to normalize the time is the large-eddy turnover time 
$\textit {T}_{ {E}}$. For a fixed integral scale of turbulence, the Kolmogorov time decreases as the extension of the inertial range increases. In a statistically isotropic Navier–Stokes setting, dimensional analysis relates this large-scale quantity to the Kolmogorov time scale 
$\textit {T}_{ {{E}}} \sim \tau _{\eta } {{\it Re}}^{1/2}$ (Pope Reference Pope2000; Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009). This relation is responsible for the later arrival of particles in the diffusive regime with increasing Reynolds number reported in Sawford (Reference Sawford1991). In the MHD case, we reproduce this later arrival of particles in the diffusive regime with increasing Reynolds number when we normalize time by 
$\tau _{\eta }$. However, due to the dependence of 
$\tau _{\eta }$ on Reynolds number, normalizing by the Kolmogorov time confuses a physical interpretation of the diffusion curves. Using 
$\textit {T}_{ {{E}}}$ instead allows for a close comparison of the diffusion curves, which collapse to show a universal trend (see figure 2). At short times, our single-particle diffusion curves demonstrate the expected ballistic scaling with time to the power of two. At long times, a diffusive regime is expected; for homogeneous isotropic turbulence, the curves scale as a power of one in this regime.
Average square separation of particles from their initial position 
$\langle \Delta ^2(t) \rangle$, for (a) separations perpendicular to the mean magnetic field and (b) separations aligned with the mean magnetic field. Each curve represents an average over at least three independent initial times. Distance is measured in units of the large-eddy length scale 
$\textit {L}_{ {{E}}}$ and time is measured in units of the large-eddy turnover time 
$\textit {T}_{ {{E}}}$. The straight black lines indicate the scaling laws that are theoretically predicted, with the scaling exponent labelled.
A close examination of the diffusive scaling in our simulations (see figure 3) reveals oscillations in the log derivative curves in the diffusive regime. This derivative is calculated using a simple forward Euler method. Fluctuations are present in homogeneous isotropic turbulence; however, the oscillations that we observe in anisotropic MHD turbulence tend to be more regular and noisier. In the perpendicular direction a scaling with an average value close to unity has been established. In the parallel direction, this is also true for the lower Reynolds number simulations. The highest Reynolds number simulation continues to have an average slope that is clearly steeper than one throughout the simulation time; at no point is this slope as low as one. This measurement of parallel diffusion is calculated parallel to 
$B_0$ rather than parallel to the local mean field; particularly in the present case where 
$B_0$ is approximately equal to the fluctuations in magnetic field, that difference may be significant, so that a deeper interpretation of this anisotropy becomes difficult.
Derivative of the log of 
$\langle \Delta ^2(t) \rangle$ as in figure 2 for (a) separation perpendicular to the mean magnetic field and (b) separation aligned with the mean magnetic field. Each diffusion curve represents an average over at least three independent initial times. Distance is measured in units of the Kolmogorov length scale 
$\eta _{ {{kol}},\perp }$ and time is measured in units of the Kolmogorov time scale 
$\tau _{\eta }$. A grey line indicates the theoretical prediction that the diffusion curve scales linearly with time.
5. Lagrangian statistics for particle pairs
Large-scale flow features, in the present case Alfvénic fluctuations (e.g. Howes Reference Howes2015), can dominate Lagrangian statistics based on single particles. We therefore focus our attention on diagnostics based on pairs of Lagrangian tracer particles to expose characteristics of the smaller scales of turbulence.
5.1. Lagrangian two-particle dispersion
The most commonly examined pair statistic is two-particle dispersion. Two-particle dispersion is the separation of a pair of particles relative to each other, and is usually expressed as a mean-square displacement. The separation of a pair of Lagrangian tracer particles, labelled 
$i$ and 
$j$, is simply 
$\boldsymbol {\xi } = \boldsymbol {r}_i(t) - \boldsymbol {r}_j(t)$, where 
$\boldsymbol {r}_i$ is the vector position of particle 
$i$ in three-dimensional space. Dispersion is typically calculated as 
$\langle (\boldsymbol {\xi } - \boldsymbol {\xi }_0)^2 \rangle$, where the angular brackets denote an average over all particle pairs that have an initial separation 
$\boldsymbol {\xi }_0$. At time 
$t=t_0$, the quantity 
$\langle (\boldsymbol {\xi } - \boldsymbol {\xi }_0)^2 \rangle$ is identically zero. Three subranges are theoretically predicted to exist for isotropic hydrodynamic turbulence,
\begin{equation} \langle (\boldsymbol{\xi} - \boldsymbol{\xi}_0)^2 \rangle \sim \begin{cases} t^2, & \text{ballistic regime}, \\ t^3, & \text{Richardson regime}, \\ t, & \text{diffusive regime}. \end{cases} \end{equation}These three predictions are relevant to short separation times, intermediate separation times and long separation times, respectively. The ballistic regime and diffusive regime have been theoretically motivated, and confirmed by simulations and experiments for hydrodynamic turbulence; simulations have also confirmed that these two regimes exist for isotropic MHD turbulence. The Richardson scaling is a plausible prediction for high-Reynolds-number hydrodynamic turbulence (for example, see figure 5 of Bourgoin Reference Bourgoin2015), and assumes that the initial separation of the pair of particles is arbitrarily small. For details of the derivations of these scaling laws and further work to improve them, we refer to the reviews of Salazar & Collins (Reference Salazar and Collins2009) and Sawford (Reference Sawford2001).
We examine dispersion curves for pairs of particles that are initially separated in either the direction perpendicular or parallel to the mean magnetic field; we call these ‘perpendicular pairs’ and ‘parallel pairs’, respectively. We also calculate the distances that the particle pairs separate in the perpendicular direction and the parallel direction. Because a pair of particles has an initial separation that is small compared with the large-scale flow features, two-particle statistics avoid the influence of large-scale sweeping. The most relevant time scale for measuring two-particle dispersion is therefore the Kolmogorov time scale 
$\tau _{\eta }$ rather than the large-eddy turnover time 
$\textit {T}_{ {{E}}}$. Dispersion curves for pairs initially separated by 
$2 \eta _{ {{kol}},\perp }$ are displayed in figure 4. This is the closest initial separation for the particle pairs that we have followed; larger initial separations and directions also display a scaling regime with a clear slope but the transitions are less sharp. The dispersion curves of our four simulations appear nearly identical at short times until approximately 
$10 \tau _{\eta }$, a range in time corresponding roughly to the ballistic regime; the length of this range appears similar in parallel and perpendicular directions. During the ballistic regime, faster relative dispersion is observed in the perpendicular direction shown in figure 4(a) than in the parallel direction shown in figure 4(b). At long times approaching 
$400 \tau _{\eta }$, the scaling of the dispersion curves is steeper than one. For the three simulations with the lowest Reynolds number (sim. no. 1–3), the single-particle diffusion curves have established a clear diffusive scaling at these late times, characterized by a mildly superdiffusive trend.
Average square separation of particle pairs initially separated by 
$\xi _0=2 \eta _{ {{kol}},\perp }$ for (a) perpendicular pairs and separation measure perpendicular to the mean magnetic field, and (b) parallel pairs and separation measured in the direction aligned with the mean magnetic field. Reference scalings for the ballistic regime, Richardson regime and diffusive regime are included as straight black lines, with the scaling exponent labelled.
Following Salazar & Collins (Reference Salazar and Collins2009), we refer to the middle range of scales, between approximately 
$20 \tau _{\eta }$ and 
$100 \tau _{\eta }$, as the inertial subrange without implying that the Richardson scaling or any other scaling with time is a correct prediction for anisotropic MHD turbulence. For simulations 1 and 2 in table 1, which have lower Reynolds numbers, the dispersion curves display a mild curvature that prevents the clear definition of a scaling during the inertial subrange of time scales. For the higher Reynolds number simulations numbered 3 and 4, the dispersion curves flatten during the inertial subrange and approach a scaling prediction that is clearer. In the perpendicular direction, for the pairs with initial separation 
$\xi _0=2 \eta _{ {{kol}},\perp }$ shown in figure 4, this scaling approaches 3, a number reminiscent of the Richardson prediction. In the parallel direction the equivalent scaling appears to be steeper than 3.
We examine the log derivative, to quantitatively analyse these scalings (see figure 5). Here the derivative of the log is calculated using a simple forward Euler method, because such a two-point method allows the very early behaviour to be seen most clearly. The log derivative shows an initial slope of 2 during the ballistic regime for both perpendicular and parallel results. The inertial subrange and diffusive regime are both characterized by chatter in the log derivative. During the inertial subrange, between roughly 
$20 \tau _{\eta }$ and 
$100 \tau _{\eta }$, the log derivative of the perpendicular dispersion curve for perpendicularly separated pairs ranges between 2.3 and 3.7, while the parallel dispersion of pairs separated in the parallel direction ranges between approximately 3.0 and 4.0. The range over which these derivatives chatter provides an estimate of the uncertainty for the scaling of the dispersion curves. Some of the noisiness here results from the fact that our two-particle statistics are constructed from a single initial time. During the diffusive regime both perpendicular and parallel dispersion curves decay to values between approximately 1.0 and 2.0. Simulations 1 and 2 have an average scaling greater than 1 and less than 1.5 between 
$300 \tau _{\eta }$ and 
$400 \tau _{\eta }$, corresponding to a slightly superdiffusive separation.
Derivative of the log of the average square separation of particle pairs for (a) perpendicular pairs and separation measured perpendicular to the mean magnetic field, and (b) parallel pairs and separation measured in the direction aligned with the mean magnetic field, as in figure 4. The initial separation of particle pairs is 
$\xi _0=2 \eta _{ {{kol}},\perp }$. A grey line indicates the theoretical prediction that the dispersion curve scale with the square of time at early times.
As we have discussed, for particle pairs that are initially separated by 
$\xi _0=2 \eta _{ {{kol}},\perp }$, the separation in the perpendicular direction appears close to the Richardson prediction of a scaling of 3. Calculating the log derivative for particle pairs with different initial separations 
$\xi _0$ in simulation 4 (see figure 6) clarifies that this scaling is indeed dependent on 
$\xi _0$ (as discussed by Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005). This figure also provides a comparison to the hydrodynamic simulation 3H. For anisotropic MHD turbulence, the log derivative reveals larger changes between the ballistic and the inertial subranges than in the hydrodynamic case. For each group of particle pairs with the same initial separation, there is first a dip before or near 
$\tau _{\eta }$ indicating a temporary slowing down of dispersion. This dip is followed by a series of peaks as the pair separation enters the inertial subrange; since the hydrodynamic simulation has a single smooth peak, the multiple peaks and the chatter during this period are likely to result from Alfvénic fluctuations. The maximum value of the log derivative for anisotropic MHD turbulence is a higher value than the hydrodynamic case. Thus, although the dispersion curves for simulation 4 point toward a Richardson-like scaling regime for the perpendicular dispersion, they do not clearly confirm Richardson scaling for anisotropic MHD turbulence. The log derivative curves indicate a larger slope is achieved during the inertial subrange for a smaller initial separation 
$\xi _0$. Therefore, in the limit where 
$\xi _0$ becomes arbitrarily small, we expect that the slope of the dispersion curves will be greater than 3 in both the perpendicular and parallel directions. We thus also expect that corrections to a Richarson-like prediction are needed for the case of anisotropic MHD turbulence.
Log derivative of the average square separation of particle pairs in the direction perpendicular to the mean magnetic field for (a) perpendicular pairs and separation measured perpendicular to the mean magnetic field, and (b) parallel pairs and separation measured in the direction aligned with the mean magnetic field. Each line is labelled by the initial separation distance of the particle pairs. Data from simulation 4, described in table 1. Equivalent curves from simulation 3H are shown as dashed lines in the background for comparison. A grey line indicates the theoretical prediction that the dispersion curve scale with the square of time at early times.
It has been shown (e.g. Yeung & Borgas Reference Yeung and Borgas2004) that the intermittency of particle-pair dispersion is larger in higher Reynolds number simulations of isotropic hydrodynamic turbulence. In hydrodynamic turbulence, increasing the Reynolds number leads to a stronger intermittency at the small scales; the separation of particle pairs provides a convenient measure for this, since they sample the velocity field on a length scale comparable to their separation distance. The skewness, a normalized third moment, indicates the asymmetry of the wings of a distribution; it is therefore one indicator of intermittency, which can also be observed in other high-order moments. A larger skewness of the distribution of particle-pair separations indicates the importance of the extremes of dispersion, and the relative importance of pairs of particles that separate much faster than the average. In our simulations of anisotropic MHD turbulence, the skewness of the separation distance provides a particularly clear result that demonstrates the Reynolds number dependence (see figure 7). In the figure the particle pairs are initially separated by 
$4 \eta _{ {{kol}},\perp }$, so that this diagnostic can be compared with figure 13 of Yeung & Borgas (Reference Yeung and Borgas2004). We also include the results from our simulation 3H in this figure to provide a direct comparison with the isotropic hydrodynamic case. We find that anisotropic MHD turbulence develops a larger skewness of the particle-pair separations than isotropic hydrodynamic turbulence, even for similar Reynolds number simulations.
Skewness of the particle-pair separations for (a) perpendicular pairs with separation measured in the perpendicular direction, and (b) parallel pairs with separation measured in the aligned direction. For comparison with isotropic hydrodynamic turbulence, we provide the blue line from simulation 3H. The initial separation of the particle pairs is 
$4 \eta _{ {{kol}},\perp }$ in each simulation.
At early times, the values of the skewness of particle-pair separations are small and negative, but quickly become positive. The skewness rises more rapidly at earlier times for perpendicular pairs than for parallel pairs. This difference may be due to the presence of Alfvénic fluctuations, oscillations on large-scale temporal and spatial scales, which become elongated in the direction of the mean magnetic field. The initialization of particle pairs allows some pairs to reside entirely inside such flow structures. Other particle pairs saddle the boundary of a flow structure, with one particle inside and another particle outside; those pairs tend to separate more rapidly than the average. Because of the elongation of flow structures, pairs that saddle such boundaries are more likely to have an initial separation in the perpendicular direction.
As the pairs of particles move further apart in space, they experience a decorrelation between their local velocities. Eventually the rate of separation of these early fast-separating pairs slows, and the skewness reaches a peak. The peaks in skewness for our anisotropic MHD simulations occur at a time of approximately 
$10 \tau _{\eta }$, a later time than for our isotropic hydrodynamic turbulence simulation. This peak is positioned near the point of transition between the ballistic regime and the inertial subrange. We observe no clear difference in the placement of this peak with Reynolds number, or between the anisotropic directions. The peak in skewness is higher for a larger Reynolds number, and is also higher for parallel pairs than for perpendicular pairs. The difference in the height of the peak indicates that intermittency is more intense along the direction of the mean magnetic field, as well as more intense for a higher Reynolds number. This difference between the parallel and perpendicular directions suggests that the current sheets that define anisotropic MHD turbulence could be responsible for more frequent large particle separations in this setting.
5.2. Two-particle velocity statistics
To expose further differences between homogeneous isotropic hydrodynamic turbulence and MHD turbulence, we examine the separation speed, i.e. the projection of the velocity difference experienced by a pair of Lagrangian tracer particles onto the line connecting the particles. The separation speed has been used to construct stochastic models and other Lagrangian statistics (e.g. Sokolov Reference Sokolov1999; Boffetta & Sokolov Reference Boffetta and Sokolov2002), and to examine intermittency and alignment in isotropic hydrodynamic flows (Biferale et al. Reference Biferale, Boffetta, Celani, Devenish, Lanotte and Toschi2005; Yeung & Borgas Reference Yeung and Borgas2004). From early times until the end of the inertial subrange, the separation speed is significantly higher for perpendicular pairs than for parallel pairs. In our simulations, the average separation speed displays a characteristic dip, first identified by Müller & Busse (Reference Müller and Busse2007), near the beginning of the inertial subrange of time scales for MHD turbulence (see figure 8). The isotropic hydrodynamic simulation 3H is provided on the plot for comparison; it does not experience a similar dip. This dip occurs between approximately 
$2 \tau _{\eta }$ and 
$10 \tau _{\eta }$, a period where pairs of particles have separated sufficiently to sense temporal fluctuations of the velocity field. At the same time, the slope of the average dispersion curve is changing between the ballistic regime and the inertial subrange; the skewness of the pair separation distance is increasing to a peak. We find that this characteristic dip in the separation speed is deeper and lasts longer for higher Reynolds number simulations.
Average separation speed for (a) perpendicular pairs, and (b) parallel pairs. Time is given in units of the Kolmogorov time scale 
$\tau _{\eta }$. These pairs of particles are initially separated by 
$2 \eta _{ {{kol}},\perp }$.
Following the dip, the separation speed rises again as the particle pairs probe ever larger and more energetic fluctuations. This comes to an end as the pairs reach separation distances comparable to the size of the largest fluctuations in the system. We find that the separation speed ultimately reaches a higher value for higher Reynolds number simulations.
As an increasing number of particle pairs reaches the diffusion regime, the average separation velocity begins a period of fluctuating decay. At this point the large-scale velocities experienced by the separating particles tend to become decorrelated, resulting in slower, asymptotically diffusive separation dynamics. The separation speed is expected to level off at a quasi-stationary diffusive separation velocity once most pairs have left the superdiffusive region of pair dispersion. The fluctuations during this decay period are more noisy for anisotropic MHD turbulence, where Alfvénic fluctuations are present at all scales, than for isotropic hydrodynamic turbulence.
We calculate the log derivative to determine whether a scaling exists with Reynolds number for the separation speed (see figure 9). The log derivative of the separation speed shows that the dip deepens, and is not simply shifted for higher Reynolds number. Instead, a higher Reynolds number causes a milder upward slope at the end of the ballistic regime, and this persists throughout the inertial subrange. For example, at 
$10 \tau _{\eta }$ the slope is clearly different for different Reynolds number simulations. This is a consequence of the longer duration of the alignment process for the separation velocity that occurs at higher Reynolds number, as well as the stronger slowing down resulting from it.
Derivative of the log of the average separation speed for (a) perpendicular pairs, and (b) parallel pairs. These curves correspond to figure 8.
5.3. Alignment statistics
Yeung & Borgas (Reference Yeung and Borgas2004) explain the behaviour of the average separation speed in isotropic hydrodynamic turbulence through an examination of the alignment statistics using the angle between the relative velocity and separation vector of pairs of particles. Concepts of alignment also prove useful in explaining differences between isotropic hydrodynamic and MHD turbulence, an idea first explored by Müller & Busse (Reference Müller and Busse2007). Here we expand on those ideas of alignment statistics; we extend the arguments of Müller & Busse (Reference Müller and Busse2007) to examine the relationship between alignment and Reynolds number in the distinct case of 
$B_0=1$ anisotropic MHD turbulence.
To quantify how the separation of pairs of particles differs in the anisotropic MHD case, we examine the average of the cosine of the angle between the relative velocity of pairs of Lagrangian tracer particles and their separation vector. Following Yeung (Reference Yeung1994) and Yeung & Borgas (Reference Yeung and Borgas2004), we describe this angle as the ‘alignment angle’, and designate it with 
$\beta$. Examining the cosine of the alignment angle is particularly useful because it should be larger when the relative velocity and separation are better aligned. Recently Malik & Hussain (Reference Malik and Hussain2021) have defined a pair diffusion coefficient 
$K(t)$ to be the average value of the scalar product of the separation vector and the relative velocity, a quantity strongly related to the alignment angle. Their work demonstrates that the alignment angle is integral to scaling laws that can be constructed for the inertial subrange.
The average cosine of the alignment angle is clearly different for particles initially separated in the directions parallel and perpendicular to the mean magnetic field. These diagnostics also display a trend with Reynolds number. In figure 10 the average cosine of the alignment angle grows from approximately zero to a peak at approximately 
$2 \tau _{\eta }$. This measure reveals that during these early times, the separation vector of a particle pair tends toward a state in which it is better aligned with the separation velocity. The growth phase takes place during the ballistic regime of dispersion, as the particles follow the initial velocity of the fluid. For perpendicular pairs, the amplitude of the peak is approximately 0.4, which is slightly lower, but comparable to, the hydrodynamic case; for parallel pairs, the amplitude of the peak is approximately 0.28. The effect of the particle separations aligning with the relative velocity is larger in the perpendicular direction. This can be related to the higher separation speed that we observed for perpendicular pairs.
Average cosine of the alignment angle 
$\beta$. This average is shown for (a) perpendicular pairs, and (b) parallel pairs. These pairs of particles are initially separated by 2 
$\eta _{ {{kol}},\perp }$. Time is given in units of the Kolmogorov time scale 
$\tau _{\eta }$.
As particle pairs separate further in the flow, they begin to probe the lower end of the inertial subrange of time scales. The signature of this in the average cosine of the alignment angle is a drop-off between 
$2 \tau _{\eta }$ and 
$10 \tau _{\eta }$, after which the average cosine enters a plateau and exhibits noisy behaviour. The time for the first drop-off, which signifies the loss of initial strong alignment of the particle pairs, correlates with the slow down in separation velocity noted in figure 8. For a higher Reynolds number, the drop-off is larger, suggesting a link between the loss of alignment and the magnitude of the slow down in separation speed. The isotropic hydrodynamic case is provided in figure 8 for comparison; in this case the alignment experiences a small and brief dip, then regains and maintains its initial high alignment throughout the inertial subrange of time scales. For particle pairs in anisotropic MHD, regardless of their initial separation direction, the value of the average cosine of the alignment angle during the plateau is between 0.2 and 0.3; for isotropic hydrodynamic turbulence, the magnitude is nearly twice as large. Finally, at long times corresponding to the diffusive range, the average cosine of the alignment angle again decreases. For anisotropic MHD turbulence, it ultimately drops to approximately 0.1; for isotropic hydrodynamic turbulence, this value is higher, approximately 0.2. In both the inertial subrange and the diffusive regime, the relative dispersion is slower for MHD turbulence, when compared with hydrodynamic turbulence, because of the weaker alignment between the relative velocity and the separation vector.
For MHD turbulence, Müller & Busse (Reference Müller and Busse2007) demonstrated that the angle between the separation vector and the local mean magnetic field needs to be considered. We call this the ‘magnetic alignment angle’ and designate it with 
$\gamma$. The local mean magnetic field is simply the average of the total magnetic field experienced by the pair of Lagrangian tracer particles at their positions at each point in time, which includes contributions of both the fluctuating field and mean field. The average values of 
$\cos {\gamma }$ should approach zero for a perfectly isotropic turbulent flow. In an anisotropic flow, where the initial magnetic alignment angle needs to be taken into account, this is not the case. For perpendicular pairs, alignments perpendicular to the local magnetic field have an increased probability because the local magnetic field tends to be dominated by the mean magnetic field. Thus, the initial value for 
$\langle \cos {\gamma } \rangle$ is close to zero. For parallel pairs, the initial value for 
$\langle \cos {\gamma } \rangle$ is finite because the separation vector between particles consistently points in the direction antiparallel to the mean magnetic field. Regardless of the initial separation direction, as the pairs of particles evolve in time, 
$\langle \cos {\gamma } \rangle$ approaches zero because the dependency on the initial alignment state disappears. The average thus does not display a clear Reynolds number dependence in our simulations.
We find that higher-order statistics, such as the standard deviation 
$\sigma [\cos {\gamma }]$ do exhibit clear trends for our simulations (see figure 11). If no magnetic alignment angles are preferred then the distribution of 
$\cos {\gamma }$ is uniform on 
$[-1,1]$ and its standard deviation is 
$\sqrt {1/3}$, a point that is marked by a horizontal grey line in the figure. For perpendicular pairs, at early times the distribution of magnetic alignment angles is concentrated around 
$\cos {\gamma }=0$ to a higher degree than for an isotropic case. A rise in 
$\sigma [\cos {\gamma }]$ proceeds between roughly 
$2 \tau _{\eta }$ and 
$20 \tau _{\eta }$, suggesting that as a pair of particles leaves the ballistic regime, the particles achieve a wider range of alignments with the local mean magnetic field, which also results in a lower alignment of pair separation with the relative velocity vector. A more uniform distribution develops in this transient phase; the standard deviation passes through the point 
$\sqrt {1/3}$ at 
$3 \tau _{\eta }$ and continues to rise. Near the beginning of the inertial subrange, the 
$\sigma [\cos {\gamma }]$ reaches a peak, at a state where the separation vector is preferentially aligned parallel or antiparallel to the local mean magnetic field. For higher Reynolds number simulations, the peak in the standard deviation of 
$\cos {\gamma }$ is higher, corresponding to particle pairs that align more thoroughly with the local mean magnetic field. Throughout the inertial subrange and diffusive regime, this alignment decreases, indicating that there is a progressive decorrelation; the value attained for long times is close, but does not reach, the reference value of 
$\sqrt {1/3}$ for an isotropic distribution. Parallel pairs are initialized so that their separation vectors are preferentially oriented to be antiparallel to the mean magnetic field. As this distribution becomes more uniform, 
$\sigma [\cos {\gamma }]$ grows. Eventually, the parallel pairs attain a mix of antiparallel and parallel alignments, and exhibit similar trends to the perpendicular pairs at longer times.
Standard deviation of the cosine of the magnetic alignment angle 
$\gamma$, for (a) perpendicular pairs and (b) parallel pairs. These pairs of particles are initially separated by 2 
$\eta _{ {{kol}},\perp }$. Time is given in units of the Kolmogorov time scale 
$\tau _{\eta }$. The grey line indicates the value for an isotropic distribution of magnetic alignment angles.
A comparison between the different alignment statistics in figures 11 and 10 shows interesting similarities. Parallel pairs reach the peak in 
$\langle \cos {\beta } \rangle$ slightly later than the perpendicular pairs, and the peak in 
$\sigma [\cos {\gamma }]$ also occurs at a slightly later time. Although there is a lower level of alignment achieved for these parallel pairs, as measured by 
$\langle \cos {\beta } \rangle$, both types of pairs reach similar values of 
$\sigma [\cos {\gamma }]$. These observations point to two competing alignment processes at work, stemming from the physical differences in the parallel and perpendicular directions. Parallel pairs are initially strongly aligned with 
$B_0$, and they do not achieve a strong alignment between the separation vector and the relative velocity. For perpendicular pairs, a strong alignment between the separation vector and relative velocity is readily developed at early times. Then an increasing alignment with the local mean magnetic field reduces the alignment between the separation vector and relative velocity. This ultimately slows the relative dispersion for both groups of particle pairs.
6. Discussion and conclusions
We have produced Lagrangian statistics for incompressible homogeneous anisotropic MHD turbulence and examined how those statistics change as a function of the Reynolds number. Several of these diagnostics were investigated in earlier studies (Yeung & Borgas Reference Yeung and Borgas2004; Yeung et al. Reference Yeung, Pope and Sawford2006; Sawford et al. Reference Sawford, Yeung and Hackl2008) of incompressible homogeneous isotropic hydrodynamic turbulence at different Reynolds numbers. From among the many diagnostics in these earlier works, we have selected those that are most significant to the physical complexities of anisotropic MHD turbulence. Additional diagnostics have also been evaluated to further clarify the Reynolds number dependence in this new physical setting.
Anisotropic MHD turbulence is distinct from isotropic hydrodynamic turbulence in that the presence of a magnetic field interacts with the velocity of the flow, and a mean magnetic field introduces a global anisotropy in the dynamics. That anisotropy necessitates the use of a more restrictive resolution criterion, which we have introduced to assure that the energy spectra in the direction perpendicular to the magnetic field are resolved at the same level as the spectra in the direction aligned with the magnetic field. Our examination of single-particle diffusion confirms that ballistic and diffusive scaling regimes apply in each of these physically distinct directions. Two-particle dispersion also produces clear ballistic and diffusive regimes. The diffusive regime exhibits a mildly superdiffusive scaling in both parallel and perpendicular directions. During the inertial subrange, a scaling emerges as we examine larger Reynolds number simulations. We quantify this scaling by examining the log derivatives of the dispersion curves. In the perpendicular direction, for particles with separation 
$2 \eta _{ {{kol}},\perp }$, this scaling is close to 3; in the direction aligned with the magnetic field, this scaling is clearly larger than 3. The skewness of the pair separation diagnostic reveals a higher level of intermittency with larger Reynolds number, and that intermittency is more pronounced in the direction aligned with the mean magnetic field.
The presence of a mean magnetic field also introduces large-scale Alfvénic fluctuations to the system, which interact with the turbulent dynamics. We observe oscillations in the average square separation of particle pairs, clearly visible in the log-derivatives, which can be attributed to Alfvénic fluctuations. Further work would be necessary to attempt to separate the effect of waves from the fundamental anisotropy.
To explore the anisotropy more deeply with a view toward theoretical modelling, we have examined the separation velocity of particle pairs. The average separation velocity of particle pairs shows a dip following the ballistic range and preceding the beginning of the inertial subrange, between 
$2 \tau _{\eta }$ and 
$10 \tau _{\eta }$. This dip deepens for a larger Reynolds number, and provides a transition to the inertial subrange that is characteristically different from that of isotropic hydrodynamic turbulence. We probe these results further by examining alignment statistics: the cosine of the alignment angle and the cosine of the magnetic alignment angle. Slightly preceding the dip, an increase in the standard deviation of the cosine of the magnetic alignment angle begins, and this increase continues throughout the dip period. The dip period also correlates with a sharp drop in the average cosine of the alignment angle. This indicates that there are two competing alignment processes at work, namely the alignment between the separation vector of a pair of particles and the relative velocity, and the alignment between the separation vector and local mean magnetic field. The balance between these alignment processes is affected by the Reynolds number and by the anisotropy of the system. These alignment statistics show clear trends with increasing Reynolds number. A detailed evaluation of the strength of the mean magnetic field on alignment processes should be evaluated in future studies.
Funding
This material is based upon work supported by the National Science Foundation under grant no. PHY-1907876. A. B. gratefully acknowledges support via a Leverhulme Trust Research Fellowship. Simulations were performed on the Konrad and Gottfried computer systems of the Norddeutsche Verbund zur Förderung des Hoch- und Höchstleistungsrechnens (HLRN) by the project bep00051 ‘Lagrangian studies of incompressible turbulence in plasmas.’ Part of this work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
Declaration of interests
The authors report no conflict of interest.
Open access
