Scale statistics of current sheets in relativistic collisionless plasma turbulence

Abstract

We analyse distributions of the spatial scales of coherent intermittent structures – current sheets – obtained from fully kinetic, two-dimensional simulations of relativistic turbulence in a collisionless pair plasma using unsupervised machine-learning data dissection. We find that the distribution functions of sheet length $\ell$ (longest scale of the analysed structure in the direction perpendicular to the dominant guide field) and curvature $r_c$ (radius of a circle fitted to the structures) can be well-approximated by power-law distributions, indicating self-similarity of the structures. The distribution for the sheet width $w$ (shortest scale of the structure) peaks at the kinetic scales and decays exponentially at larger values. The data shows little or no correlation between $w$ and $\ell$, as expected from theoretical considerations. The typical $r_c$ depends linearly on $\ell$, which indicates that the sheets all have a similar curvature relative to their sizes. We find a weak correlation between $r_c$ and $w$. Our results can be used to inform realistic magnetohydrodynamic subgrid models for plasma turbulence in high-energy astrophysics.

1. Introduction

Over the last few decades, an increase in computational power has led to larger and larger direct numerical simulations of magnetohydrodynamic (MHD) turbulence (Biskamp 2008; Beresnyak 2019; Schekochihin 2022). While these simulations generally support the theory of Goldreich & Sridhar (1995) for self-similar fluctuations in MHD turbulence, the same simulations also show the presence of transient coherent structures – spatiotemporal fluctuations known as the intermittency (see, e.g., Vlahos & Isliker (2023), for a review). Among these coherent structures are current sheets, which are ribbon-shaped regions of high electric current density. Current sheets have also been identified in simulations of collisionless plasmas where the full kinetic equations are solved (e.g., Comisso & Sironi 2018, 2019).

Several groups have performed statistical studies of current sheets in MHD turbulence to categorize their geometric distributions (Servidio et al. 2010; Uritsky et al. 2010; Zhdankin et al. 2013, 2014, 2016, 2017 a ; Ross & Latter 2018). They measured the distribution of the width, length and thickness of the current sheets. The results revealed that the current sheets are indeed self-similar structures, and their geometric components obey power-law or exponential distributions. In addition, recent works have detected and analysed current sheets in collisionless regimes (Makwana et al. 2015; Azizabadi, Jain & Büchner 2021; Sisti et al. 2021; Bussov & Nättilä, 2021). Others have also attempted to associate the properties of current sheets with the statistical distribution of turbulent field fluctuations, including recent works on relativistic MHD (Chernoglazov, Ripperda & Philippov 2021) and kinetic turbulence (Davis, Comisso & Giannios 2024).

Here, we focus on current sheets in fully kinetic simulations of magnetically dominated (relativistic) plasma turbulence. In the magnetically dominated regime, a large reservoir of magnetic energy is available to cascade down and dissipate, leading to heating and acceleration of non-thermal particles. None of the previous studies have performed the geometric statistical analysis in this regime of turbulence, which has $v_{A} \approx c$ , where $v_{A}$ is the Alfvén velocity and $c$ is the speed of light. Relativistic kinetic turbulence is an area of active study in recent years, spurred by progress with particle-in-cell (PIC) simulations (Zhdankin et al. 2017 b ; Comisso & Sironi 2018; Nättilä & Beloborodov, 2021; Vega et al. 2022; Meringolo et al. 2023; Singh et al. 2024). It is relevant for many high-energy astrophysical systems (e.g., compact objects and transients). The motivation of this study is to perform a two-dimensional (2-D) statistical analysis of current sheets and compare/contrast with previous results in the literature.

In this article, we apply the computer vision algorithm developed by Bussov & Nättilä (2021) to 2-D fully kinetic (PIC) simulations of relativistic turbulence in pair plasmas. Notably, we measure the distributions of spatial scales of current sheets in the relativistic regime of kinetic turbulence. Section 2 describes the fundamental time scales and length scales and presents analytical expectations for the statistical scalings. In § 3, we discuss the numerical set-up of our PIC simulations, review the machine-learning segmentation algorithm and define our spatial scale measurements. In § 4, we present our measured distributions and their scaling indices. In § 5, we discuss the similarities between our scaling indices and previous studies and the implications for high-energy astrophysical plasmas.

2. Plasma length and time scales

We analyse the properties of current sheets in collisionless electron–positron pair plasmas. We consider pair plasma with an average temperature $\langle T \rangle$ (angular brackets $\langle \cdot \rangle$ indicate spatial averages in the remainder of the paper). Then, the average dimensionless temperature is given by

(2.1) \begin{align} \theta \equiv \frac {k_{\mathrm {B}} \langle T \rangle }{m_{\mathrm {e}}c^{2}} \,, \end{align}

where $k_{\mathrm {B}}$ is the Boltzmann constant, $m_{\mathrm {e}}$ is the electron mass and $c$ is the speed of light.

The plasma is magnetized by threading it with a uniform background field $\mathbf {B}_0 = B_{0} {\hat{\mathbf{z}}}$ . The dimensionless parameter quantifying the strength of the magnetic field is the magnetization, given by

(2.2) \begin{align} \sigma _{\mathrm {0}} \equiv \frac {B_{\mathrm {0}}^{2}}{4\pi n_0m_{\mathrm {e}}c^{2}} \,, \end{align}

where $n_{\mathrm {0}}$ is the total particle number density.

The microscopic time scale of the plasma is characterized by the (cold) plasma frequency,

(2.3) \begin{align} \omega _{\mathrm {p}} \equiv \bigg (\frac {4\pi n_0 e^2 }{m_{\mathrm {e}}}\bigg )^{1/2} \,, \end{align}

where $e$ is the electron charge. In relativistic plasmas, the plasma frequency also defines the microscopic length scale of the system, known as the plasma skin depth,

(2.4) \begin{align} \frac {c}{\omega _{\mathrm {p}}} = \bigg (\frac {m_{\mathrm {e}} c^2}{4\pi n_0 e^2}\bigg )^{1/2} \, . \end{align}

The (non-relativistic) gyrofrequency is

(2.5) \begin{align} \omega _B = \frac {e B_0}{m_{\mathrm {e}} c} \, . \end{align}

The related gyroradius is

(2.6) \begin{align} r_g \equiv \frac {c}{\omega _B} \frac {\langle p \rangle }{m_{\mathrm {e}} c} =\sigma _0^{-1/2} \frac {c}{\omega _{\mathrm {p}} } \frac {\langle p \rangle }{m_{\mathrm {e}} c} \,, \end{align}

where $p$ is the momentum of the particle. In the following, we use the skin depth $c/\omega _{\mathrm {p}}$ as the reference scale when measuring sizes. This choice is made because $c/\omega _{\mathrm {p}}$ is the largest of the kinetic scales when $\sigma _0 \gt 1$ , thus determining the characteristic width of structures. ¹

2.1. Current sheet distributions

The predominant intermittent structures in magnetically dominated plasma turbulence are current sheets. ² In MHD turbulence, these structures are responsible for a large fraction of the overall resistive dissipation (e.g., Zhdankin, Boldyrev & Uzdensky 2016). In collisionless plasmas, they are likewise thought to be sites of energy dissipation through either magnetic reconnection or other kinetic damping mechanisms (e.g., Wan et al. 2012; TenBarge & Howes 2013; Howes 2016; Parashar & Matthaeus 2016; Loureiro & Boldyrev 2017; Mallet, Schekochihin & Chandran 2017; Ripperda et al. 2021; Nättilä & Beloborodov 2022; Borgogno et al. 2022; Chernoglazov, Hakobyan & Philippov 2023; Comisso & Jiang 2023; Davis et al. 2024).

In domains with the full three spatial dimensions, current sheets have a ribbon-like morphology and, therefore, can be ideally characterized by three scales (ignoring curvature and finer structure). For 2-D simulations with an out-of-plane background magnetic field $\mathbf {B}_0$ , the largest scale is in the direction of $\mathbf {B}_0$ , which we will ignore. The remaining two scales are in the directions perpendicular to the background magnetic field; we use the terminology length $\ell$ and width $w$ to describe these two scales, ordered such that $\ell \gt w$ . See figure 1 for a schematic view of our general definitions. We emphasize that this differs from the terminology used in three-dimensional (3-D) studies.

Figure 1. Schematic illustration of a 3-D current sheet along a background magnetic field $\mathbf {B}_0$ . The two directions perpendicular to the magnetic field are denoted with the 2-D plane. The length $\ell$ of a sheet (blue curved segment), followed by the width of the sheet $w$ (green line), and the radius of curvature $r_c$ (orange line).

Numerical simulations of MHD turbulence exhibit current sheets that have a broad distribution of possible values of $\ell$ , spanning the turbulent inertial range; in particular, the data are well-described by a power-law distribution when $\ell$ is in the inertial range,

(2.7) \begin{align} P_\ell (\ell ) \equiv \frac {\mathrm {d} N_{{sheets}}}{\mathrm {d} \ell } \propto \ell ^{-\alpha } \,, \end{align}

where $N_{{sheets}} = \int P_\ell (\ell ) \, \mathrm {d} \ell$ is the number of sheets, and $\alpha \approx 3.3$ is measured in MHD turbulence simulations (Zhdankin et al. 2016). We note that $\alpha$ is used as a general index for any power-law distribution in the following text. The widths, on the other hand, mainly reside near the dissipation scale $s$ , with a narrowly peaked distribution, which we will approximate as

(2.8) \begin{align} P_w(w) \equiv \frac {\mathrm {d} N_{{sheets}}(w)}{\mathrm {d} w} \propto \exp \left ({-} \frac {w}{s} \right ) \,, \end{align}

where $s \sim c/\omega _{\mathrm {p}}$ . Furthermore, we measure the circular curvature of the current sheets, $r_c$ , which may be expected to be a power-law distribution,

(2.9) \begin{align} P_{r_c}(r_c) \equiv \frac {\mathrm {d} N_{{sheets}}}{\mathrm {d} r_c} \propto r_c^{-\alpha } \,, \end{align}

because the structures are distorted by turbulent eddies spanning a broad range of scales. The power-law index for $P_{r_c}(r_c)$ is non-trivial to predict from theory. In addition, we will be analysing the joint 2-D distributions,

(2.10) \begin{align} P_{w,\ell }(w,\ell ) \equiv \frac {\partial ^2 N_{{sheets}}(w, \ell )}{\partial w \, \partial \ell } \, \\[-24pt] \nonumber, \end{align}

(2.11) \begin{align} P_{w, r_c}(w, r_c) \equiv \frac {\partial ^2 N_{{sheets}}(w, r_c)}{\partial w \, \partial r_c} \, \\[6pt] \nonumber \end{align}

and

(2.12) \begin{align} P_{r_c, \ell }(r_c,\ell ) \equiv \frac {\partial ^2 N_{{sheets}}(r_c, \ell )}{\partial r_c \, \partial \ell } \,, \end{align}

which capture the correlations between $\ell$ , $w$ and $r_c$ . The one-dimensional (1-D) distributions follow naturally from the marginalization of these distributions as, for example, $P_w(w)= \int P_{w,\ell }(w,\ell ) \, \mathrm {d} \ell$ . We anticipate that since $w \sim s \sim c/\omega _{\mathrm {p}}$ , the width will be weakly correlated (if at all) with the other scales. The correlation between $r_c$ and $\ell$ is non-trivial to predict, in general, since structures of varying lengths will be distorted by different populations of turbulent eddies. However, we note that $r_c \propto \ell$ if the relative curvature of the structure is scale independent (i.e., the typical shape of a curved structure is independent of its length).

3. Methods

3.1. Numerical methods and set-up

We analyse the 2-D freely evolving (decaying) turbulence simulations presented in Nättilä & Beloborodov (2021). Here, we briefly review the simulation set-up, and initial conditions and refer the reader to the original paper for more details.

The simulations are initialized with a homogeneous pair plasma with a uniform density $n_{\pm, 0}$ and an initial temperature $\theta _{\mathrm {0}} = 0.3$ . The plasma is threaded by a uniform out-of-plane magnetic field $\mathbf {B}_0 = B_0 {\hat{\mathbf{z}}}$ . To seed the turbulence, we perturb the magnetic field in the directions perpendicular to the magnetic field, $\mathbf {B}_{\perp } = (B_{\mathrm {x}}, B_{\mathrm {y}},0)$ with sinusoidal large-scale modes with wavelength $l_0 = L/8$ (where $L$ is the size of the full simulation domain). The initial perturbation has amplitude $\mathbf {B}^{{rms}}_{\perp }/\mathbf {B}_0 = 1$ , where $\mathbf {B}^{{rms}}_{\perp } = \sqrt {\langle \mathbf {B}^2_{\perp } \rangle }$ .

The simulations use the relativistic PIC module in the Runko framework (Nättilä 2022). The PIC algorithm evolves all three components of the electromagnetic fields using a second-order finite-difference scheme, while particles evolve in time using a relativistic Boris pusher. These simulations are 2.5D; i.e., all three components of the particle momenta and fields are retained, while there is variation in only two spatial coordinates. In addition, we impose doubly periodic boundary conditions on the domain and perform eight passes of digital current filtering at each time step. The current filtering is performed immediately after the particle’s current deposition using a three-point binomial filter to smooth out the numerical noise that results from a finite number of computational particles.

The 2-D domain in our simulations is a square in the x–y plane of size $L = 1024 c/\omega _{\mathrm {p}}$ , and it is covered with $5120^{2}$ cell resolution. The (initial) plasma skin depth is resolved by five grid cells. In our 2-D simulations, we initialize fluctuations on a scale of $l_{0} \approx 125 c/\omega _{\mathrm {p}}$ . We simulate the turbulence for $20 l_{0}/c$ . We focus on two simulations with $\sigma _0 = 1$ and $10$ , which correspond to the transrelativistic (magnetically dominated) regime. During the simulation, the plasma heats up and so $\theta _{0}$ increases (see figure 6 in Nättilä & Beloborodov (2021)). For our simulation with $\sigma _0 = 1$ at $t = l_{0}/c$ , we measure $\theta \approx 0.45$ which then increases to $\theta \approx 0.60$ at $t = 15l_{0}/c$ . Similarly, for the case of $\sigma _0 = 10$ , we first measure $\theta \approx 0.50$ which then increases to $\theta \approx 2.70$ . Thus, the temperatures become only mildly relativistic, and it is reasonable to use the non-relativistic expressions for kinetic scales from § 2.

We present six visualizations of one of our simulations with $\sigma _0 = 10$ at $t = 4.6\, l_0/c$ in figure 2. The magnetic eddies are visible for the fully developed turbulence as large plasma over-densities (figure 2 a). Secondly, we can visually trace out two types of structure in the out-of-the-plane current density component $J_z$ (figure 2 b): elongated stripes (current sheets) and circular spots (plasmoids). These same structures are also visible as coherent structures in the (perpendicular) magnetic field $\mathbf {B}_\perp$ (figure 2 c): the current sheets coincide with regions of alternating magnetic fields (in between colliding plasmoids). The plasmoids coincide with the location of magnetic loops. The regions with current sheets (with antiparallel fields) are prominent locations of energy dissipation. The sheet regions can be associated with areas of energy dissipation when visualizing the work done by the electric field $\propto \mathbf {J}\cdot \mathbf {E}$ (figure 2 d). We note that a significant fraction of the regions with high electric currents have magnetic fields that are parallel (rather than antiparallel) on both sides (Ha et al. 2024).

A standard method to segment the current sheets into individual structures is to apply a threshold in the current density. To illustrate this method, we trace out the outlines of the sheets by applying a threshold of $J_{\mathrm {z}}/J_{{rms}} \gt 3$ (figure 2 e). In our case, however, we use the regions originating from a more advanced region-of-interest (ROI) detection (measured at a time $4l_0/c$ ), relying on a computer vision segmentation algorithm presented in Bussov & Nättilä (2021) (figure 2 f). Such segmentation method has the advantage that it clusters the sheets based on not only the electric current but also other features; this enables us to separate the current sheets and plasmoid cores from each other. The details of the method are given in the next section.

3.2. Overview of the machine-learning algorithm

We employ the unsupervised ensemble machine-learning algorithm developed by Bussov & Nättilä (2021) to segment the current sheets from fully kinetic plasma turbulence simulations. The algorithm is an ensemble extension of the self-organizing map (SOM) algorithm (Kohonen 2001, 2013) dubbed ‘statistically combined ensemble’ (SCE), which merges multiple (independent) SOM evaluations together, hence increasing the robustness of the ROI boundaries.

Briefly, we outline how the algorithm works. We apply the SOM algorithm to our simulation snapshots and cluster the data using three features: $\boldsymbol{B}_{\perp }$ , $\boldsymbol{J}_{\boldsymbol{z}}$ and $\boldsymbol{J}_{\boldsymbol{z}} \cdot \boldsymbol{E}$ , where $\boldsymbol{B}_{\perp}$ is the (in-plane) perpendicular magnetic field, $\boldsymbol{J}_{\boldsymbol{z}}$ , is the (out-of-the-plane) parallel component of the current and $\boldsymbol{J}_{\boldsymbol{z}} \cdot \boldsymbol{E}$ is a measure of the work done by the electric field $\boldsymbol{E}$ . Each pixel on the image then defines a multidimensional feature-space data point.

Figure 2. Visualization of the analysed turbulence data measured at time $t = 4.6 l_0/c$ : plasma density $n/n_0$ , where $n_0$ is the initial plasma density (a); out-of-the-plane current density $J_{z}/n_0 e c$ (b); strength of the in-plane magnetic field component $\sqrt {B_x^2+B_y^2}/B_0$ , and the field lines (red curves), where $B_0$ is the initial guide field strength (c); a proxy of the work done by the parallel electric field $\mathbf {J} \cdot \mathbf {E}/\sqrt {\langle (\mathbf {J} \cdot \mathbf {E})^2 \rangle }$ in units of the root mean square value (d); regions of the current density with $J/J_{{rms}} \gt 3$ , where $J_{{rms}} = \sqrt {\langle J^2 \rangle }$ (e); and, current sheet regions from the SCE algorithm (f). The SCE algorithm is shown at $t = 4.0 l_0/c$ .

The SOM algorithm is based on a 2-D network of neurons that are trained to approximate the multidimensional (feature-space) data. In practice, the map can be thought of as a (finite) surface that is fitted to the input data points. Conversely, each feature-space data point can be associated with a best-matching neuron on the map. The map defines a latent (reduced) data space for the input data.

The neural map is then used to cluster the feature-space data into different groups by associating a region on the map (i.e., a set of connected neurons) as one data cluster. The (feature-space) distance between the neurons defines a convenient cost function to separate the neurons into a finite number of clusters. In theory, the number of resulting clusters is arbitrary. Still, in practice, the map often has distinct regions of closely connected neurons that help set the number of resulting clusters from the analysis. In our case, the SOM algorithm results, on average, in four distinct clusters, one of which is the current sheets (see, Bussov & Nättilä (2021) for more details). ³

In practice, we use an ensemble version of the SOM to segment the sheets. This so-called SCE method combines multiple SOM segmentations into one clustering realization to increase the signal-to-noise ratio of each ROI boundary. In general, our machine-learning-based segmentation method should be contrasted with the previously used threshold technique, where the boundaries of the current sheets are found by setting a threshold value for the magnitude of the electric current (see figure 2 f). We emphasize that our definition of a current sheet is more complex than what results, e.g., from segmenting the sheets with the threshold method: our definition of the current sheets is based on a combination of all three input features. The second reason, in addition to the more generic definition of the current sheet cluster, is that the SCE method is computationally faster, making it feasible to construct large structure catalogues. We note, however, that the algorithm is currently calibrated only on 2-D data, so we focus on that here.

3.3. Measurements

The following section presents our definitions of geometric length and width. Specifically, we analyse a time series of 2-D masks separated by $\Delta t = 1 l_0/c$ for two simulations with $\sigma _0 = 1$ and $\sigma _0 = 10$ . We analyse a time interval, $l_0/c \leqslant t \leqslant 15 l_0/c$ when the turbulence is fully developed but still sustains large amplitude fluctuations. In practice, when analysing the structures, we use the open-source geometry library `Shapely’ (Gillies et al. 2024) to automate the geometric analysis of the structures. First, we dissect each independent closed contour in a mask into $(x,y)$ coordinate lists that trace the element’s boundary. This procedure is repeated for every mask in the time series. Specifically, we use the Shapely.Polygon class, which takes a set of $N_p$ data points and connects the chain of coordinate points into a ‘ring’. We generate an instance of this polygon class for each structure in a given image. Once we generate the Polygon class to represent the current sheet, we measure its perimeter $\mathcal {S}$ ( $=\oint \mathrm {d} s$ , where the line integral is performed along the closed curve around the polygon) and use it to define a robust estimator of the sheet length $\ell$ as

(3.1) \begin{align} \mathcal {S} \approx 2\ell + 2 w \approx 2\ell, \end{align}

so that $\ell \approx \mathcal {S}/2$ , where we assumed that $w \ll \ell$ in the second part of the equation. We also tried measuring $\ell$ by taking the pairwise distance between each coordinate point and identifying the largest distance as $\ell$ . However, we forgo this definition since we find it to be a noisier estimator of $\ell$ . We also note that when accounting for the possible curvature of the structure, the natural length is the arclength of the structure, not the distance along the structure’s longest dimension; therefore, our definition $\ell$ also naturally includes the curvature effects.

To measure the width (see appendix A for the full derivation), we take each structure and calculate the major axis of each sheet, which is defined as the largest distance between any set of points. Then, we rotate the structure such that its major axis lies on the $x$ axis. Then, we translate the structure to the origin. Once in this new reference frame, we measure the width using (B1).

Finally, we introduce a measurement of the curvature of the current sheet, dubbed ‘curvature radius’, $r_{c}$ . The radius of curvature is calculated by fitting a circle through the three points on the curve: leftmost point, centre point and rightmost point (in the new reference frame described above). Then, a circle with a radius $r_c$ can be defined to pass through each point. In the limit that the curvature radius $r_c \rightarrow 0$ , the structure is highly curved, and as $r_c\rightarrow \infty$ , the structure becomes flat. We detail our simple circle fitting algorithm in appendix C.

Figure 3. Examples of current sheets in our catalogue. Individual current sheets of varying length, $\ell$ , and radius of curvature, $r_{c}$ , are visualized in the transformed position space $(x',y')$ using (A1) and (A2). The three red points mark the locations used to calculate the curvature radius using (C2). The green line shows the arc of the circle passing through the three points, with the centre calculated using (C6).

In figure 3, we display representative examples from the resulting current sheet catalogue for the simulation with $\sigma _0 = 10$ . These examples show that most of the sheets have approximately constant width $w$ .

Lastly, the approximation of the sheet curvature as part of a circle’s arc seems appropriate in almost all cases; in total, we find that only approximately $0.5\,\%$ of the sheets can not be adequately modelled as arcs (but are, e.g., $S$ - or $L$ -shaped). Such non-standard morphologies can sometimes result from the interactions and mergers of sheets.

4. Results

Our catalogue contains $\approx \!13\,000$ individual objects for $\sigma _0 = 1$ and $\approx \!10\,000$ for $\sigma _0 = 10$ simulations. Each time snapshot contains, on average, $\approx \!1\,000$ current sheets and provide detailed time-dependent statistics of the structures. The total volume-integrated magnitude of the electric current closely traces the evolution of the decaying turbulence: for $t = 0$ to $\sim 3 \,l_0/c$ , we observe an initial transient period where the turbulent structures form and the magnitude of the electric current rapidly grows; for $t \gtrsim 3 \,l_0/c$ , the turbulence decays and the magnitude of the current decreases slowly (see also, e.g., Nättilä & Beloborodov 2021). The analysed current sheets are long-lived intermittent structures and their time evolution is analysed and discussed in more detail, e.g., in Imbrogno et al. 2024.

Figure 4. The filling fraction of current sheets, $f$ , of our simulation of $\sigma _0 = 10$ measured at $4.6l_0/c$ , as function of $J_{{thr}}/J_{{rms}}$ (blue line), where $J_{{thr}}$ is the threshold factor, and where $J_{{rms}} = \sqrt {\langle J^2 \rangle }$ is the root mean square of the current density. The star symbol denotes the filling fraction measured using our machine-learning algorithm. We also fit the MHD data from Zhdankin et al. (2016) (green line).

We calibrate our current sheet identification method with the previous method of ‘thresholding’ used in Zhdankin et al. (2014) and Zhdankin et al. (2016). In figure 4, we plot the volume filling fraction $f$ of the structures occupying pixels having current densities $|\mathbf {J}| \gt J_{{rms}}$ , where $J_{{rms}}$ is a threshold parameter that we vary. The filling fraction of the structures with current densities above the threshold exhibits a steep exponential decay at low $J_{{thr}}/J_{{rms}} \lt 3$ , which then flattens to a less steep decay for $J_{{thr}}/J_{{rms}} \gt 3$ . We also compare our results with the MHD analysis by Zhdankin et al. (2016) which shows a very similar behaviour for $J_{{thr}}/J_{{rms}} \gt 3$ . We find that the filling fraction of structures detected by the SCE algorithm agrees with the thresholding method when $J_{{thr}}/J_{{rms}} \approx 3$ (the red star in figure 4). This indicates the ‘effective’ threshold of the SCE algorithm.

Figure 5. Width distributions, ${\rm d}N/{\rm d}w$ , for simulations with initial $\sigma _0 = (1, 10)$ . Two reference exponential fits with the corresponding index are shown in the legend.

Figure 6. Length distribution, ${\rm d}N/{\rm d}\ell$ , for simulations with initial $\sigma _0 = (1, 10)$ . Three reference power laws with the corresponding index are shown in the legend.

Figure 7. Radius of curvature distributions, ${\rm d}N/{\rm d}r_{c}$ , for simulations with initial $\sigma _0 = (1, 10)$ . Two reference power laws with the corresponding index are shown in the legend. This distribution is restricted to current sheets with $\ell \gt 7.5 c/\omega _p$ , to avoid spurious measurements from small structures.

In figure 5, we plot the width distribution, $P_w(w)$ , of our current sheets for both values of $\sigma _0 = 1$ and $10$ . The distribution resembles an exponential distribution, peaking at around $3 \, c/\omega _{\mathrm {p}}$ . Qualitatively, the distributions are best modelled with an exponential function, $\exp (-w\omega _{\mathrm {p}}/c)$ . The width distribution does not greatly vary over time.

In figure 6, we visualize the length distribution of our current sheets, $P_{\ell }(\ell )$ , and similarly to the width distribution, we find that the distribution remains roughly constant over time. The distribution can be modelled as a power-law distribution with an index, $\alpha = 2.5$ . The distribution peaks at $5 c/\omega _{\mathrm {p}}$ , and extends to $\sim \!100 \, c/\omega _{\mathrm {p}}$ ; i.e., close to the initial energy-carrying scale $l_0$ of the simulation.

Next, in figure 7, we visualize the curvature distribution $P_{r_c}(r_c)$ (although we omit structures with $\ell \lt 7.5 c/\omega _{\mathrm {p}}$ because the curvature measurement becomes unreliable for very short structures). The distribution is a clear power-law with an index $\alpha = 2.0$ . The distribution peaks at $3 \, c/\omega _{\mathrm {p}}$ , and extends to approximately $1000 \, c/\omega _{\mathrm {p}}$ . The distribution is steady in time, with no significant evolution.

Furthermore, we study the 2-D joint distributions of the width and length, $P_{w, \ell }(w,\ell )$ , in figure 8. We find no clear scaling between the two quantities. The distribution is bounded such that the upper limits of the distribution are fit by $w \propto \ell ^{\alpha }$ , where $\alpha = 1.0$ . We fit $2\log _{10}(\ell )$ for the lower bound.

Figure 8. The 2-D histogram correlating the average width and length distributions for simulations with initial $\sigma _0 = 1, 10$ . The power law index fits the distributions shown in the legend.

Figure 9. The 2-D histogram correlating the radius of curvature and length distributions for simulations with initial $\sigma _0 = 1, 10$ . The power law index fits the distributions shown in the legend.

Figure 10. The 2-D histogram correlating the width and radius of curvature distributions for simulations with initial $\sigma _0 = 1, 10$ . The power law index fits the distributions shown in the legend.

Similarly, in figure 9, we visualize the 2-D distributions of the curvature and length, $P_{r_c, \ell }(r_c,\ell )$ . The mean value of the sheets linearly scales with $\ell$ (dark orange line). This suggests that the relative curvature to size is the same for all $\ell$ . Additionally, we find that the distribution can be traced with $r_{c} \propto \ell /\pi$ . We also show a fit that corresponds to $r_{c} \propto \ell ^{\alpha }/2$ , where $\alpha = 1$ . One final note is that at $\ell \lt 7 c/\omega _{\mathrm {p}}$ the curvature measurement becomes unreliable and one can see a small ‘cluster’ of points at $r_c \approx 100 c/\omega _{\mathrm {p}}$ and $r_c \approx 300c/\omega _{\mathrm {p}} - 400 c/\omega _{\mathrm {p}}$ .

Finally, in figure 10, we visualize the 2-D distribution of the width and curvature, $P_{w, r_c}(w,r_c)$ . Here, no strong correlation is seen. The upper end of the distribution is traced with $w \propto \ell ^{\alpha }$ , where $\alpha = 1$ . Since the structures have widths close to the dissipation scale, it seems reasonable that they are independent of the curvature, which is more likely sensitive to MHD-scale fluctuations.

5. Discussion

We have analysed various length statistics of tens of thousands of current sheets found in freely evolving (decaying) plasma turbulence. In particular, we have found the following.

1. (i) The current sheet length distribution has a power law tail with index $\alpha = 3.3$ ; it is shallower ( $\alpha \approx 2$ ) at small lengths ( $\ell \lesssim 30 \, c/\omega _{\mathrm {p}}$ ).

2. (ii) The width distribution is sharply peaked at $w \sim 3\, c/\omega _{\mathrm {p}}$ ; it resembles an exponential decay at higher $w$ . There is a slight logarithmic increase of $w$ as a function of $\ell$ .

3. (iii) The characteristic curvature radius is $r_c \sim \ell /2$ , hence the sheets are slightly curved. However, we find that short ( $\ell \sim 10\, c/\omega _{\mathrm {p}}$ ) sheets are often flatter (larger $r_{c}$ ); hence, the curvature distribution is a shallower distribution than the length distribution (as seen by the 1-D histogram indices), with $\alpha = 2.0$ .

4. (iv) All distributions have weak (or no) time dependence as the turbulence evolves during the studied period of fully developed turbulence.

5. (v) All distributions have weak (or no) magnetization dependency in the considered range $1 \leqslant \sigma _0 \leqslant 10$ .

6. (vi) In addition, we calibrated the machine-learning-based identification algorithm of Bussov & Nättilä (2021) against the thresholding method. We found that the filling fraction of current sheets in our algorithm matches the thresholding method at a value of $J_{{thr}}/J_{{rms}} = 3$ .

The limitations of our study are mainly due to our algorithm being 2-D. This is a simplification since, in this case, the sheets have a simpler topology to characterize. In three dimensions, the current sheets can have much more complex shapes. Ultimately, however, the analysis must be generalized to full three dimensions to obtain physically realistic sheet distributions.

Many possible extensions of the work can also be devised. For example, it is possible to measure plasma parameters (e.g., plasma- $\beta$ , real magnetization $\sigma$ , temperature $T$ , the ratio of reversing field to guide field $\delta B/B_{0}$ , etc.) upstream from the current sheet using the SCE algorithm, as done in Nurisso et al. (2023), to gather information about the particle acceleration conditions. The sheet statistics should also be studied in different parameter regimes by varying, e.g., the initial amplitude $\delta B/B$ , helicity and plasma composition (e.g., using an electron–ion rather than pair plasma composition). In addition, different driving methods should be explored, such as via wave collisions or purely compressive/solenoidal driving. These are expected to change the turbulence intermittency and affect the sheet statistics. The presented work is, therefore, only the first step in this direction. The differences between the non-relativistic ( $\sigma _0 \ll 1$ ) and ultrarelativistic ( $\sigma _0 \gg 1$ ) regimes are not obvious and will require more parameter studies.

Our findings are relevant to a broad range of astrophysical systems because kinetic/MHD turbulence is also ubiquitous in those systems. Some problems where intermittency of kinetic/MHD turbulence plays a crucial role are as follows.

1. (i) Black hole accretion flows and jets in low luminosity active galactic nuclei. Their environments are expected to be magnetized and turbulent (e.g., Nättilä 2024). Most of the current efforts in supermassive black-hole accretion flow modelling employ general relativistic MHD simulations (Event Horizon Telescope Collaboration et al. 2019). However, ideal fluid simulations cannot self-consistently model the temperature of electrons so the electron temperature profile is often ‘painted on’ via postprocessing (e.g., Mościbrodzka et al. (2016); however, see Galishnikova et al. (2023) for a recent fully kinetic general relativistic PIC simulations with realistic electron plasma physics). Alternatively, some general relativistic MHD simulations employ a variety of subgrid models to compute the electron and ion temperatures based on local PIC simulations of magnetic reconnection (Ball, Sironi & Özel 2018; Werner et al. 2018) or kinetic/MHD turbulent heating based models (Howes 2010; Kawazura, Barnes & Schekochihin 2019; Meringolo et al. 2023).

2. (ii) The Crab pulsar wind nebula (PWN) is another example of a highly magnetized $(\sigma _0 \gtrsim 1)$ environment that can accelerate particles to very high energies (Rees & Gunn 1974; Kennel & Coroniti 1984 a , b ; Atoyan 1999). Recently, in an attempt to resolve observational discrepancies of the Crab PWN with theoretical models, Luo et al. (2020) proposed a turbulent origin of the non-thermal emission in the Crab PWN.

3. (iii) Cosmic ray diffusion in galactic environments, such as interstellar and intercluster mediums, strongly depends on the turbulence physics (Ruszkowski & Pfrommer 2023). It is still not understood what exactly scatters the cosmic rays. However, one hypothesis is the intermittent structures in the turbulence (Kempski et al. 2023; Lemoine 2023).

In summary, we have obtained the scaling distributions for the geometric properties of current sheets, including the length, width and radius of curvature. We find that the length and radius of curvature distributions are well-approximated as power-law distributions and the width distribution with an exponential function. This information can be used to develop realistic subgrid models for turbulent heating in the accretion flows of compact objects and other turbulent environments.