2. Model for GME distribution

Consider a collisionless plasma in a closed system. The evolution of the fine-grained particle distribution for a given species can be represented by the (relativistic) Vlasov equation,

(2.1)\begin{equation} \partial_t f + \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f + \boldsymbol{F} \boldsymbol{\cdot} \partial_{\boldsymbol{p}} f= 0 , \end{equation}

where $f(\boldsymbol {x},\boldsymbol {p},t)$ is the particle momentum distribution function (normalized such that $\int \, {{\rm d}}^{3}p \, {{\rm d}}^{3}x f = N$ is the total number of particles), $\boldsymbol {v} = \boldsymbol {p}c/ (m^{2}c^{2}+p^{2})^{1/2}$ is the particle velocity (with $m$ the particle mass) and $\boldsymbol {F}(\boldsymbol {x},\boldsymbol {p},t)$ is a phase-space conserving force field ($\partial _{\boldsymbol {p}} \boldsymbol {\cdot } \boldsymbol {F} = 0$), containing the electromagnetic force and external forces. Here, $\boldsymbol{x}$ is the position vector, $\boldsymbol{p}$ is the momentum vector, and c is the speed of light. Equation (2.1) can be applied to any particle species, with appropriate $\boldsymbol {F}$. We denote particle kinetic energy by $E(p) = (m^{2}c^{4}+p^{2}c^{2})^{1/2} - m c^{2}$ and the system-averaged kinetic energy by $\bar {E}$.

The Vlasov equation formally conserves the BG entropy $S=-\int \, {{\rm d}}^{3}x \, {{\rm d}}^{3}p f \log {f}$ as well as an infinite set of quantities known as the Casimir invariants. The latter can be manipulated to yield quantities with dimensions of momentum, introduced in Zhdankin (Reference Zhdankin2021a) as the Casimir momenta

(2.2)\begin{equation} p_{c,\chi}(f) \equiv n_0^{1/3} \left( \frac{1}{N} \int \, {{\rm d}}^{3}x \, {{\rm d}}^{3}p f^{\chi} \right)^{{-}1/3(\chi-1)} ,\end{equation}

where $n_0$ is the mean particle number density and $\chi > 0$ is a free index that parameterizes the weight toward different regions of phase space: large (small) values of $\chi$ are sensitive to low (high) energies. The phase-space integral in (2.2) resembles those used in the non-extensive entropies of Rényi (Reference Rényi1961) and Tsallis (Reference Tsallis1988). The Casimir momenta, however, manipulate this integral into a dimensional form that is interpretable physically. In particular, the anomalous growth of $p_{c,\chi }$ is indicative of irreversible entropy production at the corresponding momentum scale in phase space (with $\chi \to 0$ corresponding to momenta far in the tail, and $\chi \to \infty$ corresponding to momenta near the mode).

As described in Zhdankin (Reference Zhdankin2021a), $p_{c,\chi }$ share many properties with the BG entropy $S$: (i) they reduce to a dimensionalized version of the BG entropy when $\chi \to 1$, as $p_{c,\chi \to 1} = n_0^{1/3} \, {\rm e}^{S/3N}$; (ii) they are maximized when $f$ is isotropic and spatially uniform; and (iii) while ideally conserved by the Vlasov equation, the formation of fine-scale structure breaks conservation of $p_{c,\chi }$ for $f$ measured at coarse-grained scales. Zhdankin (Reference Zhdankin2021a) also argued that $p_{c,\chi }$ associated with coarse-grained $f$ will tend to increase (irreversibly) when energy is injected into the system, for generic complex processes; this was demonstrated by two-dimensional kinetic simulations of relativistic turbulence. Phenomena such as the entropy cascade may lead to anomalous entropy production through finite collisionality (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009; Eyink Reference Eyink2018).

The infinite number of generalized entropies represented by $p_{c,\chi }$ complicates the application of a maximum-entropy principle. Only when dissipation occurs collisionally or at small enough energy scales ($\chi \sim 1$) is the BG entropy maximized. In general, mechanisms of anomalous entropy production can operate over a spectrum of scales, so a scale-by-scale understanding of the plasma physical processes is necessary to model the system.

In this Letter, we consider the idealized situation where entropy is maximized at a characteristic momentum scale represented by $p_{c,\chi _d}$ with a given index $\chi _d$ where the subscript $d$ denotes ‘dissipation’. Physically, particles are energized up to this scale (on average) while mixing causes the distribution to smooth out to the equilibrium state.

Suppose that the system evolves to maximize $p_{c,\chi _d}$. The GME distribution is isotropic and spatially uniform $f(\boldsymbol {p},\boldsymbol {x})=f(p)$, and can be derived from the functional

(2.3)\begin{align} {\mathcal{L}} = N^{1/3} \left( \int \, {{\rm d}}^{3}p f^{\chi_d} / N \right)^{{-}1/3(\chi_d-1)} - \lambda_1 \left( \int \, {{\rm d}}^{3}p f - N \right) - \lambda_2 \left[ \int \, {{\rm d}}^{3}p E(p) f - N \bar{E}\right] , \end{align}

where $\lambda _i$ are Lagrange multipliers enforcing number and energy constraints. By requiring $\delta {\mathcal {L}} = 0$ upon variations of the distribution $\delta f$, we obtain

(2.4)\begin{equation} \frac{p_{c,\chi_d}^{3\chi_d-2} \chi_d f^{\chi_d-1}}{3(1-\chi_d) N^{2/3}} - \lambda_1 - \lambda_2E(p) = 0 , \end{equation}

which leads to the GME distribution

(2.5)\begin{equation} f = C [ E(p)/E_b+ 1 ]^{{-}1/(1-\chi_d)} , \end{equation}

where $C$ and $E_b$ are the normalization factor and characteristic energy, determined by requiring $4{\rm \pi} \int dp p^{2} f = N$ and $4{\rm \pi} \int dp p^{2} E(p) f = N \bar {E}$. Note that (2.5) is operationally equivalent to the Tsallis distribution (Tsallis Reference Tsallis1988); this equivalence is due to the fact that the Tsallis entropy and Casimir momenta are both obtained from the same fundamental phase-space integral (involving powers of $f$). We will restrict our attention to $\chi _d < 1$, in which case there is a power-law tail (whereas $\chi _d > 1$ would lead to a narrow distribution with sharp cutoff). The derivation of (2.5) from maximizing a dimensional representation of generalized entropy is the first main result of this work.

In the ultra-relativistic (UR) limit, $\bar {E} \gg m c^{2}$, the GME distribution ((2.5)) becomes

(2.6)\begin{equation} f \xrightarrow[]{\text{UR}} C (p/p_b + 1)^{-\alpha-2}, \end{equation}

where $\alpha = (2\chi _d-1)/(1-\chi _d)$, $C = N (\alpha - 1) \alpha (\alpha + 1)/8{\rm \pi} p_b^{3}$ and $p_b = (\alpha - 2) \bar {E}/3 c$. In the non-relativistic (NR) limit, $\bar {E} \ll m c^{2}$, (2.5) becomes

(2.7)\begin{equation} f \xrightarrow[]{\text{NR}} C ( p^{2}/p_b^{2} + 1 )^{-\alpha-1/2}, \end{equation}

where $\alpha = (1+\chi _d)/2(1-\chi _d)$, $C = N \varGamma (\alpha +1/2)/ {\rm \pi}^{3/2} p_b^{3} \varGamma (\alpha -1)$ and $p_b = [4 (\alpha -2) m \bar {E}/3]^{1/2}$. The NR expression ((2.7)) is equivalent to the kappa distribution.

In both limits, we used $\alpha$ to denote the power-law index of the corresponding energy distribution,

(2.8)\begin{equation} F(E) = \frac{{{\rm d}}p}{{{\rm d}}E} 4{\rm \pi} p^{2} f(p) |_{p=[E (E + 2\,m c^{2})]^{1/2}/c} , \end{equation}

such that $F(E) \propto E^{-\alpha }$ at high energies. Also note that $\chi _d \to 1$ ($\alpha \to \infty$) recovers the thermal (Maxwell–Jüttner) distribution, using the identity $(A/x + 1)^{-x} = {\rm e}^{-A}$ as $x \to \infty$ for any $A$.

Since the GME distribution has an infinite extent in energy, $\alpha > 2$ is necessary for finite $\bar {E}$. Thus, the domain is $3/4 < \chi _d < 1$ for the UR case and $3/5 < \chi _d < 1$ for the NR case. We note that the GME framework can be extended to allow $1 < \alpha < 2$ if an additional constraint is imposed to make the distribution vanish at a maximum momentum $p_{\rm max}$ (which may be related to the system confinement scale, for example). This would be implemented by adding a third Lagrange multiplier to ${\mathcal {L}}$ that enforces $p_{c,\chi \to 0} = p_{\rm max}$. However, the resulting equation does not have an analytically tractable solution for $f$, so we defer such an extension to future work. This extension may be necessary to accurately model particle distributions in relativistic magnetic reconnection (Sironi & Spitkovsky Reference Sironi and Spitkovsky2014) or turbulence (Zhdankin et al. Reference Zhdankin, Werner, Uzdensky and Begelman2017) at high magnetization, where $\alpha < 2$ has been measured in simulations.

3. Model for power-law index

Suppose that the dynamics is sufficiently complex to cause the initial distribution (which is arbitrary) to evolve into the GME state. One can then compare the momentum at which entropy is maximized, $p_{c,\chi _d}$, with the momentum of the typical particle given by $p_{c,\infty }$ (note that $p_{c,\infty }$ lies close to $p_b$). Evaluating $p_{c,\chi _d}/p_{c,\infty }$ using (2.2) with the GME distribution ((2.5)), one obtains in the UR limit

(3.1)\begin{equation} \frac{p_{c,\chi_d} }{p_{c,\infty}} \xrightarrow[]{\text{UR}} \left(\frac{\alpha+1}{\alpha-2} \right)^{(\alpha+2)/3} , \end{equation}

and in the NR limit

(3.2)\begin{equation} \frac{p_{c,\chi_d} }{p_{c,\infty}} \xrightarrow[]{\text{NR}} \left( \frac{\alpha-1/2}{\alpha-2} \right)^{(2\alpha+1)/6} . \end{equation}

This relates the power-law index $\alpha$ to the maximum-entropy scale, which can be modelled phenomenologically (as considered below). In figure 1, we show $\alpha$ versus $p_{c,\chi _d}/p_{c,\infty }$, separately for the UR ((3.1)) and NR ((3.2)) limits. Note the divergence $\alpha \to \infty$ when $p_{c,\chi _d}/p_{c,\infty } \to e \approx 2.72$ (UR case) or $p_{c,\chi _d}/p_{c,\infty } \to \, {\rm e}^{1/2} \approx 1.65$ (NR case). Thus, if entropy is maximized at momentum scales sufficiently close to the peak of the distribution, then a thermal distribution is recovered (similar to a collisional plasma). When $p_{c,\chi _d}/p_{c,\infty }$ becomes larger than a factor of few, the non-thermal state is obtained, with $\alpha \to 2$ for $p_{c,\chi _d}/p_{c,\infty } \gg 1$. Thus, in both the UR and the NR limit, the distribution will relax to the non-thermal state if entropy is maximized at momentum scales in the tail of the distribution.

Figure 1. The energy power-law index $\alpha$ of the GME distribution versus the ratio between the entropy-maximizing momentum $p_{c,\chi _d}$ and the typical momentum $p_{c,\infty }$ ((3.1) and (3.2)). The UR (red) and NR (blue) limits are shown separately, with dashed lines indicating singularities.

Physical considerations are necessary to determine $p_{c,\chi _d}/p_{c,\infty }$ as a function of system parameters, from which one can extract $\alpha$. In general, this will need to be informed by numerical simulations and analytical considerations for the given process.

For this Letter, we consider a simplified scenario to estimate the momentum scale of maximum entropy that arises from the dissipation of magnetic energy in complex field topologies (via magnetic reconnection, turbulence or instabilities). We suppose that, rather than being energized at the thermal energy scale, the typical particles are energized by an amount comparable to the free magnetic energy per particle, $E_{\rm free} = \delta B^{2}/8{\rm \pi} n_0$, over a dynamical time scale, before equilibrating to the GME state. Here, $\delta B$ is the characteristic magnetic field fluctuation (prior to dissipation), while we denote the background (non-dissipating) component by $B_0$. We denote the energy corresponding to the Casimir momenta by $E_{c,{\chi }} = E(p_{c,\chi })$ and the typical particle energy as $E_0$ (prior to dissipation), noting that the thermal dissipation energy scale is $e E_0$. The model posits that $E_{c,\chi _d} \sim e E_0 + \eta E_{\rm free}$ where $\eta$ is an order-unity coefficient describing the portion of free energy converted. We can then write

(3.3)\begin{align} \frac{p_{c,\chi_d}}{p_{c,\infty}} & = \left[ \frac{E_{c,\chi_d} ( E_{c,\chi_d} + 2\,m c^{2})}{E_{c,\infty} (E_{c,\infty} + 2\,m c^{2})} \right]^{1/2} \nonumber\\ & \sim \left[ \frac{(e E_0 + \eta E_{\rm free}) (e E_0 + \eta E_{\rm free}+2\,m c^{2})}{E_0 (E_0+2\,m c^{2})} \right]^{1/2} \nonumber\\ & \sim \left[ \frac{[e + \eta (\delta B/B_0)^{2}/\beta_c] [e+ \eta (\delta B/B_0)^{2}/\beta_c+2/\theta_c]}{1+2/\theta_c} \right]^{1/2}, \end{align}

where $\theta _c = E_0/m c^{2}$ is a characteristic dimensionless temperature and $\beta _c = 8{\rm \pi} n_0 E_0/B_0^{2}$ is a characteristic plasma beta for the particle species (which may differ from the standard plasma beta, $\beta _0=8{\rm \pi} n_0 T/B_0^{2}$ where $T$ is species temperature, by a factor of order unity). Equating (3.3) with either (3.1) or (3.2) yields an implicit equation for $\alpha$ as a function of $\beta _c$, $\delta B/B_0$ and $\theta _c$ in the appropriate limit. The physical parameters required to achieve a given value of $\alpha$ can then be expressed in the UR limit ($\theta _c \gg 1$) as

(3.4)

and in the NR limit ($\theta _c \ll 1$) as

(3.5)

The predicted scaling of $\alpha$ given by (3.4) and (3.5) is the second main result of this work. The right-hand side of both equations becomes zero when $\alpha \to \infty$, indicating that the thermal distribution is recovered for high beta or weak fluctuations, $(\delta B/B_0)^{2} / \beta _c \ll 1$. On the other hand, the non-thermal state is obtained when $(\delta B/B_0)^{2}/\beta _c \gtrsim 1$, for both UR and NR regimes. For $(\delta B/B_0)^{2}/\beta _c \gg 1$, the index approaches $\alpha \to 2$ (but recall that the model may be extended, in principle, to allow $1 < \alpha < 2$). The scaling is plotted in figure 2.

Figure 2. The energy power-law index $\alpha$ of the GME distribution versus (pre-dissipation) physical parameters $\eta (\delta B/B_0)^{2} / \beta _c$ for the magnetic dissipation model. The UR (red; (3.4)) and NR (blue; (3.5)) limits are shown separately.

4. Comparison with simulations

To validate the GME model, we remark on how the predictions compare with existing results from kinetic simulations of relativistic turbulence and magnetic reconnection in the literature.

In figure 3, we show the global particle energy distribution $F(E)$ arising in a $1536^{3}$-cell particle-in-cell (PIC) simulation of driven relativistic turbulence (with $\delta B/B_0 \approx 1$) studied in  Zhdankin et al. (Reference Zhdankin, Uzdensky, Werner and Begelman2018) and Wong et al. (Reference Wong, Zhdankin, Uzdensky, Werner and Begelman2020). The simulation begins with a Maxwell–Jüttner distribution of electrons and positrons with UR temperature $\theta =T/m_e c^{2}=100$ and initial magnetization $\sigma _0 = 3/8$. The magnetization is defined as the ratio of the magnetic enthalpy to plasma enthalpy, and is related to species plasma beta by $\sigma _0 = 1/(4\beta _0)$ in the UR regime; thus $\beta _0 = 2/3$. The simulation develops a non-thermal tail with index $\alpha \approx 3$. We find that the GME distribution of (2.6) provides a fair fit to the fully developed state when we choose $\chi _d = 0.815$, as shown by the dashed line in figure 3. The fit over-predicts the number of particles at energies below the peak, indicating that relaxation to the GME state is incomplete (possible reasons for this will be described in the conclusions). The PIC simulations of decaying, magnetically dominated turbulence by Comisso & Sironi (Reference Comisso and Sironi2019) also appear to resemble the GME state. Thus, we believe that the GME model provides a reasonable (if imperfect) representation of available numerical data on relativistic turbulence.

Figure 3. Energy distribution $F(E)$ in PIC simulation of relativistic turbulence for various times, taken from Zhdankin et al. (Reference Zhdankin, Uzdensky, Werner and Begelman2018), compared with the GME distribution (dashed; (2.6)) with $\chi _d = 0.815$.

We next consider the model for the power-law index $\alpha$ from magnetic dissipation. In figure 4, we compare the predicted $\alpha$ versus $\sigma$ scaling ((3.4) with $\beta _c = 1/4\sigma$, $\delta B/B_0 = 1$, $\eta = 1$) with results in the literature on relativistic turbulence in a pair plasma. PIC simulations of driven relativistic turbulence indicate that the power-law index is well described by the empirical formula $\alpha \approx \alpha _\infty + C_0 \sigma ^{-0.5}$, with $\alpha _\infty \approx 1$ and $C_0 \approx 1.5$ for large sizes (Zhdankin et al. Reference Zhdankin, Werner, Uzdensky and Begelman2017, Reference Zhdankin, Uzdensky, Werner and Begelman2018), shown in figure 4 (blue); note that a similar formula with different coefficients was also suggested for relativistic magnetic reconnection (Ball, Sironi & Özel Reference Ball, Sironi and Özel2018; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018; Uzdensky Reference Uzdensky2022). We also show approximate data points from the two-dimensional decaying relativistic turbulence simulations of Comisso & Sironi (Reference Comisso and Sironi2019) (red). The model is able to explain the trends in the numerical simulations fairly well, up to a factor of order unity in $\sigma$. Fits to the simulation data can be improved by adjusting $\eta$, noting that driven turbulence would effectively have a larger $\eta$ than decaying turbulence. Additionally, we note that Comisso & Sironi (Reference Comisso and Sironi2019) find that $\alpha$ increases with decreasing $\delta B/B_0$, consistent with the GME prediction.

Figure 4. Energy power-law index $\alpha$ versus magnetization $\sigma$ from the GME model in the UR limit (black; (3.4) with $\beta _c = 1/4\sigma$, $\delta B/B_0 = 1$, and $\eta = 1$) compared with empirical fitting formula $\alpha \approx \alpha _\infty + C_0 \sigma ^{-0.5}$ from PIC simulations of driven relativistic turbulence in Zhdankin et al. (Reference Zhdankin, Werner, Uzdensky and Begelman2017) (blue). Also shown is the approximate range of indices from PIC simulations of decaying relativistic turbulence from Comisso & Sironi (Reference Comisso and Sironi2019) (their figure 5 inset; red).

In addition to these quantitative comparisons, the GME model provides a resolution to several mysterious findings from kinetic simulations in the literature. Kinetic simulations of disparate processes (turbulence, magnetic reconnection and instabilities) often exhibit very similar power-law distributions for given plasma parameters. For example, PIC simulations find comparable non-thermal particle acceleration from magnetic dissipation with different current sheet geometries and ensuing dynamics (e.g. Werner & Uzdensky Reference Werner and Uzdensky2021). PIC simulations of relativistic turbulence find that non-thermal particle distributions have a similar shape for different driving mechanisms (electromagnetic, solenoidal, compressive, imbalanced), despite different time scales to arrive at those distributions (Zhdankin Reference Zhdankin2021b; Hankla et al. Reference Hankla, Zhdankin, Werner, Uzdensky and Begelman2022). The universality revealed by these findings may be explained by all of the processes having sufficient complexity to attain a GME state at similar energy scales.

Kinetic simulations also indicate that non-thermal particle distributions formed by relativistic magnetic reconnection (Werner & Uzdensky Reference Werner and Uzdensky2017; Guo et al. Reference Guo, Li, Daughton, Li, Kilian, Liu, Zhang and Zhang2021) and turbulence (Comisso & Sironi Reference Comisso and Sironi2019) are insensitive to the number of spatial dimensions (two dimensions vs three dimensions), despite different secondary instabilities, cascades, and trapping mechanisms (e.g. long-lived plasmoids in two dimensions). The GME framework predicts that the distributions are insensitive to the number of spatial dimensions, as long as there are sufficient degrees of freedom to attain such a state.

The GME model predicts similar acceleration efficiency in the NR regime as in the UR regime, as long as no factors arise that suppress entropy production at high energies. PIC simulations in the NR regime are generally constrained in scale separation, which may limit power-law formation. However, recent PIC simulations of NR magnetic reconnection provide some evidence for (steep) power-law distributions at low $\beta$ (Li et al. Reference Li, Guo, Li, Stanier and Kilian2019). Recent simulations of reduced kinetic models indicate efficient electron acceleration by NR magnetic reconnection at macroscopic scales when $\delta B/B_0$ is large enough (Arnold et al. Reference Arnold, Drake, Swisdak, Guo, Dahlin, Chen, Fleishman, Glesener, Kontar and Phan2021). Hybrid kinetic simulations of turbulence driven by the magnetorotational instability described by Kunz et al. (Reference Kunz, Stone and Quataert2016) exhibit an ion distribution that is well fit by a kappa distribution, as predicted by the GME model ((2.7)), although the index appears to be harder than predicted by (3.5) for the high values of plasma beta (possibly a consequence of non-magnetic sources of free energy). Hybrid kinetic simulations of Alfvénic turbulence may provide further tests of the model; published cases with $\delta B/B_0 \ll 1$ and moderate beta do not exhibit significant particle acceleration, consistent with (3.5) (e.g. Arzamasskiy et al. Reference Arzamasskiy, Kunz, Chandran and Quataert2019; Cerri, Arzamasskiy & Kunz Reference Cerri, Arzamasskiy and Kunz2021). Further benchmarking of the model in the NR regime is deferred to future work.