The stochastic behavior of magnetic field lines in turbulence is explored analytically and numerically. This problem is a fundamental aspect of turbulence research but also highly relevant in the theory of energetic particles. In the current paper, previous approaches are reviewed and some simple heuristic arguments are provided helping the reader to understand the reason for the form of analytical results. The importance of the so-called Kubo number in field line random walk theory is also discussed. Furthermore, analytical results for a position-dependent field line diffusion coefficient are provided. For more realistic turbulence configurations, the field line diffusion coefficients are computed numerically. This includes quasi-slab, quasi-2D, two-component, and three-dimensional turbulence. Specific aspects of the field line random walk in each model are also discussed. Results based on a diffusion approximation are compared with numerical results obtained without employing this approximation with the aim to explore its validity and accuracy. Numerical results based on simulations for incompressible and compressible turbulence are also discussed.

## I. INTRODUCTION

A fundamental problem of modern space and astrophysics is the understanding of turbulence (see Ref. 1 for a review). Turbulence is a very complicated topic due to the fact that it is a complex non-linear dynamical phenomenon. Turbulent magnetic fields exist on the large scales of the energy containing range where energy is injected. The energy is then transferred through the inertial range into the dissipation range where the turbulence energy in converted into other energy forms such as heat. In applications, all three ranges of the von-Karman-Kolmogorov cascade can be essential and none of those ranges must be neglected in the general case. In theories of field line random walk as well as in the description of perpendicular transport of energetic particles, the large scales of the energy range are particularly important.

In the considered scenario, we assume a superposition of a constant mean magnetic field $B\u21920=B0z\u0302$ and magnetic fluctuations with $\delta B\u2192(x\u2192)=\delta Bx(x\u2192)x\u0302+\delta By(x\u2192)y\u0302+\delta Bz(x\u2192)z\u0302$. In our own solar system, for instance, the mean field, sometimes called background field, is the solar magnetic field. The fluctuations are sometimes associated with magnetohydrodynamical waves such as Alfvén or magnetosonic waves. However, there are clearly indications that the solar wind is turbulent meaning that simplified wave pictures cannot be applied (see again Ref. 1). Very often one considers the incompressible or nearly incompressible case where we have a negligible parallel component of the magnetic fluctuations $\delta Bz\u22480$. Especially in analytical treatments of field line random walk, the assumption of nearly incompressible turbulence allows for a significant simplification of fundamental equations. Furthermore, there are observational indications that the solar wind and the very local interstellar medium are indeed nearly incompressible (see Ref. 2 for compressible and incompressible turbulence observed by the *Voyager 1* space probe). However, there are also situations where one is interested in the compressible case ($\delta Bz\u22600$) where one finds that magnetic field lines behave very differently. In the case $\delta Bz>B0$, magnetic field lines can form loops and the diffusion coefficient sometimes follows what is called *percolation scaling* (see, e.g., Refs. 3–6).

As explained in detail in this review, one has great interest in the structure of magnetic field lines in magnetic turbulence. This topic is relevant in space physics (see, e.g., Refs. 7–9) as well as in the physics of magnetically confined fusion plasmas (see, e.g., Refs. 10–15). Field lines cannot simply be computed due to the fact that they are stochastic curves related to the properties of the magnetic fluctuations. Therefore, one needs to employ methods of statistical physics in order to explore the properties of the aforementioned field lines. Alternatively, magnetic field lines and their random walk can be explored via numerical tools usually referred to as *simulations*.

The understanding of the stochastic behavior of magnetic field lines is an essential ingredient in the theory of energetic particles. This is, in particular, the case for the transport of particles across the mean magnetic field (see Ref. 16 for a review). The simplest picture of perpendicular particle transport is that of the *field line random walk limit* where it is assumed that particles follow field lines while their motion along the mean field occurs with constant velocity. In this case, the resulting perpendicular particle diffusion parameter is given by $\kappa \u22a5=v\kappa FL/2$ where we have used the particle speed *v* and the field line diffusion coefficient $\kappa FL$. However, pitch-angle scattering leads to a diffusive motion of particles in the parallel direction. This effect suppresses perpendicular transport to a sub-diffusive level. The latter process is also known as compound sub-diffusion (see Ref. 17 for a comprehensive description of this type of transport). In turbulence with transverse structure, pitch-angle scattering can scatter the energetic particles off the original magnetic field line they were tied to. Because in this type of turbulence magnetic field lines separate, one obtains an increase in the complexity of the particle motion across the mean field leading to the recovery of Markovian diffusion. In this scenario, one finds for the ratio of perpendicular and parallel particle diffusion coefficients (see Ref. 18)

where $\u2113\u22a5$ is the perpendicular correlation scale of the turbulence. In the latter case, in the FLRW limit, as well as for compound sub-diffusion, the perpendicular diffusion coefficient of the energetic particles depends on the field line diffusion coefficient. Very recently, a new analytical theory for perpendicular transport was proposed (see Ref. 19). In this theory, the field line diffusion coefficient explicitly enters fundamental equations for the particle diffusion coefficient. It has to be noted that the work discussed here focuses on magnetic field lines and energetic particles in space plasmas and astrophysical plasmas. However, a lot of work has been done in the context of magnetically confined fusion plasmas (see Refs. 10–15).

In turbulence without transverse complexity, in the literature this is usually referred to as *slab turbulence*, the theory of field line random walk is exact. However, in real turbulence there should be significant transverse structure. In this case, the calculation of field line properties, such as field line diffusion coefficients, becomes a non-trivial problem. One way to tackle this problem is the application of perturbation theory. The perturbing parameter is in this case the so-called *Kubo number*. The latter number was introduced by Kubo (see Ref. 20) and is defined via

where we have used the turbulence correlation scales in parallel and perpendicular directions, the turbulent magnetic field component $\delta Bx$, and the mean field *B*_{0}. In turbulence with large Kubo numbers, one is in a completely non-linear regime where perturbation theory fails. In this case, one is required to employ a different set of approximations. In particular, the so-called *Corrsin approximation* (see Ref. 21) sometimes also called *random phase approximation*, becomes a powerful tool. In combination with a diffusion approximation, this leads to a non-linear integral equation for the field line diffusion coefficient (see, e.g., Refs. 7 and 12). Alternatively, the diffusion approximation can be dropped leading to a second-order integrodifferential equation (see, e.g., Refs. 8 and 9).

In addition to analytical approaches in order to describe the random walk of magnetic field lines, one can also develop simulations. In this case, magnetic fields are created with the help of random number generators. Like in analytical theories, one needs to specify the properties of the magnetic fluctuations via the spectrum. Those simulations do not require to use approximations such as perturbation theory or Corrsin's hypothesis and, therefore, they are quasi-exact. Of course one has to carefully choose numerical parameters to ensure that the simulations provide the correct results. Furthermore, the simulations can be computationally expensive and, in most cases, parallel computing becomes a required tool. In Fig. 1, we show magnetic field lines generated numerically for slab turbulence corresponding to turbulence without transverse structure, and field lines for two-component turbulence where transverse complexity is significant.

It is the purpose of the current paper to review existing analytical approaches for describing the random walk of magnetic field lines. Simple heuristic arguments are presented as well. Another focus of this paper is to provide a full position-dependent description of the random walk. Furthermore, for two-component turbulence, we explore the validity of the diffusion approximation and explore in detail how a position-dependent description deviates from a diffusive description. Analytical results are also compared with simulations in order to understand the accuracy of certain approximations used in analytical treatments of the random walk. Field line random walk in various three-dimensional turbulence models is also discussed. This includes incompressible as well as compressible turbulence.

The remainder of this paper is organized as follows. In Sec. II, we briefly discuss a heuristic approach for field line random walk. This is useful for an interpretation and understanding of analytical results obtained later in this paper. In Sec. III, we discuss analytical descriptions and focus on the position dependence of different quantities. In Sec. IV, we outline numerical simulations of the field line random walk, and in Sec. V, we summarize and conclude.

## II. BASIC EQUATIONS AND HEURISTIC ARGUMENTS

In the current article, we study the behavior of stochastic magnetic field lines in turbulence. This problem is commonly referred to as the field line random walk (FLRW). In the current section, we focus on the basic physics of FLRW, the importance of the Kubo number, and present some simple heuristic arguments with the aim to find analytical forms for the field line diffusion coefficient.

### A. More on the physics of FLRW

In the literature, the term *field line diffusion* can also be found, but this is confusing. In field line random walk theory, we do not study the time-evolution of magnetic fields or field lines. Instead, we explore the probability to find a magnetic field line at a certain position of space. The statistics of field lines is usually described via the second moment $\u27e8(\Delta x)2\u27e9$ also known as *mean square displacement*. The latter quantity, in turn, can at least in some cases be expressed by a field line diffusion coefficient $\kappa FL$. The latter quantity, however, is different compared to usual diffusion coefficients. For instance, the parameter $\kappa FL$ has length dimensions $[L]$ rather than the usual $[L]2/[T]$ of particle diffusion coefficients. A related problem is the separation of two field lines in turbulence with transverse complexity. The latter problem is not studied in the current article, but it can also be important (see, e.g., Refs. 22 and 23). Both problems, field line random walk and field line separation, are depicted in Fig. 2.

In realistic scenarios of plasma physics, space science, and astrophysics, one expects that magnetic turbulence is three-dimensional. This includes so-called two-component turbulence where a superposition of slab and two-dimensional modes is considered. The latter model is commonly used to approximate solar wind turbulence where observations very clearly show a cross-distribution of magnetic fluctuations (see, e.g., Ref. 24). Although real turbulence is very complex, one can consider extreme limits to understand the basic physics of magnetic field lines. This is also useful for the interpretation of obtained results.

Normal or Markovian diffusion is defined via $\u27e8(\Delta x)2\u27e9=2\kappa FL|z|$. A generalization of the latter relation is $\u27e8(\Delta x)2\u27e9\u221d|z|\alpha $ where the special case *α* = 1 corresponds to normal diffusion, $0<\alpha <1$ to sub-diffusion, $1<\alpha <2$ to superdiffusion, and *α* = 2 corresponds to ballistic transport. Anomalous FLRW has been discussed extensively in the literature (see, e.g., Refs. 3–6, 8, 9, and 25).

### B. Field line equations and the Kubo number

As we shall see later in this paper, the starting point for any analytical or numerical investigation of field line random walk is the field line equations

where, in the general case, the turbulent magnetic field depends on all three coordinates of space $\delta Bn=\delta Bn(x,y,z)$. Therefore, the random walk is a non-trivial problem. Note in the field line equations above we have assumed $\delta Bz\u226aB0$ meaning that we focus on the nearly incompressible case one can find in the solar wind and the very local interstellar medium (see again Ref. 2). We shall come back to the compressible case ($\delta Bz\u22600$) in Sec. IV of this review. Throughout this paper, we focus on the axi-symmetric case where the symmetry axis is the mean magnetic field (the *z*-axis in our case). FLRW in non-axi-symmetric turbulence is discussed in Ref. 26.

The field line equations listed above are given in natural physical units. Alternatively, we can measure all distances along the mean magnetic field with respect to the correlation scale $\u2113\u2225$ and all distances in the perpendicular direction with respect to the scale $\u2113\u22a5$. Furthermore, we can write the magnetic field via $\delta Bx(x,y,z)=\delta Bxf(x,y,z)$. This allows us to write the first field line equation in (3) as

where we have used the Kubo number *K* as given by Eq. (2). With the help of Eq. (4), we can see in a very simple way how the Kubo number occurs naturally in field line random walk theory. We shall see throughout this paper how *K* controls the random walk. In the incompressible case, for instance, the Kubo number indicates the turnover from the quasi-linear regime ($K\u226a1$) to the non-linear regime ($K\u226b1$). For compressible turbulence, on the other hand, large values of the Kubo number can lead to *percolation scaling* (see Refs. 3–6).

### C. The initial ballistic regime

One can consider the case of very short distances where the magnetic field is still very close to its initial value. In such cases, we can approximate $\delta Bn(x,y,z)\u2248\delta Bn(0)$ in Eq. (3). In this asymptotic limit, one finds the trivial result

corresponding to the *initial ballistic regime*. A running diffusion coefficient is usually defined via

and for the normal diffusive case, we have

where $\kappa FL$ is the (constant) field line diffusion coefficient. For the initial ballistic regime, Eq. (6) becomes

corresponding to a linear increase in the diffusion coefficient with distance *z*. In any plot where a *z*-dependent field line diffusion coefficient is shown, we can see this linear increase regardless whether we show analytical or numerical results (see, for instance, Figs. 3, 5, 7, and 11 of the current paper).

### D. Heuristic arguments

First, we consider the case of turbulence with very small Kubo numbers. One can also refer to this case as quasi-slab or pure slab turbulence. We expect that in the formal limit $z\u2192\u221e$ the initial ballistic regime has passed and diffusion is restored. In a simplified description of slab turbulence, there is only one characteristic scale of turbulence that we call $\u2113\u2225$. In Fourier space, this is the so-called *bendover scale* indicating the turnover from the large scales of the energy range to the inertial range of the spectrum (see Fig. 4 for different spectra with energy range). In configuration space, on the other hand, the same scale indicates the correlation of magnetic fields (see, e.g., Ref. 27). Therefore, it is valid to either call $\u2113\u2225$ the parallel bendover scale or the parallel correlation scale. In slab turbulence, one assumes that at $z\u2248\u2113\u2225$ the turnover from the ballistic regime to the diffusive regime occurs because this is the only existing turbulence scale in the considered case. According to Eq. (8), we can then estimate $\kappa FL=\u2113\u2225\delta Bx2/B02$. The latter formula was obtained by using very simple heuristic arguments. The correct and exact analytical calculation yields

where $L\u2225$ is the *parallel integral scale*. This result is derived systematically in Sec. III of this review. However, the parallel integral scale is usually similar or at least directly proportional to the bendover scale $\u2113\u2225$. Thus, we can easily see that our very simple heuristic arguments yield the correct scaling law for the diffusion coefficient. Equation (9) is usually called the quasi-linear scaling.

We now consider the opposite case, namely, turbulence with large Kubo numbers, which can also be called quasi two-dimensional turbulence. As in the previous case, there is only one characteristic scale, namely, the bendover or correlation scale $\u2113\u22a5$. While we move along the mean field, there is an increase in $\u27e8(\Delta x)2\u27e9$. We expect that the turnover from the ballistic regime to the diffusive regime occurs as soon as $\u27e8(\Delta x)2\u27e9\u22482\u2113\u22a52$. The latter condition can be called the *transverse complexity condition* (see, e.g., Ref. 18). Note the factor 2 therein is a bit arbitrary and can be omitted. Let us assume that we need to travel a parallel distance corresponding to *L _{K}*, which is sometimes called the

*Kolmogorov–Lyapunov length*. Using this assumption in Eqs. (5) and (8) yields

and

From the first of these two relations, we can easily obtain $LK=\u2113\u22a5B0/\delta Bx$. We can clearly see that for a greater magnetic field ratio the Kolmogorov–Lyapunov length is shorter. This is because for a larger field ratio $\delta Bx/B0$ the field line spread is greater and, thus, the turnover to the diffusive regime happens earlier. Using this formula for *L _{K}* in Eq. (11) yields $\kappa FL\u2248\u2113\u22a5\delta Bx/B0$. As in the slab case, this result was obtained by employing simple heuristic arguments. Analytical theory based on Corrsin's approximation provides a slightly different result, namely,

where we have used the so-called *ultra-scale L _{U}*. It will be demonstrated in this review that there is actually a factor $2$ in front of the ultra-scale if the diffusion approximation is dropped which is possible for two-dimensional turbulence. However, Eq. (12) is often considered to be the standard result and, therefore, we use this form here. It has to be emphasized that Eq. (52) derived in the current paper is more accurate. For spectra commonly used in this field of research (see, e.g., Refs. 25 and 28), the ultra-scale is directly proportional to the perpendicular bendover scale. We call Eq. (12) the

*non-linear*or

*Bohmian*scaling.

As shown above, we can find analytical forms for the field line diffusion coefficient by using very simple heuristic arguments. Of course this approach is incomplete and cannot replace analytical and numerical treatments of field line random walk. First of all, the aforementioned formulas were obtained in the limit of very small and large Kubo numbers. Real turbulence is usually assumed to be different compared to that or at least a superposition of quasi-slab (where $K\u22480$) and quasi-2D (where $K\u2248\u221e$) modes. Furthermore, the details of the spectrum control the ratios $L\u2225/\u2113\u2225$ and $LU/\u2113\u22a5$.

## III. ANALYTICAL THEORIES

### A. The non-linear field line random walk equation

Running diffusion coefficients, mean square displacements, and *velocity auto-correlation functions* are related to each other via the *Taylor-Green-Kubo formulation* (see Refs. 29–31). This relation is given by

and the running diffusion coefficient $dFL(z)$ is related to the mean square displacement via Eq. (6). The applicability of the Taylor–Green–Kubo formula in particle transport theory was discussed in detail in Ref. 32. However, it can also be used in field line random walk theory. In this case, the *velocities* are given by $vx=dx/dz=\delta Bx/B0$. Combining Eqs. (6) and (13), therefore, yields

relating mean square displacements to magnetic field correlations along the field lines. For the magnetic field, we use the Fourier representation

To work in Fourier space is the common approach in the theory of field line random walk because it allows for the inclusion of realistic spectra of turbulence (see below for some examples). An alternative approach based entirely on magnetic correlation functions was developed in Ref. 33. With this Fourier representation, Eq. (14) can be written as

which contains a higher-order correlation function, which is difficult to handle without using further assumptions and approximations. To achieve analytical tractability, we now employ *Corrsin's approximation* (see Ref. 21)

to write this as

This means that with the help of Corrsin's hypothesis, we were able to write the correlations occurring in Eq. (16) as a product of two correlations. As shown in the following, these two correlations can be dealt with. If the turbulence is homogeneous, we have

where we have used the Dirac delta and introduced the components of the *spectral tensor*$Pnm(k\u2192)$. Combining Eqs. (18) and (19) yields

In the theory of field line random walk, *z* is the variable. For the magnetic field lines, we assume a Gaussian distribution with vanishing mean so that we can use

The assumption of a Gaussian distribution of magnetic field lines is usually tied to the assumption of normal diffusion. The reason is that a Gaussian function is the solution of a usual diffusion equation (see, e.g., Ref. 34). For non-diffusive cases, one would, therefore, expect that this assumption is no longer valid. It will be demonstrated later, however, that even for superdiffusive field line random walk the theory discussed here works very well. After using the Gaussian model in Eq. (20), we finally find the second-order differential equation for the field line mean square displacement as derived in Ref. 8, namely,

It needs to be emphasized that the only two approximations, which have been used to obtain this result, are Corrsin's approximation and the assumption of a Gaussian distribution of field lines. However, the assumption of a Gaussian statistics can easily be replaced (see, e.g., Ref. 35). Corrsin's approximation, on the other hand, can be tested by comparing results obtained by analytical theory and simulations of field line random walk (see, e.g., Refs. 26, 36, and 37). This is done throughout this review. The widespread use of the Corrsin hypothesis in diffusion theories was discussed in Ref. 38. To use this approximation in field line random walk theory is common since several decades (see, e.g., Ref. 39). We call Eq. (22) the *non-linear field line random walk equation*.

For pure slab turbulence, we have by definition $k\u22a5=0$ and the exponential in Eq. (22) vanishes. Therefore, the non-linear field line random walk equation turns into a linear equation in this very special case. This linear equation can easily be derived directly without employing Corrsin's approximation and assuming a Gaussian field line distribution. Therefore, for slab turbulence, Eq. (22) becomes an exact relation. Simplifications of Eq. (22) for a variety of turbulence models as well as a diffusion approximation are discussed in Subsections III B–III H.

### B. Slab turbulence

The slab model corresponds to a model without any transverse complexity and, therefore, it does not contain enough physics for several applications. Examples are the perpendicular transport of energetic particles where transverse structure is needed in order to restore normal diffusion (see, e.g., Refs. 16, 18, 40, and 41 for more details). Furthermore, the effect of field line separation (see, e.g., Refs. 22 and 23) is non-existent in slab turbulence. However, there is evidence that solar wind fluctuations can be approximated by a superposition of slab or slab-like perturbations and a quasi-2D component (see Secs. III D and III F for more details).

By definition, the components of the slab spectral tensor are given by

for $n,m=x,y$. All other tensor components are zero due to the solenoidal constraint $\u2207\u2192\xb7B\u2192=0$. The function $gslab(k\u2225)$ corresponds to the spectrum, which will be specified below. Slab fluctuations correspond to the extreme case where the Kubo number is exactly equal to zero. Using this spectral tensor in Eq. (22) gives us

corresponding to an exact result. There are different ways to rewrite the latter formula. One example is

where we have use the *xx*-component of the magnetic correlation tensor defined via $Rnm(z)=\u27e8\delta Bn(z)\delta Bm(0)\u27e9$. The most accurate approach to describe field line random walk in slab turbulence is to specify the spectrum $gslab(k\u2225)$ and to solve Eq. (24). However, a much simpler description can be achieved if we use the basic model

In Ref. 16, it was, for instance, shown that this is indeed correct for spectra where the inertial range spectral index is equal to 2. Using Eq. (26) in (25) gives us

Integrating this over *z* and using the running diffusion coefficient defined via Eq. (6) yields

After solving the remaining integral, we find

The latter formula provides a *z*-dependent description of the field lines. For $|z|\u226a\u2113\u2225$, one can easily obtain the initial ballistic regime as given by Eq. (8). For the case $|z|\u226b\u2113\u2225$, on the other hand, we get

corresponding to normal diffusion. One can easily see that these limits and the behavior of field line random walk in slab turbulence are in perfect agreement with the heuristic approach discussed in Sec. II of this paper.

Without assuming the form given by Eq. (26), the field line diffusion coefficient can be obtained via

to find

The integral over the Dirac delta can easily be evaluated to find

The standard model for the slab spectrum is given by (see Ref. 43)

This spectrum contains two ranges, namely, an energy range as well as a Kolmogorov (see Ref. 44) type of inertial range. The spectrum used here is perfectly flat in the energy range to ensure that the parallel integral scale is non-zero and finite (see, e.g., Refs. 16 and 28). The scale $\u2113\u2225$ denotes the bendover from the energy to the inertial range. The parameter *s* corresponds to the inertial range spectral index for which we always assume *s *=* *5/3 throughout this paper. This value corresponds to the Kolmogorov value obtained for isotropic hydrodynamic turbulence (see again Ref. 44). However, this value can also be found for magnetohydrodynamic turbulence in the solar wind (see, e.g., Refs. 1 and 45). $\delta B2$ corresponds to the total magnetic energy density so that $\delta B2=\delta Bx2+\delta By2$. The function *C*(*s*) depends only on the inertial range spectral index and follows from normalization so that

Alternatively, this can be written as

where we have used the parallel integral scale $L\u2225=2\pi C(s)\u2113\u2225$. Note this relation between the scales $\u2113\u2225$ and $L\u2225$ is based on spectrum (35). For a different spectrum, one would obtain a different relation between those two scales. In general, the parallel integral scale is defined via

For slab turbulence, this becomes (see, e.g., Ref. 16 for the detailed mathematical steps)

### C. The quasi-linear approximation

Equation (22) has to be understood as the fundamental equation in field line random walk theory. It is based on the smallest amount of approximations, namely, on Corrsin's approximation and the assumption of a Gaussian distribution of magnetic field lines. The latter assumption can easily be replaced by a different statistics. Historically, the first description of FLRW was based on quasi-linear theory (see Refs. 46 and 47) The quasi-linear limit is obtained from Eq. (22) by neglecting the exponential therein so that

If one is just interested in the diffusion limit, one can integrate the latter equation to find the quasi-linear field line diffusion coefficient

where we have used Eq. (32). However, in field line random walk theory, the quasi-linear approximation is only valid for small Kubo number turbulence. For large Kubo numbers, this approximation fails.

For small Kubo numbers, one can also develop a higher-order perturbation approach as presented in Ref. 23. In lowest-order one finds, of course, the quasi-linear scaling where the field line diffusion coefficient is directly proportional to *K*^{2}. The next-order term, we refer to it as second-order term, is then proportional to *K*^{4}. Interestingly, this term does not depend on the assumed statistics nor does it require to employ the diffusion approximation. Therefore, this next-order term is exact assuming validity of Corrsin's hypothesis. The reader can find all the mathematical details of this approach in Appendix D where the higher-order approach is combined with a *noisy slab model* for the turbulence. The corresponding field line diffusion coefficient is, in this case, given by (see again Appendix D for the derivation)

This corresponds to a reduction of the quasi-linear form. Of course the latter formula can only be used if the Kubo number is indeed small.

### D. Two-dimensional turbulence

As originally discovered by Matthaeus *et al.* (see Ref. 24), solar wind fluctuations contain a subpopulation with wave vector orientations nearly perpendicular to the mean magnetic field. This type of turbulence is usually called two-dimensional (2D) turbulence and corresponds to turbulence with a Kubo number of $K=\u221e$. As also emphasized in Ref. 24, the corresponding density fluctuations are small corresponding to the nearly incompressible case. Furthermore, two-dimensional turbulence is familiar in laboratory plasma studies such as *tokamaks* (see, e.g., Refs. 48 and 49). The theoretical framework for this type of turbulence is provided by *reduced magnetohydrodynamics* (see, for instance, Ref. 50). To understand and describe the random walk of magnetic field lines in quasi-2D turbulence is more difficult compared to the slab case because perturbation methods such as quasi-linear theory fail in this case. Therefore, one has to rely on Corrsin's approximation and equations based on this approach such as Eq. (22).

For two-dimensional turbulence, the corresponding spectral tensor has the components

where we still have $n,m=x,y$. All other tensor components are again assumed to be zero. Whereas for slab modes this is a consequence of the solenoidal constraint, the two-dimensional tensor has this property due to the fact that the considered modes are assumed to be incompressible. In principle, one could formulate a two-dimensional model for the compressible case, but this is not done in the current paper. The function $g2D(k\u22a5)$ in Eq. (44) corresponds to the two-dimensional spectrum. Possible analytical forms will be discussed below.

For the two-dimensional case, Eq. (22) simplifies to

In the following, we shall solve the latter equation without employing a diffusion approximation. This is possible as originally demonstrated in Ref. 51. First, we use the short notation $\sigma \u2261\u27e8(\Delta x)2\u27e9$ and multiply Eq. (45) by $\sigma \u2032$ to find

Very easily, we can rewrite this via

This can be integrated in order to find

One can easily see that in the formal limit $\sigma \u2192\u221e$, diffusion is restored as long as the $k\u22a5$-integral is convergent. For the field line diffusion coefficient, we then find

The *ultra-scale L _{U}* is defined via

Therewith, our field line diffusion coefficient can be written as

or

We would like to emphasize that the latter equation was derived from Eq. (45) without further assumptions or approximations. In particular, we did not employ a diffusion approximation (Subsection III E).

In order to determine the ultra-scale, we need to specify the two-dimensional spectrum. A good choice is the spectrum proposed in Ref. 25, namely,

Similar compared to our previous discussions, we have used the perpendicular bendover scale $\u2113\u22a5$ and the inertial range spectral index *s*. Compared to the slab spectrum used above, our 2D spectrum contains the additional parameter *q* controlling the behavior of the spectrum at the small wave numbers of the energy range. We refer to this parameter as the *energy range spectral index*. In Fig. 4, we have visualized the spectrum given by Eq. (53) for different values of *q*. However, the standard value used in the current paper is *q *=* *3. This choice of *q* is a bit arbitrary, but we want to ensure that fundamental length scales of turbulence such as integral scales and the ultra-scale are finite and non-zero. This forbids certain values such as *q *=* *0 or *q *=* *1 (see Refs. 16 and 28) and, therefore, *q *=* *3 is used throughout this paper. Furthermore, in Eq. (53) we have used the normalization function

where we have used again gamma functions. For this form of the spectrum, the ultra-scale (50) is computed via

The wave number integral therein can be worked out with the help of Ref. 52, and we get

where we have assumed that *q *>* *1. In order to simplify this further, we employ the relation $\Gamma (z+1)=z\Gamma (z)$ (see, e.g., Ref. 53) and use the normalization function given by Eq. (54) to obtain

Therewith, the diffusion coefficient given by Eq. (52) can be written as

In Fig. 5, we show the numerical solution of Eq. (45) for the spectrum given by Eq. (53). For those calculations, we have set *s *=* *5/3 and *q *=* *3. Furthermore, we have also visualized the diffusion coefficient given by Eq. (58) and the complete analytical solution for a Gaussian spectrum as derived in Appendix A [see Eq. (A11)]. Typical for field line random walk in two-dimensional turbulence is the slow turnover from the ballistic to the diffusive regime. This is different compared to the slab result as for instance, visualized in Fig. 3. The reason for this is that for slab turbulence the transition from the initial ballistic regime to the normal diffusive regime occurs at $z=\u2113\u2225$. For two-dimensional turbulence, on the other hand, this transition occurs at $\u27e8(\Delta x)2\u27e9\u22482\u2113\u22a52$. More details about this behavior are provided in Sec. II D, but it also follows from the non-linear field line random walk equation [see Eq. (22)].

Figure 6 shows the field line diffusion coefficient vs the magnetic field ratio $\delta Bx2/B02$. Compared is the analytical result given by Eq. (58) and the simulations described in Sec. IV of this paper. One can easily see that there is very good agreement confirming our understanding of FLRW in this specific turbulence model. In particular, the simulations confirm the non-linear scaling $\kappa FL\u221d\delta Bx/B0$.

So far we have solved Eq. (45) for the spectrum given by Eq. (53) and for an energy range spectral index of *q *=* *3. Of course the question arises what one would obtain for other values of *q*. In particular, we can explore values for which the ultra-scale is not finite. In Fig. 7, we show results for the values *q *=* *0 and *q *=* *1 in comparison with the standard value of *q *=* *3. Very clearly, we can see that the field line random walk becomes superdiffusive for smaller values of *q* as originally discovered in Ref. 25. In Ref. 51, a detailed analytical investigation of field line random walk in two-dimensional turbulence has been presented. It was shown there that the field lines become eventually diffusive for *q *>* *1. For $q\u22641$, on the other hand, one always obtains superdiffusive random walk. Assuming the form $\u27e8(\Delta x)2\u27e9\u221d|z|\alpha $, one can determine the anomalous diffusion exponent *α* by employing analytical theory. In Ref. 51, it was demonstrated that

This means, for instance, that one obtains $\alpha =4/3$ for *q *=* *0. This corresponds to weakly superdiffusive behavior and is one of the examples shown in Fig. 7. It is debatable whether this is a physical result and whether the ultra-scale needs to be finite and non-zero or not. Again the obtained analytical results are confirmed via simulations. This is particularly interesting since this agreement is found for a case where the magnetic field line random walk is not diffusive but superdiffusive. It is usually assumed that at least the assumption of a Gaussian distribution [see Eq. (21) of the current paper] breaks down if the random walk is no longer diffusive. However, Fig. 7 shows good agreement between analytical theory and simulations for the considered cases. This also confirms the validity of Corrsin's approximation.

### E. The diffusion approximation

In Subsections III A–III D, we considered extreme cases, namely, turbulence models with reduced dimensionality. However, real turbulence is usually assumed to be three-dimensional or a superposition of slab and two-dimensional modes. In such cases, a pure analytical evaluation of Eq. (22) can be difficult without employing further approximations. One way of achieving a strong simplification of the problem is to employ a diffusion approximation. Assuming that diffusion is restored in the formal limit $z\u2192\u221e$, Eq. (22) can be integrated so that

However, the *z*-integral therein can only be evaluated if one makes assumptions about the analytical form of the mean square displacement $\u27e8(\Delta x)2\u27e9$. Within the diffusion approximation considered in the current subsection, we replace the latter quantity by

meaning that we neglect the initial ballistic regime of field line random walk. Using this approximation in Eq. (60) yields

The integral therein can now easily be solved via

and we finally find the non-linear integral equation as derived in Ref. 7, namely,

It has to be noted, however, that Kadomtsev and Pogutse (see Ref. 12) derived the following equation for the field line diffusion coefficient [see Eq. (13) of their paper]

Employing the notation used in the current paper, this becomes

Taking the real part and assuming that the components of the spectral tensor are real, yields Eq. (64). We conclude that the theories of Kadomtsev and Pogutse and Matthaeus *et al.* are equivalent. However, the equation obtained by those authors was derived in different ways.

From Eq. (64), we can deduce the slab result by considering the limiting process $k\u22a5\u21920$ and using the relation (see, e.g., Ref. 42)

Very easily, one can see that we find Eq. (33). Similar considerations can be made for the case of two-dimensional turbulence. For $k\u2225\u21920$, and after using the 2D spectral tensor as given by Eq. (44), Eq. (64) simplifies to

If we compare this result with Eq. (49), we notice a factor 2 difference. This is a consequence of the diffusion approximation used in order to derive Eq. (68). Of course Eq. (49) is the more accurate result because it does not rely on this approximation. We conclude that for two-dimensional turbulence, the diffusion approximation provides a result which is a factor $2$ too small. For cases where the field line random walk never becomes diffusive (see, e.g., Fig. 7), the diffusion approximation cannot be used. For such spectra, the corresponding ultra-scale would not be finite. In such cases, Eq. (22) needs to be solved directly as it was done before in this paper.

### F. Two-component turbulence

Over the past three decades or so, the so-called two-component model of magnetic fluctuations became a standard model for approximating solar wind turbulence. The main motivation for using this model is coming from *in situ* measurements performed with various space probes (see Ref. 24 for details) where one observes a cross-type distribution of magnetic fluctuations. The interpretation is that there are slab-like (Alfvénic) fluctuations and quasi two-dimensional modes. Both modes, and the associated field line random walk, have been discussed in previous paragraphs of this paper. Bieber *et al.* (see Ref. 54) used *Helios* magnetometer data to demonstrate that solar wind turbulence possesses a dominant (roughly $80%$ by magnetic energy) two-dimensional component and a much smaller (roughly $20%$) slab contribution. On the theory side, the suggestion of a strong two-dimensional component comes from the theory of nearly incompressible magnetohydrodynamics (see Ref. 55) Based on the aforementioned papers, more research was done in recent years confirming those conclusions and providing much more details (see, e.g., Refs. 45 and 56).

The simplest spectral tensor describing the aforementioned scenario is obtained via a superposition of the slab tensor as given by Eq. (23) and the spectral tensor of the two-dimensional modes as given by Eq. (44). The resulting two-component turbulence model is, therefore, represented by

Using this model in Eq. (64) yields

where $\kappa FL$ is the total field line diffusion coefficient based on the diffusion approximation. The latter formula can easily be compared with Eqs. (33) and (68). This comparison allows us to write the above formula as

where *κ _{slab}* is the field line diffusion coefficient associated with the slab modes and $\kappa 2D$ is the diffusion coefficient related to the two-dimensional modes. Very easily, we can solve the quadratic equation (71) via

as derived the first time in Ref. 7. Thus, the recipe to obtain the total field line diffusion coefficient is simple. We employ Eqs. (33) and (68) to compute slab and two-dimensional diffusion coefficients, respectively. After that Eq. (72) is used to find the total diffusion coefficient. However, we need to keep in mind that all results obtained so far for two-component turbulence are based on the diffusion approximation. For pure slab turbulence, the FLRW theory is exact, but for two-dimensional turbulence, the diffusion approximation leads to a result, which is off by a factor $2$ as demonstrated above. At this point, it is unclear how well Eq. (72) really works and how accurate the diffusion approximation is. In order to answer this question, we have solved Eq. (22) numerically for a wide range of parameters and compared solutions with Eq. (72). Furthermore, the question arises whether Eq. (22) is even accurate since it is based on two approximations. The first one is the assumption of a Gaussian statistics, and the second one is the application of Corrsin's hypothesis. In order to explore the validity of Eq. (22) for two-component turbulence, we have compared all our analytical results with the simulations described in Sec. IV of this review.

The results are visualized in Figs. 8–10. In each plot, we have compared simulations, analytical theory represented by Eq. (22), and the results based on the diffusion approximation [see Eq. (72)]. It has to be emphasized that the simulations provide the most accurate results since they are quasi-exact. Equation (22) is based on two approximations as explained above and, thus, it is less accurate and reliable than the simulations. Equation (72) is based on the same approximations as Eq. (22), but additionally the diffusion approximation has been used. The figures show results obtained for six different turbulence parameter sets. This was done to ensure that we test our understanding of FLRW for a wide range of parameters. We can clearly see that there is overall good agreement between simulations and the different analytical theories confirming our understanding of field line random walk in two-component turbulence.

### G. Three-dimensional turbulence

In addition to the two-component model discussed in Subsection III F, full three-dimensional spectral tensors are commonly used in the literature to approximate solar wind or interstellar turbulence.

First, we consider the anisotropic model used in Ref. 57 in order to determine the parallel mean free path of cosmic rays. In this case, it is more convenient to work with spherical coordinates, which are related to cylindrical coordinates used so far via $k\u2225=\eta k$ and $k\u22a5=1\u2212\eta 2k$ where $\u22121\u2264\eta \u2264+1$. The spectral tensor for the anisotropic model has the components

where $n,m=x,y$ and all other tensor components are zero as before. It needs to be emphasized that this corresponds to the incompressible case where $\delta Bz=0$. Although the assumption of incompressible turbulence can be justified in the solar wind as explained above, we consider the compressible case later in this paper as well. For the spectral function used in Eq. (73), we follow Ref. 57 and use

However, the inertial range modeled here cannot extend to arbitrary small wave numbers and, thus, we only use this spectrum if $k\u2265kmin$ where the smallest possible wave number can be related to the largest existing turbulence scale *L* via $kmin=1/L$. In Eq. (74), we have used an anisotropy parameter Λ as well as the normalization constant *g*_{0}. Further details of the three-dimensional anisotropic model discussed here can be found in Appendix B. It has to be emphasized that the anisotropy parameter could easily be replaced by the scale ratio $\u2113\u22a52/\u2113\u22252$. Furthermore, the formal limit $\Lambda \u21920$ corresponds to quasi-2D turbulence, Λ = 1 to isotropic turbulence, and $\Lambda \u2192\u221e$ to quasi-slab turbulence. The spectrum given by Eq. (74) is similar compared to the one used in Refs. 3–6. However, in the latter papers, like in Ref. 57, a compressible version of this model has been employed. This is different compared to what we discuss in the current paragraph where we only consider the incompressible case meaning that $\delta Bz=0$.

The solutions of Eq. (22) for this spectrum are shown in Fig. 11 for different values of the anisotropy parameter Λ. One can easily see that some of the obtained results show strange behavior. Whereas the result for quasi-2D turbulence looks reasonable, the other two curves show oscillations. Mathematically, this comes due to the cosine function in Eq. (22). Physically, the observed behavior is due to the negligence of the energy range in the used spectrum. As explained above, the large scales of turbulence are important in the theory of FLRW.

Figure 12 shows the results for incompressible isotropic turbulence (Λ = 1). Although the used spectrum does not have an energy range, we clearly find diffusive behavior of field line random walk in all considered cases. Furthermore, there is almost perfect agreement between analytical theory and simulations. In Fig. 12, we can clearly see a turnover from the quasi-linear scaling $\kappa FL\u221d\delta Bx2/B02$ to the non-linear scaling $\kappa FL\u221d\delta Bx/B0$. It needs to be emphasized that the considered turbulence model is incompressible. Compressible isotropic turbulence is considered in Sec. IV C of this review.

The spectrum discussed above has the limitation that is does not contain an energy range, which is important in FLRW theory. An alternative model, which was used before in the literature, is the Gaussian model (see, e.g., Refs. 58 and 59). In this case, the spectral tensor is given by

and the model spectrum has the form

with the normalization function

Here, we have used the energy range spectral index *q*, the total magnetic field strength in the *x*-direction $\delta Bx$, and the bendover scales in the parallel and perpendicular directions $\u2113\u2225$ and $\u2113\u22a5$, respectively. The spectrum given by Eq. (76) is unrealistic in the context of space and astrophysical plasmas since it does not contain an inertial range. The latter range, however, is important in several applications. Energetic particles, for instance, can interact resonantly with the inertial range and, therefore, they experience pitch-angle scattering leading to parallel diffusion (see, e.g., Ref. 34 for a review). Therefore, the spectrum discussed in the current paragraph cannot be used for studying such scenarios. However, we can still use it in field line random walk theory where the energy range is more relevant. The mathematical details associated with using spectrum (76) in Eq. (22) can be found in Appendix C. The solutions are visualized in Figs. 13 and 14 for Kubo numbers of *K *=* *1 and *K *=* *10 and for different values of the energy range spectral index *q*. For a Kubo number of *K *=* *0.1, no significant influence of *q* on the field line diffusion coefficient was observed.

Field line random walk in other three-dimensional turbulence models was explored in the literature. Ruffolo and Matthaeus (see Ref. 60) started to explore the random walk in so-called noisy reduced magnetohydrodynamic turbulence and a Gaussian model for spectral anisotropy has been used in Ref. 61. A popular model to approximate turbulence in the interstellar medium was proposed by Goldreich and Sridhar (see Ref. 62). The stochastic behavior of magnetic field lines in this type of turbulence was explored analytically in Ref. 63.

### H. Magnetic correlations along the field lines

In advanced analytical theories for perpendicular diffusion of energetic particles, one has to deal with magnetic correlations along field lines (see Ref. 19). According to Eq. (14), such correlations are directly proportional to the second derivative of the mean square displacement. If the field line diffusion coefficient $\kappa FL$ is known, we also know the field line mean square displacement, namely, $\u27e8(\Delta x)2\u27e9=2\kappa FL|z|$ assuming the field lines behave diffusively. However, from the latter quantity one cannot obtain magnetic field correlations since one would simply get zero after calculating the second derivative of $\u27e8(\Delta x)2\u27e9=2\kappa FL|z|$. Combining Eq. (14) with Eq. (22), on the other hand, yields

In order to obtain the desired field correlations, one can use the numerically computed mean square displacement therein or, alternatively, one can employ a diffusion approximation to write this as

Without diffusion approximation, but for two-dimensional turbulence, Eq. (78) becomes

In Fig. 15, we visualize the magnetic field correlations obtained for two-dimensional turbulence and the spectrum given by Eq. (53) together with some results derived in Appendix A for a Gaussian spectrum. Clearly, we can see the rapid decay of the correlations, which is typical for normal diffusive behavior.

## IV. SIMULATIONS OF FLRW

In Secs. II and III of this review, we discussed heuristic arguments as well as analytical theories in order to describe the random walk of magnetic field lines. An alternative is provided by a pure numerical treatment of the random walk usually referred to as *simulations*. Such simulations were mainly developed in order to study the propagation of cosmic rays and solar energetic particles (see, e.g., Refs. 40, 41, and 64). However, the same tools can be used to compute magnetic field lines (see, e.g., Refs. 3–6, 26, 36, and 37). In the following, we present a detailed discussion of the simulations including the used methods, possible challenges, and, of course, some results.

### A. Magnetic field generation for two-component turbulence

Like in analytical theories, one needs to specify the spectral tensor including the different ranges of the spectrum. In simulations, the fluctuating magnetic field is usually generated by employing the formula

where we have used the magnetic fluctuations $\delta B\u2192$ as a vector. Whereas in analytical theories one needs to solve wave number integrals [see, e.g., Eq. (22)], in the simulations those are replaced by a sum over *N* wave modes. The function $A(kn)$ corresponds to the spectrum used in analytical theories and will be discussed below. The unit vectors $\xi \u0302n$ describe the orientation (polarization) of the magnetic field. For two-component turbulence, where $\delta Bz\u22480$, we usually generate those unit vectors via

In order to ensure that our system is chaotic (remember that the generated magnetic field has to be a stochastic field), the angles $\varphi n$ are obtained via a *random number generator*. The simulations presented in the current paper are based on a *MATLAB* code, which was specifically written for this review article. All random numbers therein have been created by using the *MATLAB* built-in *Multiplicative Lagged Fibonacci Generator*. Of course it is required that $0\u2264\varphi n<2\pi $. The phase of the plane waves *β _{n}* is generated with the same random number generator, and we have again $0\u2264\beta n<2\pi $. The more general case, allowing for compressible turbulence ($\delta Bz\u22600$), is discussed in Sec. IV C. As correctly noted in Ref. 64, the sets of the two random angles $\varphi n$ and

*β*are generated once for every

_{n}*n*but are then kept fixed. As a consequence, one generates the same magnetic field $\delta B\u2192(x\u2192)$ for the same set of coordinates $x\u2192=(x,y,z)$. One would expected that a large number of modes

*N*is needed in order to get an accurate result. However, over years it has been realized that this number is lower than what one would naturally expect. In fact, a typical number used in order to generate the magnetic field via Eq. (81) is

*N*=

*256. A further increase in the latter number does not alter the result. For all simulations shown in the current paper,*

*N*=

*512 was used. The parameter*

*δB*in Eq. (81) is the constant total magnetic field amplitude so that $\delta B2=\delta Bx2+\delta By2$. In Appendix E, the reader can find some more details concerning the magnetic field generation including normalization and a discussion of the solenoidal constraint.

What remains is the specification of the function $A(kn)$ in Eq. (81). First, it follows from the latter equation that this function needs to satisfy the constraint

as demonstrated in detail in Appendix E of the current paper [see Eq. (e7)]. For the function $A(kn)$, we assume that

The reason for this form can easily be understood. First of all the function $A(kn)$ needs to satisfy Eq. (83) so that our field is correctly normalized. Furthermore, $A2(kn)$ has to be directly proportional to the spectrum for which we usually employ

in agreement with the forms we have used in analytical theories [see, e.g., Eqs. (35) and (53) of the current paper]. The values used in the simulations for the energy range spectral index *q* are identical compared to the values used in analytical theory. This means, for instance, that we have set *q *=* *0 for slab turbulence but used *q *=* *3 for the two-dimensional case. Furthermore, $A2(kn)$ needs to be directly proportional to the spacing between two neighboring wave numbers $\Delta kn=kn+1\u2212kn$. In most simulations, a logarithmic spacing has been implemented so that

which is constant. In simulations, one usually specifies the minimal and maximal wave numbers *k*_{min} and *k*_{max} as well as the number of wave modes *N*. Therefore, the parameter *D* is determined and one can compute the next wave mode via

For the simulations shown in the current paper, $kmin=10\u22123$ and $kmax=10+3$ were used for all computations involving slab and two-dimensional modes in order to ensure that the used spectra have proper energy and inertial ranges.

Furthermore, we need to discuss the argument of the cosine function in Eq. (81). In numerical approaches, one works with dimensionless quantities. This means that all wave numbers shown in the current section are the real wave numbers multiplied by the corresponding bendover scale. This means that if we simulate magnetic fields in slab turbulence, we have

meaning that the wave numbers *k _{n}* in the simulations are the real wave numbers multiplied by the parallel bendover scale $\u2113\u2225$. Consequently, all distances along the mean field are normalized with respect to this scale. If we consider magnetic field lines in two-dimensional turbulence, on the other hand, we have

Here, we can also see that the angle $\varphi n$ occurs in the argument of the cosine function in Eq. (81), but it also controls the magnetic field orientations. It has to be like this in order to satisfy the solenoidal constraint $\u2207\u2192\xb7B\u2192=0$ (see the Appendix E).

### B. Obtaining field lines and their statistics

The magnetic field lines are obtained by solving Eq. (3) numerically. This can be done via basic integration methods such as an *Euler integrator* or higher-order *Runge–Kutta methods*. For a simple Euler integrator, for instance, the field line equations are solved via

where $\Delta z=zn+1\u2212zn$ is the constant step size. Of course much better numerical integration methods can be found in the literature (see, e.g., Ref. 65) such as methods with adaptive step size. An alternative method, which is equally easy to program, is the semi-implicit Euler method where we solve

This semi-implicit Euler method is a symplectic integrator. In the context of particle motion, this would lead to much better energy conservation. However, in the current paper we do not consider particles but magnetic field lines. Furthermore, we have used the semi-implicit method to redo some of the plots shown in the current paper (see, e.g., Fig. 16). No significant improvement compared to standard Euler integration was found.

For the initial values, one usually sets $x1=y1=0$, meaning that we assume that the field line goes through the origin of the Cartesian coordinate system. The magnetic fields at the field line position in Eq. (90) are computed via Eq. (81). As in the magnetic field generation, we work with dimensionless quantities so that one actually solves

In those two equations, we can nicely see how the Kubo number [see Eq. (2) of the current paper] occurs again naturally. The simulations shown in the current paper have been performed for 2500 Euler steps, but less steps can be used with superior integrators such as *Runge–Kutta methods*. In Fig. 16, we show a comparison between analytical results and simulations. We have done this for a slab model where the analytical theory of field line random walk is exact. Therefore, the slab case is perfectly suited to test our simulations. According to Fig. 16, the mean square displacement converges quickly if the number of Euler steps is increased. Therefore, to use a simple Euler integrator can easily be justified.

After solving Eq. (92), one obtains magnetic field lines (see Fig. 1 for two example plots). In order to compute mean square displacements and running diffusion coefficients, one has to average over an ensemble of computed field lines. This is computationally the most expensive aspect of the simulations because in order to obtain accurate and smooth results, one needs to work with thousands of field lines. Most results in the current paper were obtained by using 10^{4} field lines. However, in cases of large $\delta Bx/B0$ we have considered up to 10^{5} field lines. Computing such a huge amount of field lines can easily be accelerated by using parallel computing. For the simulations presented in the current paper, the *MATLAB Parallel Computing Toolbox* has been used which provides a straightforward tool for using multiple cores. Alternatively, one could even use GPUs of modern video cards.

In the literature, one can find two ways of obtaining the running field line diffusion coefficient via simulations. One is by using the definition given by Eq. (6), which is used in analytical work. However, in this case one needs to compute the derivative of a noisy result numerically which produces even more noise as commonly known. Therefore, in this type of numerical work, one often defines a running diffusion coefficient via

meaning that we divide by *z* instead of computing the derivative with respect to *z*. Of course the latter method has the advantage that it produces data with much less noise. However, both definitions are only identical for the normal diffusive case. This makes it more difficult to compare with analytical results. In Fig. 16, the exact analytical result is compared with the simulations where the running diffusion coefficient is computed via the standard definition (6) as well as a diffusion coefficient calculated via Eq. (93). Very clearly, we can see that in the latter case the data show almost no noise compared to the former case. Then, on the other hand, one can clearly see that the different diffusion coefficients agree with each other only in the asymptotic limit where the random walk becomes diffusive. This difference has to be kept in mind if simulations and different analytical results are compared to each other. For all other simulations shown in the current paper, we have used Eq. (6) so that we can compare analytical results and simulations. This explains the noise in the data despite using a huge number of field lines.

### C. Compressible isotropic turbulence

So far in this review, we have focused on incompressible turbulence meaning that we assumed $\delta Bz\u22480$. As explained before, this is motivated by measurements, indicating that real turbulence in the solar wind and the very local interstellar medium is indeed nearly incompressible (see again Ref. 2). Furthermore, analytical work is usually based on this assumption as well. In the following, we consider simulations and the case of compressible turbulence. There are indications that terrestrial magnetosheath turbulence is indeed compressible (see Ref. 66).

First, we can easily understand that field line random walk is very different for $\delta Bz>B0$. The main problem in the context of theoretical work is that in such cases magnetic field lines can form loops as visualized in Fig. 17. This is not so much a problem in simulations but a bigger problem in pure analytical work. A simple analytical approach based on magnetic field correlation models was developed in Ref. 67, but this approach is not as accurate as the analytical theories developed for the incompressible case.

For compressible turbulence models, there are different ways to generate the turbulent magnetic field in simulations. One way of approaching this problem is to generalize the method described in Sec. IV B. First, we replace Eq. (81) by

where now

and

In the three equations listed here, the parameters $\varphi n$ and *β _{n}* are random numbers (angles) as before. In order to create isotropic turbulence, for instance, the parameters

*α*and

_{n}*η*are random numbers as well. Note we have $0\u2264\alpha n<2\pi $ and $\u22121\u2264\eta n\u2264+1$. For three-dimensional anisotropic turbulence, the parameters

_{n}*α*and

_{n}*η*need to be adjusted accordingly (see, e.g., Ref. 16). For some plots shown in the current paper, we have considered a quasi-isotropic incompressible model. This means that the position dependence of the magnetic field is isotropic, but we set $\delta Bz=0$. In the simulations, this can easily be achieved by setting $\alpha n=0$ instead of using a random number for the latter quantity.

_{n}Another difference concerns the calculation of field line diffusion coefficients. Since magnetic field lines can now form loops as visualized in Fig. 17, using the parallel position *z* as variable is no longer possible. Instead, we use a *trace parameter s* and aim to compute the magnetic field line as a function of this parameter so that $x=x(s),\u2009y=y(s)$, and $z=z(s)$. This was for instance done to generate the field lines shown in Fig. 17. The magnetic field line equations we need to solve numerically are then given by

where $|B\u2192|$ is the absolute value of the total magnetic field vector. The field line diffusion coefficient is then computed via

instead of Eq. (6). Overall, this approach has been used in numerous papers with the aim to explore the magnetic field line statistics numerically for the compressible case. Very good examples can be found in Refs. 3–6.

In the current paragraph, we focus on compressible isotropic turbulence as an example. The used spectrum is the one given by Eq. (74) but for the special case Λ = 1. This also means that we neglect the energy range of the turbulence completely and set *k _{min}* = 1. This was done so that we can relate our numerical findings to previous work (see again Refs. 3–6) where exactly this type of spectrum was used.

In Fig. 18, we show the field line diffusion coefficients computed by using the method explained above. All results were obtained for compressible isotropic turbulence and the diffusion coefficients were calculated for different values of the magnetic field ratio $\delta B/B0$. First of all, we can see that for the considered cases, the field lines always behave diffusively. Furthermore, for a small magnetic field ratio, the diffusion coefficient agrees with the quasi-linear scaling

as one would expect from analytical theories [see, e.g., Eq. (9) of the current paper]. For a large field ratio, however, one finds a flatter scaling, which is clearly in disagreement with analytical theories where one finds for large Kubo numbers that the field line diffusion coefficient is directly proportional to the magnetic field ratio. However, analytical theories were developed for the nearly incompressible case. The results visualized in Fig. 18 were obtained for compressible turbulence where the analytical theories discussed in Sec. III cannot be applied.

We also like to emphasize that the flatter dependence on the ratio $\delta B/B0$ depicted in Fig. 18 agrees very nicely with the findings of Ref. 5 obtained for compressible isotropic turbulence (see Fig. 6 of the latter article). A scaling, which was discussed before in the literature, is the so-called *percolation scaling* for which we have

The percolation theory developed by Gruzinov *et al.* (see Ref. 68) was applied to the FLRW by Isichenko (see Ref. 69). Percolation scaling has been seen numerically (see Refs. 3–6). It is often stated in the literature (see, e.g., Refs. 16 and 70) that the percolation model applies to the limit of extremely high Kubo numbers. It is certainly true that one needs to consider turbulence with large Kubo number to find percolation scaling. However, one needs to consider the compressible case in order to find this type of transport. For incompressible turbulence, percolation scaling cannot be found regardless of what the Kubo number is. This can clearly be seen by considering Fig. 6 of the current paper where two-dimensional turbulence (where $K=\u221e$) is considered and one clearly finds the non-linear scaling $\kappa FL\u221d\delta Bx/B0$ as predicted by analytical theory. Furthermore, for incompressible isotropic turbulence, we find a turnover from the quasi-linear regime to the non-linear regime for an increasing value of $\delta Bx/B0$ (see Fig. 12 of the current paper). Figure 18, on the other hand, shows the results for compressible isotropic turbulence and one can clearly see a flatter dependence of the diffusion coefficient for large Kubo numbers, which is somewhat similar compared to percolation scaling.

## V. SUMMARY AND CONCLUSION

In this article, we have reviewed the random walk of magnetic field lines. To understand the stochastic properties of field lines is a fundamental part of turbulence research. Over the past few decades, significant progress has been achieved due to the development of non-linear tools allowing for an analytical description of the magnetic field lines. The key element for achieving this was and still is Corrsin's approximation [see Eq. (17) of the current paper]. The most famous and commonly used formulas for the field line diffusion coefficient are based on this approximation in combination with a diffusion approximation. The outcomes are the non-linear integral equations originally derived in Refs. 7 and 12. A more general formulation is provided by Eq. (22), which does not rely on the diffusion approximation. The latter approach has the advantage that a full position-dependent description of the random walk can be obtained including the initial ballistic regime. Furthermore, the aforementioned diffusion theories are sometimes off by a factor $2$. Equation (22) is, therefore, more accurate.

For certain applications, Eq. (64), which is based on the diffusion approximation, is still preferred since it allows for a simpler and computationally cheaper description of magnetic field lines. In the current review, we have presented a detailed comparison between analytical results obtained for two-component turbulence [see Eq. (72)] originally derived in Ref. 7 and numerical solutions of Eq. (22). The corresponding results are visualized via Figs. 8–10. In the worst case, the diffusion theory is off by a factor $2$ as predicted by analytical solutions of Eq. (22) for two-dimensional turbulence.

All analytical results obtained in this paper depend sensitively on the behavior of the turbulence spectrum in the energy range. The best example is two-dimensional turbulence for which the results are depicted in Fig. 7. In this case, certain values of the energy range spectral index lead to superdiffusive field line random walk. Therefore, if an analytical model of magnetic turbulence is developed, a proper form of the energy range should be included. This can, for instance, be seen in Fig. 11 where an anisotropic model without energy range has been used. Figures 13 and 14, on the other hand, show a well-behaving running diffusion coefficient. This was obtained for a very simple Gaussian spectrum without proper inertial range. Of course this type of spectrum cannot be used for other applications (e.g., pitch-angle scattering and parallel diffusion) because in the general case the inertial range must not be neglected either.

An alternative compared with analytical theories is provided by simulations. In this case, the magnetic fluctuations are generated with the help of random numbers. Of course one still needs to pick a certain turbulence model including the spectrum in such pure numerical work. After solving field line equations, one can average over thousands of magnetic field lines in order to obtain their statistics. This alternative approach does not require to use approximations such as Corrsin's hypothesis or the assumption of a Gaussian statistics. The downside is, however, that such simulations are computationally much more expensive. This is, in particular, the case if compared with analytical theories based on the diffusion approximation where sometimes very simple formulas can be obtained. Nevertheless, the simulations are, in particular, useful in order to test analytical theories. This was also done throughout the current review. According to Figs. 8–10, analytical theories and simulations agree well for the considered two-component turbulence model.

Analytical theories were developed with having the incompressible or nearly incompressible case in mind. This is motivated by the measurements of magnetic fluctuations in the solar wind and the very local interstellar medium. However, one can explore the random walk of magnetic field lines for the compressible case. Simulations can easily be performed for this case, but one needs to work with a trace parameter and redefine the field line diffusion coefficient. This is due to the fact that for $\delta Bz>B0$ magnetic field lines can form loops as visualized in Fig. 17 of the current paper. Field line behavior in compressible turbulence was explore numerically in Refs. 3–6. In the current review, the only considered compressible turbulence model is the isotropic model. We found that for a weaker turbulent field the diffusion coefficient is in the quasi-linear regime (see Fig. 18), but for a larger turbulent field, a flatter dependence on the Kubo number compared to the non-linear scaling was found. This dependence obtained for large Kubo numbers is somewhat similar compared to percolation scaling.

Of course the question arises how reliable analytical and numerical results are in the context of real physical scenarios. We have compared analytical theory with simulations for a variety of turbulence configurations and parameter values. Overall, we found very good agreement between the different results. This confirms our understanding of FLRW in general, but it also supports the application of Corrsin's approximation and the assumption of a Gaussian statistics. However, there are also some exceptions. As pointed out already above, there could be cases where the turbulence is compressible, but analytical theories were developed for the nearly incompressible case. Furthermore, there are flux tube-like structures in the solar wind where the magnetic field changes significantly. In this case, the magnetic field lines are confined to individual flux tubes. The aforementioned drastic changes of magnetic field directions introduce a source of turbulence intermittency and can affect the field line random walk as well as the transport of energetic particles (see, e.g., Ref. 71).

The knowledge of field line diffusion coefficients is an important ingredient in theories of energetic particle transport. Such particles experience transport across the mean magnetic field, and this, in turn, can be very important in the theory of diffusive shock acceleration (see, e.g., Refs. 72–79) and studies of solar modulation (see, e.g., Refs. 80–82). Perpendicular particle transport sensitively depends on the field line diffusion coefficient (see, e.g., Refs. 16–19 and 83). In the limit of high particle energies, the perpendicular diffusion coefficient of the particles is directly proportional to the field line diffusion coefficient. For low particle energies, on the other hand, perpendicular transport is first sub-diffusive because particles follow field lines, while the parallel motion is diffusive. As a consequence, one finds so-called compound sub-diffusion as well described in Ref. 17. In this regime, the transport is again controlled by the field line diffusion coefficient. For later times, however, the perpendicular spread of the particles has reached a point where the transverse complexity of the turbulence becomes important. Particles start to leave the original field lines they were tied to and normal diffusion is restored. In this case, the ratio of perpendicular and parallel diffusion coefficients is constant and directly proportional to the field line diffusion coefficient squared [see Eq. (1) of the current paper]. Therefore, understanding the random walk of magnetic field lines is essential in order to understand the propagation and acceleration of energetic particles such as cosmic rays.

## ACKNOWLEDGMENTS

Support by the Natural Sciences and Engineering Research Council (NSERC) of Canada is acknowledged. Furthermore, I would like to thank G. P. Zank for some very useful comments concerning magnetic turbulence in the solar wind and the very local interstellar medium. I am also grateful to G. Zimbardo for a very interesting email discussion about field line random walk in three-dimensional compressible turbulence.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: POSITION-DEPENDENT FLRW IN TWO-DIMENSIONAL TURBULENCE

For turbulence models other than the slab case, it is difficult to obtain a pure analytical position-dependent description of field line random walk. A strong mathematical simplification of FLRW in two-dimensional turbulence can be achieved if we use the Gaussian spectrum

With this spectrum, Eq. (45) becomes

The latter integral corresponds to a Gaussian integral having the solution

where we have used $\sigma =\u27e8(\Delta x)2\u27e9/\u2113\u22a52$. For the position along the mean field, we use from now on $z\u0303=z/\u2113\u22a5$. Therefore, $\sigma \u2033$ in Eq. (A3) corresponds to the second derivative of *σ* with respect to $z\u0303$. The ordinary differential equation (A3) can be solved by using standard methods. First, we multiply Eq. (A3) by $\sigma \u2032$ so that

This equation can be integrated once, and with the initial conditions $\sigma (0)=0$ and $\sigma \u2032(0)=0$, we derive

The used initial conditions are a consequence of the initial ballistic regime [see Eqs. (5) and (8) of the current paper]. We can rewrite Eq. (A5) as

Now we integrate this equation to deduce

In order to solve the latter integral, we employ the substitution $y=sinh2(x)$ to find

Using the relation

therein gives us

and we finally obtain

This formula can be used to plot the mean square displacement $\sigma \u2261\u27e8(\Delta x)2\u27e9/\u2113\u22a52$ vs position $z\u0303=z/\u2113\u22a5$. An example can be found in Fig. 5. The running diffusion coefficient can be computed by using Eq. (A5), which can be written as

We can simplify Eq. (A11) by considering asymptotic limits. For small *σ,* we find

corresponding to the initial free streaming regime as given by Eq. (5). For $\sigma \u2192\u221e$, on the other hand, Eq. (A11) can be approximated by

Practically, this means that the turnover from the ballistic regime to the diffusive regime occurs at $\u27e8(\Delta x)2\u27e9\u2248\u2113\u22a52$. In physical quantities, Eq. (A14) can be written as

and the diffusion coefficient is in this case

which is perfectly in agreement with the heuristic formula obtained in Sec. II of this paper. In the general case, the scale $\u2113\u22a5$ has to be replaced by the ultra-scale *L _{U}* as discussed in the main part of this paper [see, e.g., Eq. (12)]. Furthermore, one can combine Eqs. (A3) and (14) to find for the magnetic field correlations along a magnetic field line

For $\sigma \u21920$, we find

as expected. For larger distances, one can achieve a strong simplification by using a diffusion approximation of the form

to find for the correlations

The latter equation provides a very simple formula for the magnetic field correlations along a field line. Such correlations are visualized in Fig. 15.

### APPENDIX B: MATHEMATICAL DETAILS OF THE ANISOTROPIC TURBULENCE MODEL

In this part of the appendix, we deal with the mathematical details of the three-dimensional anisotropic model considered in the main part of this article. For the anisotropic model, the spectral tensor is given by Eq. (73) and the spectral function by Eq. (74) as long as $k\u2265kmin$. For all other wave numbers, the spectrum is assumed to be zero. The smallest possible wave number can be related to the largest existing turbulence scale via $kmin=1/L$. The normalization constant or function *g*_{0} can be computed via

According to Ref. 52, the remaining integral can be expressed by a *Gaussian hypergeometric function* so that we can write

From this, one can easily read off the normalization parameter *g*_{0}. The three-dimensional spectrum given by Eq. (74) can now be used in Eq. (22) to deduce

For numerical investigations, it is useful to define $\sigma :=\u27e8(\Delta x)2\u27e9/L2$ and $z\u0303:=z/L$. After employing the integral transformation $y=kL\u2261k/kmin$, Eq. (B3) can be written as

where we have used

### APPENDIX C: MATHEMATICAL DETAILS OF THE GAUSSIAN TURBULENCE MODEL

A simple model for magnetic turbulence is the Gaussian spectrum given by Eq. (76). It is mostly used to achieve analytical tractability. Using this spectrum in Eq. (22) allows us to find

Interestingly, the two wave number integrals therein can be solved analytically without the need of further assumptions or approximations. With the help of Ref. 52, we obtain

In order to simplify this, we employ the normalization function given by Eq. (77) to derive

Alternatively, we can use the dimensionless quantities $z\u0303:=z/\u2113\u2225$ and $\sigma :=\u27e8\Delta x)2\u27e9/\u2113\u22a52$ to write this as

where we have used the Kubo number *K* as given by Eq. (2). Very clearly, we can see that there is only one number controlling the solution of differential equation (C4) and this is the aforementioned Kubo number.

One can also combine the turbulence model discussed here with quasi-linear theory. Using the correlation tensor (75) with spectrum (76) in Eq. (42) yields

The remaining integral can be solved, and we obtain

Finally, we can employ Eq. (77) to derive

It needs to be emphasized that this result is only valid for very small Kubo numbers. However, an interesting result found here is that the energy range spectral index *q* does not impact the field line diffusion coefficient in this limit. This is, of course, a special feature of the Gaussian model considered here.

### APPENDIX D: HIGHER-ORDER PERTURBATION THEORY

The quasi-linear approximation works well in the small Kubo number regime. However, the quasi-linear approach has to be understood as a lowest-order approximation. In the following, we compute second-order corrections by following Ref. 23. We start our investigations with Eq. (20) where we even did not use the assumption of a Gaussian distribution. Regardless of what the distribution is, in the small Kubo number limit we can Taylor-expand

so that Eq. (20) becomes

Here, we have used the lowest-order (quasi-linear) mean square displacement $\u27e8(\Delta x)2\u27e90$ as well as the next-order correction $\u27e8(\Delta x)2\u27e92$. In lowest order, we have [see Eq. (41) of the current paper]

The recipe is the following. We use the latter equation to compute the lowest-order mean square displacement $\u27e8(\Delta x)2\u27e90$. Thereafter, we employ Eq. (D2) to compute the next-order contribution. In order to use all this, we need a turbulence model corresponding to small Kubo numbers. The simplest model is the noisy slab model formulated the first time in Ref. 84. Within this model, the spectral tensor is given by

where $\u2113\u22a5$ denotes a characteristic scale describing the correlation of magnetic fields across the mean field as before. Furthermore, we have employed the *Heaviside step function*$\Theta (x)$, which is defined so that

Therefore, $\Theta (1\u2212k\u22a5\u2113\u22a5)=0$ for $k\u22a5\u2113\u22a5>1$, meaning that there is a cutoff at $k\u22a5=1/\u2113\u22a5$. For the model spectrum $gslab(k\u2225)$, one can employ the same form as for slab turbulence. Due to the step function in Eq. (D4), the wave vectors are no longer aligned perfectly parallel with respect to the mean field. However, the idea here is that the parameter $\u2113\u22a5$ is large so that the width of the distribution of wave numbers in the perpendicular direction is small. This also ensures that the Kubo number is small.

One can easily show that in lowest-order perturbation theory, one finds the same result as for slab turbulence. This means, for instance, that the running diffusion coefficient is given by Eq. (29) if *s *=* *2. The latter formula can be integrated again to find for the lowest-order mean square displacement

In order to determine the second-order mean square displacement, we use Eq. (D2), which can be written as

In the latter equation, we now employ the spectral tensor given by Eq. (D4) to find

Using Eq. (24) therein allows us to write this as

Therein, we now employ Eqs. (27) and (D6) to replace the lowest-order mean square displacements. We obtain

To get the field line diffusion coefficient, we need to integrate the latter equation and employ Eq. (6). We derive after some straightforward algebra

Using the Kubo number as given by Eq. (2) allows us to write this as

We can easily see that compared to the quasi-linear result there is a reduction of the field line diffusion coefficient and that this reduction depends only on the Kubo number. It needs to be emphasized that this type of perturbation approach only works for turbulence with small Kubo numbers such as the noisy slab model used here.

### APPENDIX E: MAGNETIC FLRW VIA SIMULATIONS IN TWO-COMPONENT TURBULENCE

In the following, we discuss some details of the computer simulations described in the main part of this paper. In the simulations, the magnetic field is generated via Eq. (81). In this part of the appendix, we focus on the simulations for two-component turbulence, which is a special case of incompressible turbulence. In this case, the polarization is given by Eq. (82) where the angles $\varphi n$ are random numbers. In the following, we demonstrate how the spectrum is normalized in the simulations, and that for the considered case, the solenoidal constraint $\u2207\u2192\xb7B\u2192=0$ is satisfied.

At the initial position $x\u2192=0$, Eq. (81) becomes

The *x*-component of the latter equation is

From this, on the other hand, it follows that

If we average over all random numbers therein, we find

The two averages yield

where we have used the Kronecker delta *δ _{nm}*. Therewith, Eq. (d4) can be simplified to

Since we have $\delta B2=\delta Bx2+\delta By2$, it follows from Eq. (e6) that we need to satisfy

The latter equation corresponds to the normalization condition used in the simulations.

For general position $x\u2192$, we have for the three components of the magnetic field vector

The divergence of the fluctuating magnetic field becomes in this case

so that we satisfy the solenoidal constraint $\u2207\u2192\xb7B\u2192=0$ as required.

### APPENDIX F: SIMULATIONS FOR THE COMPRESSIBLE ISOTROPIC CASE

In the following, we discuss the computer simulations for the compressible isotropic case. In this case, the magnetic field is generated via Eq. (94). Since we consider the isotropic case, the quantities $\varphi n$, *α _{n}*,

*η*, and

_{n}*β*are random numbers. Below we discuss again normalization as well as the solenoidal constraint.

_{n}At $x\u2192=0$, Eq. (94) becomes

If we average over all random numbers, we find

Therefore, we can easily derive

Since we have $\delta B2=\delta Bx2+\delta By2+\delta Bz2$, it follows from Eq. (F4) that the spectrum needs to satisfy

as before [see Eq. (e7) of the current paper].

For general position $x\u2192$, we have for the three components of the magnetic field

The divergence of this vector field is

which is equal to zero so that we satisfy the solenoidal constraint $\u2207\u2192\xb7B\u2192=0$ as needed.