## Abstract

The results of three-dimensional fully kinetic simulations of decaying turbulence with the amplitude of the fluctuating magnetic field comparable to that of the mean field are presented. Coherent structures in the form of localized depressions in the magnitude of the magnetic field are observed to form self-consistently in the simulations. These depressions bear considerable resemblance to the so-called magnetic holes frequently reported in spacecraft observations. The structures are pressure-balanced and tend to be aligned with the local magnetic field. In the smallest structures observed, the decrease in the magnetic field strength is compensated by an increase in the electron perpendicular pressure, such that the transverse size of these structures is comparable to the electron gyroradius inside the depression. It is suggested that the structures evolve self-consistently out of the depressions in the fluctuating magnetic field, rather than being the consequence of instability growth and saturation. This is confirmed by additional, small-scale simulations, including those with realistic mass ratio between protons and electrons.

## 1. Introduction

Localized depressions in the magnitude of the magnetic field, frequently referred to as magnetic holes (e.g. [1,2], see also review in [3]), have attracted considerable interest since their initial discovery. Despite their significant history, a number of important questions remain unanswered. For example, limited information is available from observations on the three-dimensional shape of the structures or the plasma distributions. More significantly, the physical mechanisms leading to the generation of the magnetic holes are not entirely understood. These structures are often observed in high-*β* plasmas (where *β* is the ratio of plasma thermal energy to the magnetic energy density) and are frequently associated with proton anisotropy [4–6]. This has led to their interpretation as remnants of mirror instability [4]. Alternatively, magnetic holes have been proposed to result from steepening of large-amplitude Alfvén waves [7]. Setting aside the dynamics of magnetic hole formation, a number of soliton-like solutions with properties resembling magnetic holes have been introduced [8]. The majority of the observations reported magnetic holes of scales larger than the proton gyroradius, but smaller-scale structures have also been reported [9].

In this paper, we present the results of three-dimensional fully kinetic simulations of decaying plasma turbulence where long-lived coherent structures bearing many characteristics of magnetic holes are generated. The structures form at relatively early stages of the simulations and some of them persist for several Alfvén transit times. The simulations evolve a full Vlasov–Maxwell system of equations for both electrons and ions, thus providing a self-consistent model corresponding to the most complete description of collisionless plasmas. Many of the observed structures have scales comparable to the electron gyroradius within the depression, such that the electron kinetic physics plays an important role. The fully kinetic formalism provides an accurate description of such magnetic holes and enables their internal structure to be studied in detail. Additional small-scale three-dimensional simulations, including those with realistic proton-to-electron mass ratio, demonstrate that the structures can evolve self-consistently out of localized depressions in the magnetic field associated with turbulent fluctuations.

We note that the short-wavelength range of collisionless plasma turbulence has attracted considerable interest owing to its role in the dissipation of cascading energy (see e.g. the review by Howes [10]). Several of the various dissipation mechanisms discussed in the literature are associated with coherent structures, such as current sheets. In this context, the results presented in this paper are of note because they demonstrate the possibility of formation of another class of coherent structures, distinct from the planar quasi-two-dimensional current sheets that are most commonly considered.

The paper is organized as follows. The simulation parameters and the initial conditions are introduced in §2. The main results are presented in §3, including the general properties of the decaying turbulence (§3*a*) and the structure of the magnetic holes (§3*b*). The conclusions of the paper are summarized in §4.

## 2. Fully kinetic simulations

The simulations described in the paper were performed using the general-purpose plasma simulation code VPIC [11], which solves the relativistic Vlasov–Maxwell system of equations using a particle-in-cell (PIC) algorithm [12]. Turbulence was seeded by imposing an initial perturbation on uniform magnetized plasma with density *n*_{0}. The unperturbed distribution for both electrons and ions is a Maxwellian with temperature *T*_{0}. The ratio of total plasma thermal energy to the reference magnetic energy is , where *B*_{0} is the strength of the uniform magnetic field applied in the *z*-direction, ** B**=

*B*

_{0}

*e*_{z}. The simulations were conducted in a fully periodic three-dimensional domain of size

*L*

^{3}with the resolution of 2048

^{3}cells, where

*L*≈42

*d*

_{i},

*d*

_{s}=

*c*/

*ω*

_{ps}is the inertial length for species

*s*with mass

*m*

_{s}, and

*ω*

_{ps}=(4

*πn*

_{0}

*e*

^{2}/

*m*

_{s})

^{1/2}is the corresponding plasma frequency. The length of the domain is chosen in such a way that the maximum wavelength in each direction is

*k*

_{min}

*ρ*

_{i}=0.075, where

*ρ*

_{s}=(2

*T*

_{s}/

*m*

_{s})

^{1/2}/

*Ω*

_{s}is the gyroradius and

*Ω*

_{s}=

*eB*

_{0}/(

*m*

_{s}

*c*) is the gyrofrequency of species

*s*. The ion-to-electron mass ratio is

*m*

_{i}/

*m*

_{e}=50 and the ratio

*ω*

_{pe}/

*Ω*

_{ce}=2. The average number of particles per cell per species is 150, corresponding to approximately 2.6×10

^{12}total simulation particles. The time step is Δ

*tω*

_{pe}≈0.08.

The initial perturbation is of the form . Here, *p* denotes two orthogonal polarizations chosen such that *δ**B*_{1,k}⋅*B*_{0}=0, *k*⋅*δ**B*_{p,k}=0 and *δ**B*_{1,k}⋅*δ**B*_{2,k}=0. The amplitudes |*δ**B*_{p,k}|^{2}∝*k*^{−3} are chosen to yield an overall power spectrum decaying as *k*^{−1} for a range of wavenumbers with equal power in both polarizations and the root-mean-squared amplitude where 〈…〉 is the average over the domain. The phase angles *ϕ*^{x,y,z}_{p,k} are random, but are chosen to satisfy the above-listed constraints on *δ*** B**. In addition to the perturbation of the magnetic field, a velocity perturbation is loaded at time

*t*=0 by off-setting the Maxwellian particle distributions by velocity

*δ*

**that has the same form as the magnetic perturbation. The relative phase angle between the velocity and magnetic field is equal to zero for all but the six modes corresponding to the lowest**

*U**k*in each direction (figure 1

*a*). For the latter modes, the velocity and magnetic field perturbations are randomly phased. This mixture of Alfvénic and randomly phased perturbations yields the initial value of normalized cross-helicity

*σ*

_{c}≈0.44.

These initial conditions bear some resemblance to the large-scale perturbations observed in the solar wind, but also differ from them in several important ways. First, the *k*^{−1} portion of the spectrum extends up to *k*∼*d*^{−1}_{i}, much further than in the solar wind. Second, the solar wind magnetic field is typically characterized by the condition |*B*_{0}+*δ*** B**|≈const. The initial conditions used in the simulations do not enforce a similar constraint. We note that methods for constructing such magnetic fields have been proposed [13]. The structure of the magnetic field at

*t*=0 is illustrated in figure 1

*b*–

*d*, which show profiles and spectrum of

**along a randomly, chosen one-dimensional cut through the simulation. It is immediately apparent that the imposed large perturbation creates multiple regions where |**

*B***/**

*B**B*

_{0}| is rather small. As we discuss below, these imposed depressions in the magnetic field amplitude may play a role in the formation of the magnetic holes observed in the simulations.

## 3. Results

### (a) Evolution of the initial perturbation

After the simulation is initialized, the subsequent evolution is entirely self-consistent, with no applied external drive. Because the initial perturbation is a generic strong perturbation of the equilibrium that does not correspond to a solution of the Vlasov–Maxwell system, it quickly decays in time, generating fluctuations across all scales. The evolution of the spectrum of magnetic field fluctuations is illustrated in figure 2*a*,*b*. The parallel (with respect to *B*_{0}) spectrum is defined as , where is the three-dimensional fast Fourier transform (FFT) of ** B**(

*x*,

*y*,

*z*). The perpendicular spectrum is defined by , where the sum is restricted to the modes with . The shape of the perpendicular spectra remains mostly stationary after a relatively short transient interval

*t*≈0.2

*τ*

_{A}, where . The parallel spectrum

*S*

_{∥}undergoes a more pronounced evolution, becoming somewhat steeper in time. This is not surprising given that spectra in magnetized turbulence are typically anisotropic, whereas the spectrum of the imposed initial perturbation was isotropic. The shape of the spectra is further illustrated in figure 2

*c*, which demonstrates the parallel and perpendicular spectra at time

*t*≈4

*τ*

_{A}. The spectra obtained from instantaneous magnetic field

**(**

*B**x*,

*y*,

*z*,

*t*) are characterized by an upturn at high

*k*, which is a cumulative effect of the discrete particle noise inherent in the PIC algorithm. Several techniques for reducing the effects of the noise have been proposed, including spatial filtering and time averaging. Throughout this paper, we use the fields smoothed by both methods. The time-averaged fields were obtained by performing in situ averaging over 100 steps, corresponding to a time interval of approximately 0.65 times the electron gyroperiod in the reference magnetic field

*B*

_{0}. This procedure averages out high-frequency oscillations and PIC sampling noise, but should have a minimal effect on the relatively low-frequency phenomena that are of interest here. The time-averaged fields are rather expensive to compute and are available only at later times in the simulation. In order to perform uniform analysis, some of the quantities below (e.g. kurtosis) were computed from the fields that were spatially filtered using a low-pass filter with a cut-off at

*k*=3

*d*

^{−1}

_{e}[14]. The spectra of the smoothed magnetic fields obtained using both methods are included in figure 2

*c*for comparison. It is interesting to note that filtering high-

*k*noise dramatically improves the spectrum at lower

*k*, suggesting that the dominant pollution of parallel spectra is associated with sampling PIC noise at very high perpendicular wavenumbers. Both the parallel and the perpendicular spectra are well approximated by the expression , with

*α*=2.8 for

*S*

_{∥}and

*α*≈2.1 for

*S*

_{⊥}. The respective values of

*k*

_{0}are

*k*

_{0}

*d*

_{e}≈0.3 and

*k*

_{0}

*d*

_{e}≈0.5. In general, spectra of this shape are expected in a strongly dissipative system. Given the relatively low mass ratio employed by the present simulations, it is not surprising that almost the entire available range of scales could be fitted with a single function . Solar wind observations, in general, support the idea of dissipation at electron scales, although the details remain controversial [15–19].

Figure 3 demonstrates the global energy balance in the simulations by tracking in time the changes of the magnetic energy , flow energy in each species , and the ‘random’ particle energy , where *P*_{ij,s} is the pressure tensor. We note that ions receive the largest portion of the dissipated energy. This is in contrast to what is observed in a similar simulation with perturbation of lower energy and probably indicates that the partition of the dissipated energy depends on the strength of the initial perturbation, as first observed in two-dimensional simulations [20]. In figure 3, an estimate of the corresponding numerical heating was subtracted from all of the quantities. This estimate was obtained by performing a control simulation with parameters identical to the case presented here, but of smaller size and without the initial perturbation. We note that the estimate of total energy *E*_{T}=*E*_{B}+*E*_{U,i}+*E*_{U,e}+*E*_{P,i}+*E*_{P,e} remains constant in time with a good accuracy, indicating the validity of the approach.

The turbulence that develops during the decay of the initial perturbation is characterized by the presence of coherent volume-filling current structures at multiple scales. This is illustrated in figure 4, which shows volume rendering of the current density in the simulation domain at four different times. An animation showing full time evolution is available in the electronic supplementary material. At early times, when the dissipation is relatively strong, the current density distribution is characterized by the presence of intense two-dimensional current sheets. At later times, the current density structures tend to be quasi-one-dimensional elongated strands. The intermittent nature of the magnetic field is further illustrated by the typical values of the kurtosis of magnetic field elements *χ*(*s*)=〈[** B**(

**+**

*x**s*

**)⋅**

*e***]**

*e*^{4}〉/〈[

**(**

*B***+**

*x**s*

**)⋅**

*e***]**

*e*^{2}〉

^{2}. As illustrated in figure 5, the kurtosis peaks at small scales in the manner typical of both observations and other simulations. The peak value of kurtosis is observed at early times in the simulations, but only after the fluctuations develop sufficiently. In order to compute the kurtosis, the spatially filtered magnetic field was sampled along multiple trajectories in the simulation domain oriented in the

*y*-direction. The filtering is crucial to recover correct behaviour across scales [14].

### (b) Magnetic holes

An intriguing result of the simulations described here is the existence of long-lived isolated current structures (see figure 4 and the electronic supplementary material, movie S1). These structures are associated with strong depressions in the magnetic field and bear considerable resemblance to the magnetic holes frequently reported in spacecraft observations. Figure 6 demonstrates the geometrical shape and the distribution of magnetic field depressions at the later stage of the simulation *t*/*τ*_{A}≈3.8. All of the depressions are field-aligned. Depressions with relatively small perpendicular size tend to be cylindrical, whereas the larger depressions have more complicated shapes (e.g. structures marked A and B in figure 6*a*). The profiles of various quantities along the cut passing through the two largest structures are shown in figure 7. The magnetic field does not appreciably change its direction through the structure A, whereas the angle of rotation through the structure B is *θ*≈20°. Here, is the angle of rotation of the magnetic field along the trajectory passing through a reference point *x*_{0} and parametrized by unit vector ** e**. Similar to the magnetic holes observed in space, the structures are pressure-balanced, with both ion and electron perpendicular temperatures rising to compensate the depression in

*B*

^{2}. Owing to the rapid increase in plasma

*β*, the anisotropy

*T*

_{⊥}/

*T*

_{∥}inside the structures is close to the magnetohydrodynamic (MHD) threshold of mirror instability. This is highly surprising, given that the transverse size of the smallest structures is comparable to the electron gyroradius inside the structure, so that they are clearly outside of the regime where MHD considerations should apply.

In what follows, we focus on the smaller structures, which apart from their size are distinguished by nearly symmetric shapes and relatively long lifetimes. A typical example of a structure with perpendicular size of approximately one ion inertial length is shown in figure 6*b* and figure 8. Most of the profiles show qualitative behaviour similar to those observed in larger structures, with the magnetic field rotation angle *θ*<5°. However, the decrease of the magnetic field is mostly balanced by the electron perpendicular pressure in this case, whereas the ion temperature remains constant through the structure. The total anisotropy is still near marginal stability for the mirror threshold. We note that the perpendicular size of the structure corresponds to approximately one electron gyroradius defined with the peak temperature and the minimum magnetic field. This is in contrast to the larger structures described above, where the size is several electron gyroradii, but still below the ion gyroradius.

In order to gain further insights into the structure of the magnetic holes observed in the simulations, we have conducted several additional three-dimensional simulations with smaller domain sizes, including those with realistic proton-to-electron mass ratio *m*_{i}/*m*_{e}=1836. The electron-scale depressions were also observed to form in these simulations. Here, we discuss one such simulation with the size of the spatial domain reduced by a factor of 4 in each direction compared to the simulation described in §2, but other parameters (including *m*_{i}/*m*_{e}=50) remaining the same. The much smaller size of the reduced simulation made it possible to obtain field and full particle data with a much higher temporal resolution, thus allowing detailed analysis of the distribution function inside the magnetic holes. Figure 9 demonstrates some of the salient features of the electron distribution. Figure 9*a*,*d* show the *x*–*y* and *y*–*z* plane cuts through a magnetic hole. The distribution function was measured at two locations indicated by white squares in panel (*a*). Figure 9*b*,*c* show the electron distribution function at the edge of the hole (location 1 in panel (*a*)). The distribution is very well approximated by a Maxwellian at *v*_{y}<0 (which corresponds to the negative angular direction in a cylindrical coordinate system located at the hole centre). In contrast, the distribution deviates from the Maxwellian at positive *v*_{y}, so that the azimuthal current is predominantly carried by particles with , where is the initial electron thermal velocity. In order to interpret the observed features of the distribution, consider a static axisymmetric magnetic field, which in cylindrical coordinates can be written as ** B**=

*B*

_{z}

*e*_{z}. The field can be described by a single component of the vector potential

*A*

_{θ}(

*r*), so that

*B*

_{z}=(1/

*r*) d(

*rA*

_{θ})/d

*r*. The Hamiltonian for a particle of species

*s*is 3.1where

*p*

_{r}=

*m*

_{s}

*v*

_{r},

*p*

_{z}=

*m*

_{s}

*v*

_{z}, and

*Ψ*(

*r*) is the electrostatic potential. Because

*H*,

*p*

_{z}and

*p*

_{θ}are constants of motion, the problem of finding the particle orbits reduces to solving the one-dimensional equation for radial motion in an effective potential

*U*

_{eff}, 3.2and 3.3The angular motion is then determined by . Figure 10

*a*shows the profile of the magnetic field

*B*

_{z}, the vector potential

*A*

_{θ}and a quasi-potential along the cut in the

*x*-direction originating at the centre of the hole. Because the magnetic field inside the structure is mostly in the

*z*-direction, . We also assume that the electric field is mostly electrostatic, such that

*Φ*≈

*Ψ*. Because the electric field does not appear to play an important role in the radial confinement of particles inside the hole, the latter assumption is not essential. Figure 10

*b*shows the effective potential

*U*

_{eff}for particles located at

*r*≈4

*d*

_{e}with

*v*

_{r}=0 and several values of

*v*

_{θ}=

*v*

_{y}. It is clear that the particles with

*v*

_{θ}≈2

*v*

_{0}have

*p*

_{θ}≈0, and execute relatively wide orbits (note that, ignoring the electrostatic potential, only particles with

*p*

_{θ}≈0 can reach

*r*=0). Moreover, it is easy to see that particles with

*p*

_{θ}>0 never change sign of the azimuthal velocity , which indicates that their orbits encircle the centre of the hole. We conclude that the diamagnetic current supporting the magnetic hole is carried predominantly by particles executing wide orbits inside the hole with characteristic size comparable to the transverse size of the hole (the characteristic length scale for the magnetic field gradient).

The mechanism of the longitudinal confinement of particles in the hole is considerably more difficult to analyse, because this is essentially a two-dimensional problem. The small-scale magnetic holes have a mirror-like configuration of the magnetic field. If the magnetic moment was conserved for all electrons, the distribution function at mid-plane would be expected to show characteristic loss cones defined by . However, no obvious loss cones are observed in the distribution functions. This is not entirely surprising, given that a significant fraction of the electrons move on wide orbits of complicated shapes such that *μ*_{e} is not conserved. We leave detailed analysis of the longitudinal structure of the magnetic holes observed in our simulation to a future exploration, but point out that the ion distribution function *f*(*v*_{z}) shows a distinct two-beam structure throughout the hole (figure 11). The peaks appear to be related to the presence of an axial electric field *E*_{∥}. In principle, such a distribution can become unstable if the peaks are sufficiently separated.

## 4. Discussion and conclusions

In this paper, we reported on the formation of magnetic holes in fully kinetic three-dimensional simulations of decaying turbulence. The main findings are as follows:

— The simulation demonstrates that long-lived coherent structures associated with strong local depressions of the magnetic field can evolve self-consistently in a turbulent environment.

— The structures have many characteristics similar to the magnetic holes reported in spacecraft observations. They are pressure-balanced, field-aligned significant depressions of the magnitude of the magnetic field.

— The smallest structures observed in the simulations are distinguished by nearly cylindrical shapes and relatively long lifetimes. They correspond to electron-scale toroidal current sheets where the decrease of the magnetic field is supported predominantly by the increase in the electron perpendicular pressure. The transverse size of these structures is comparable to the size of electron orbits inside the structure.

— Owing to the increase in perpendicular pressure, the central regions of the depressions are slightly anisotropic.

— The ion distribution function

*f*(*v*_{∥}) inside the small-scale structures is found to have a distinct two-peak shape resulting from the presence of axial electric fields.— The magnetic holes in the simulation evolve out of slightly anisotropic ambient plasma with

*T*_{∥}>*T*_{⊥}for both electrons and ions, suggesting that mirror instability is not the origin of the observed structures. Instead, the structures evolve self-consistently out of magnetic field depressions associated with the strongly fluctuating magnetic field.

We note that the simulations used in this work solve a Vlasov equation for each plasma species coupled to a full set of Maxwell's equations, thus providing a self-consistent model with a physically accurate description of various kinetic effects. At the same time, the large-scale simulations necessarily use reduced parameters, such as rather small mass ratio *m*_{i}/*m*_{e}=50 and a modest value of *ω*_{pe}/*Ω*_{ce}=2. The presented simulations are close to the limits of current computing capabilities, and the full investigation of how the results change with the above parameters will only be possible with the advent of the next generation of supercomputers and/or advanced algorithms that can circumvent the limitations of the fully explicit PIC simulations presented here. Because the formation of the electron-scale structures can be clearly traced to diamagnetic effects associated with wide electron orbits inside the hole (see §3*b*), we expect that the main results concerning the existence and the structure of electron-scale magnetic holes will persist at realistic values of *m*_{i}/*m*_{e} and *ω*_{pe}/*Ω*_{ce}. Indeed, our small-scale three-dimensional simulations with realistic values of *m*_{i}/*m*_{e} demonstrated the formation of electron-scale holes with nearly identical structure. While this paper was under review, we also became aware of recent two-dimensional simulations with realistic mass ratios demonstrating the formation of small-scale magnetic holes [21] in a manner quite similar to that discussed here.

## Data accessibility

The full simulation data are available upon request.

## Author contributions

H.K. proposed the concept of the simulation; A.R. proposed the initial conditions; V.R. executed the simulation and wrote the paper; A.R., V.R. and H.K. extensively discussed all aspects of the project and of the paper.

## Conflict of interests

We have no competing interests.

## Acknowledgements

We gratefully acknowledge support from NASA grant NNX14AI63G at SSI. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. Additional simulations were performed on the Pleiades supercomputer provided by the NASA HEC program.

## Footnotes

One contribution of 11 to a theme issue ‘Dissipation and heating in solar wind turbulence’.

- Accepted February 9, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.