Relativistic Alfvén Turbulence at Kinetic Scales

In a strongly magnetized, magnetically dominated relativistic plasma, Alfvénic turbulence can extend to scales much smaller than the particle inertial scales. It leads to an energy cascade somewhat analogous to inertial- or kinetic-Alfvén turbulent cascades existing in nonrelativistic space and astrophysical plasmas. Based on phenomenological modeling and particle-in-cell numerical simulations, we propose that the energy spectrum of such relativistic kinetic-scale Alfvénic turbulence is close to k −3 or slightly steeper than that due to intermittency corrections or Landau damping. We note the analogy of this spectrum with the Kraichnan spectrum corresponding to the enstrophy cascade in 2D incompressible fluid turbulence. Such turbulence strongly energizes particles in the direction parallel to the background magnetic field, leading to nearly one-dimensional particle momentum distributions. We find that these distributions have universal log-normal statistics.


INTRODUCTION
Large-scale astrophysical flows are often hydrodynamically unstable, which leads to the excitation of velocity, density, and electromagnetic fluctuations spanning a broad range of spatial and temporal scales.A significant fraction of energy associated with large-scale motions can then be converted into random collective plasma fluctuations and eventually dissipated at significantly smaller scales due to non-ideal processes, such as particle collisions or wave-particle interactions.In a weakly collisional plasma, the dissipated turbulent energy is converted into heat or non-thermally accelerated particles.The particle distribution functions can significantly deviate from a Maxwellian, which affects plasma dynamics and thermodynamics, as well as the radiative signatures of astrophysical objects (e.g., Drake et al. 2013;Sironi & Spitkovsky 2014;Zhdankin et al. 2017;Comisso & Sironi 2019;Wong et al. 2020;Nättilä & Beloborodov 2021;Demidem et al. 2020;Trotta et al. 2020;Pezzi et al. 2022;Vega et al. 2022a;Bresci et al. 2022;Comisso & Sironi 2022;Vega et al. 2023).
Since astrophysical plasmas are typically magnetized, large-scale fluctuations are often dominated by lowfrequency shear-Alfvén modes.The energy dissipation occurs at much smaller, kinetic scales that are com-Corresponding author: Cristian Vega csvega@wisc.eduparable to plasma microscales such as particle inertial and gyro scales.At such scales, the shear-Alfvén modes transform into kinetic-Alfvén or inertial-Alfvén modes.Energy cascade in the kinetic range arguably governs the energy dissipation and particle heating in a collisionless turbulent plasma, and it may be relevant for nonthermal particle acceleration.Kinetic-range modes play an important role in space and astrophysical plasmas where they have been studied analytically, numerically, and where possible, in in situ spacecraft measurements (e.g., Chen et al. 2013;Sahraoui et al. 2013;Chen et al. 2014;Told et al. 2015;He et al. 2020;Chen & Boldyrev 2017;Roytershteyn et al. 2019;Zhou et al. 2023;Vega et al. 2023;Mallet et al. 2023).
Relativistic plasma turbulence has several features that are qualitatively different from the non-relativistic case.First, if the plasma temperature is relativistic, the Alfvén speed as well as plasma frequency and gyrofrequency depend on temperature, since instead of the rest-mass density, their expressions include particle energy.Second, in such a plasma the speed of sound is comparable to the speed of light.Even when the magnetic energy exceeds the relativistic kinetic energy of the particles, which is a situation somewhat analogous to a non-relativistic low-beta case, thermal corrections to the wave dispersion relation remain significant (e.g., Ten-Barge et al. 2021;Vega et al. 2022a).This complicates the analytical consideration of relativistic turbulence.
It has been observed in numerical simulations of strongly magnetized relativistic plasma that the Alfvénic turbulence cascade indeed continues to the scales smaller than the particle inertial scales (e.g., Zhdankin et al. 2018;Comisso & Sironi 2019;Vega et al. 2022a).However, as such scales are typically not well resolved in studies of large-scale Alfvénic turbulence, the measurements of the energy spectrum and other statistical properties of kinetic-scale fluctuations remain inconclusive.In this work, we study numerically and analytically kinetic-range Alfvénic turbulence in ultra-relativistic electron-positron plasma immersed in a strong background magnetic field, B 0 ≫ δB 0 .We discuss the spectrum of magnetic and electric fluctuations, as well as the intermittency properties of turbulence.As relativistic turbulence leads to efficient particle acceleration, we also discuss the resulting nonthermal distribution function of the accelerated particles.
In what follows, we consider a pair plasma and reserve the standard notations, d e = c/ω pe , ω pe = 4πn 0 e 2 /m e , and Ω ce = |e|B 0 /m e c for the nonrelativistic electron inertial scale, electron plasma frequency, and electron cyclotron frequency, correspondingly.Their relativistic counterparts are not unique but rather depend on the wave mode considered.Their corresponding definitions will be given in the text.

ANALYTICAL CONSIDERATION
Consider a relativistic electron-positron plasma with an imposed uniform background magnetic field B 0 = B 0 ẑ.The plasma magnetization with respect to the background field is defined as where n 0 denotes the unperturbed number density of electron or positron species, and n 0 wm e c 2 is the corresponding enthalpy density.Assuming that the plasma particle distribution is an isotropic Maxwell-Jüttner function with temperature T , the enthalpy is calculated as , where K ν is the modified Bessel function of the second kind.In this formula, θ = k B T /m e c 2 is the normalized temperature.Similarly, one can define the magnetization with respect to the fluctuating part of the field, σ = (δB) 2 4πn 0 wm e c 2 . (2) In this work, we assume that plasma is magnetically dominated, σ ≫ 1, and the magnetic fluctuations are smaller than the background field, σ ≪ σ.In our numerical simulations, we initialize the runs with magnetic fluctuations satisfying σ0 ≫ 1.However, as turbulence evolves, turbulent fluctuations efficiently heat the plasma, so plasma temperature becomes ultrarelativistic while simultaneously, σ decreases.This reflects the fact that relativistic turbulent motion is inherently compressible, which allows colliding fluid elements to convert their kinetic energy into heat rapidly (Zhdankin et al. 2018;Nättilä & Beloborodov 2021;Vega et al. 2022bVega et al. , 2023).We will therefore analyze the case when plasma bulk fluctuations are non-relativistic (or mildly relativistic), while plasma temperature is ultrarelativistic. 1ecently, a two-fluid analytic model was proposed to describe Alfvénic turbulence in a magnetically dominated relativistic plasma (Vega et al. 2022a).This model assumes that the fluctuations are anisotropic in the Fourier space with respect to the background magnetic field, k z ≪ k ⊥ , and the magnetic and electric fluctuations are dominated by their field-perpendicular components.The equations require a closure for the plasma pressure tensor.There is no rigorous prescription for such a closure in a collisionless plasma.In the magnetically dominated case considered here, only the field-parallel component of the pressure tensor contributes.We may then introduce the acoustic velocity as v 2 s = c 2 (∂P ∥ /∂u)| 0 , where u is the thermal energy density.(For relativistic plasma temperatures, the acoustic velocity is on the order of the speed of light.For example, in the case of the 3D ultrarelativistic isotropic adiabatic equation of state, one has P ≈ u/3 and v 2 s ≈ c 2 /3.)The model equations govern Alfvénic turbulence in an ultra-relativistic (θ ≫ 1) magnetically dominated (σ ≫ 1) pair plasma.It is convenient to normalize the electric potential and the z-component of the magnetic vector potential as φ = ϕc/B 0 and Ãz = A z c/(B 0 1 + 2/σ).The perpendicular components of the magnetic and velocity fluctuations are defined as δ b⊥ = − ẑ × ∇ Ãz and δ ṽ⊥ = ẑ ×∇ φ, and they both have the dimensions of velocity.Below we will use only these variables and omit the overtilde sign.The equations then have the form (Vega et al. 2022a): where is the relativistic Alfvén speed in a pair plasma, is the relativistic inertial scale, and the magnetic-fieldparallel gradient is given by The extra factor of 2 in the expressions for v A and d rel reflects the fact that the total plasma density is twice the electron or positron density.The very last term in the right-hand side of Eq. ( 4), proportional to v 2 s /c 2 , arises from the hydrodynamic closure for the parallel-pressure term discussed above.
In the linear case, these equations describe the waves with the dispersion relation (Vega et al. 2022a): Similarly to the kinetic-Alfvén waves, the numerator of this expression involves the contribution of the thermal effects.Similarly to the inertial Alfvén waves, the denominator includes the contribution of the electron (and positron) inertia.This is somewhat analogous to the inertial kinetic-Alfvén modes in a non-relativistic plasma (e.g., Streltsov & Lotko 1995;Lysak & Lotko 1996;Chen & Boldyrev 2017;Roytershteyn et al. 2019;Loureiro & Boldyrev 2018;Boldyrev et al. 2021).In contrast with the non-relativistic case, however, in a relativistic plasma, we have v s ∼ c, so that the thermal contribution in the numerator is never negligible.Rather, the thermal and inertial effects in Eq. ( 8) are necessarily of the same order.Moreover, since the speed of sound is comparable to the thermal speed, the Landau damping of the linear modes is also generally strong at k 2 ⊥ d 2 rel ≳ 1.As discussed in the Appendix, the applicability of the linear dispersion relation in Eq. ( 8) at such scales depends on the particle distribution function.Such a function may, in general, be strongly non-thermal.
The spectrum of turbulence at small, kinetic scales, rel ≫ 1, is currently not well understood.The description of relativistic turbulence at such scales faces both numerical and analytical challenges.On the numerical side, available particle-in-cell kinetic simulations do not typically have a strong enough guide field B 0 in order to separate the inertial and gyro scales, or large enough numerical resolution (number of cells and particles) in order to reliably measure the spectra in the sub The reported electromagnetic spectra at the sub-inertial scales, varied depending on the strength of the guide field.The spectrum was consistent with W k ∝ k −4.5 for relatively weak guide fields, B 0 /δB 0 ≲ 1, and flattened to approximately W k ∝ k −3.5 for stronger fields, B 0 /δB 0 ∼ 3 (e.g., Zhdankin et al. 2018;Comisso & Sironi 2018;Nättilä & Beloborodov 2021;Vega et al. 2022a,b).We will be interested in the limit of strong guide field, B 0 /δB 0 ≫ 1; our numerical simulations discussed below will correspond to B 0 /δB 0 = 10 and σ 0 = 4000.On the analytical side, at sub-inertial scales, k 2 ⊥ d 2 e ≫ 1, the model equations (3, 4) formally conserve the energy integral (Vega et al. 2022a), This quantity is expected to cascade in a Fourier space toward large wave numbers, somewhat analogously to the enstrophy cascade in 2D incompressible hydrodynamic turbulence (Kraichnan 1967).Dimensional arguments then predict the electromagnetic energy spectrum Similarly to the hydrodynamic case, the energy cascade is expected to have intermittency corrections that lead to a steeper energy spectrum, Here, k 0 is the large-scale boundary of the spectrum, which approximately corresponds to the inverse electron inertial scale.The intermittency correction reflects the non-locality of turbulence.The electromagnetic energy spectrum close to k −3 ⊥ implies that vorticity and current structures at scales k ⊥ ≫ k 0 are strained most efficiently by turbulent eddies at scales k 0 (e.g.Boffetta & Ecke 2012).
As discussed above, significant Landau damping may affect relativistic Alfvén turbulence at kinetic scales.We, however, conjecture that as a consequence of the non-locality of turbulence, the spectrum should exhibit a near power-law behavior, close to that given by Eqs.(11,11).Indeed, the Kraichnan spectrum is established due to gradient-stretching of small-scale structures by turbulent eddies of the scale k 0 ∼ 1/d rel .As a result, all small-scale modes have the same evolution time and the same parallel phase velocity.A particle resonating with structures of scales k ⊥ ≫ 1/d rel then essentially resonates with an entire eddy of scale k 0 ∼ 1/d rel .Landau damping may therefore regulate the overall intensity of kinetic-scale fluctuations, while not significantly affecting their spectrum.Our numerical results analyzed below seem to be consistent with this prediction.
Finally, we mention a formal restriction on a turbulent cascade in a magnetically dominated plasma, which applies in both relativistic and non-relativistic cases.When the magnetization parameter σ 0 is large, 2 magnetic fluctuations at small enough scale λ s may enter the so-called "charge-starved" regime when the scale satisfies (e.g., Thompson 2006Thompson , 2008;;Blaes et al. 1990;Boldyrev et al. 2021;Chen et al. 2022;Nättilä & Beloborodov 2022): At such scales, the field-parallel electric current would correspond to the electrons moving with the speed of light.In this regime, the non-relativistic equation ( 4) for parallel electron dynamics cannot be used.Since the current is restricted by the speed of light, its autocorrelation function is limited by a constant, We then conclude that in the asymptotic limit kλ s ≫ 1, the Fourier spectrum of the current cannot be flatter than k −1 ⊥ and, correspondingly, the spectrum of the magnetic field cannot be flatter than k −3 ⊥ .The spectra given by Eqs.(11,12) satisfy this restriction.
To address this problem numerically, we performed two 2.5D particle-in-cell simulations of decaying turbulence in a pair plasma with the fully relativistic code VPIC (Bowers et al. 2008).In 2.5D simulations, vector fields have three components but depend on only two spatial coordinates x and y.Similarly, the particle distribution function depends on three velocity components.A uniform background magnetic field is imposed in the z-direction, B 0 = B 0 ẑ.The simplified "2.5"dimensional setup allows us to afford a relatively high numerical resolution of kinetic-range turbulence.Since all the vector components of the electromagnetic field and particle momenta are preserved, it is expected to capture some essential nonlinear interactions existing in the 3D case.Previous numerical studies involving 2.5D and 3D runs seem to produce similar energy spectra of fields and particles (e.g., Zhdankin et al. 2017Zhdankin et al. , 2018;;Comisso & Sironi 2018, 2019;Vega et al. 2023).
Table 1 summarizes the parameters of the simulations.Run I is a large-box, high-resolution simulation that spans both hydrodynamic and kinetic scales, with about 15 cells per nonrelativistic d e .Run II is a small-box simulation where the number of cells per nonrelativistic d e was increased to 80, drastically improving the resolution of the sub-d e fluctuations while decreasing the hydrodynamic range.Both simulation domains are double periodic L × L squares with 100 particles per cell per species.The turbulence is seeded by randomly phased magnetic perturbations of the Alfvénic type where the unit polarization vectors are normal to the background field, ξk = k × B 0 /|k × B 0 |.The wave vectors of the modes are given by k = {2πn x /L, 2πn y /L}, where n x , n y = 1, ..., n max .All modes have the same amplitudes δB k .The time is normalized to the largescale crossing time l/c, where c is the speed of light and the outer scale of turbulence is evaluated as l = L/n max , with n max = 4 or 8 depending on the considered run, see Table 1.
The initial distribution of plasma particles is an isotropic Maxwell-Jüttner distribution with the temperature parameter θ 0 = 0.3.Here, θ 0 = k B T 0 /m e c 2 is the  normalized initial temperature.For θ 0 = 0.3, one has w 0 ≈ 1.88, for ultrarelativistic temperatures, θ 0 ≫ 1, one has w 0 ≈ 4θ 0 , while for non-relativistic plasma, θ 0 ≪ 1, w 0 ≈ 1.In the large-box run I, an initial plasma current is also added to the system with the goal to compensate for the curl of the initial magnetic perturbations, J z = (c/4π)∇ × δB 0 .This helps to avoid the generation of high-frequency ordinary modes with nonzero E z in addition to the low-frequency Alfvén modes.To add the current, the initial plasma density n 0 is kept uniform, and velocity U s z = J z /(2q s n 0 ) is added (in a Newtonian way) to each particle of species s with charge q s = ±e (positrons and electrons) sampled from the Maxwell-Jüttner distribution, provided |v s z + U s z | < c.The distribution is unchanged in the regions where such a velocity increase would lead to |v s z + U s z | ≥ c.The addition of such a current does not change the core of the particle distribution function but modifies its highenergy tail, as will be seen below.
In the small-box run II, the addition of a compensating current is less practical, since it would formally require the electron velocities to exceed the speed of light in more cells.We do observe the generation of a weak ordinary mode in this case.The presence or absence of the initial compensating current, however, does not qualitatively affect the particle distribution function and the spectra of Alfvénic fluctuations eventually generated by the developed turbulence.
Figure 1 illustrates the time evolution of the electromagnetic energies in both runs.Our statistical analysis of field and particle distribution functions is per-  formed after the initial relaxation is completed, that is, after quasi-steady states of electric and magnetic fluctuations are reached.Figure 2 presents a visualization of magnetic fluctuations in both runs, which shows similar morphologies of the magnetic structures generated in the large-and small-box simulations.
In both runs, turbulent fluctuations lead to significant particle energization.The particle energy distribution functions evolve fast until about half of the initial magnetic energy has been converted into the kinetic energy of the particles.After that, the distributions reach quasi-saturated states, as shown in Figure 3.In agreement with previous studies (Nättilä & Beloborodov 2022;Vega et al. 2022b), the particle momentum distribution functions are strongly non-thermal (non-Maxwell-Jüttner).Moreover, they are strongly anisotropic with respect to the background magnetic field, see Figure 4.The electrons are energized mostly in the field-parallel direction, leading to a nearly onedimensional particle distribution function.
We find that for γ > 1, the particle energy distribution functions can be well approximated by the log-normal distribution: where A is a normalization constant.Here, we use the fact that the particle energy can be represented as γmc 2 , so that the relativistic energy distribution can be characterized by the distribution of the corresponding Lorenz factors γ.The parameters µ and σ 2 s of the log-normal fits for each run are given in Table 2, which also shows the measured average Lorenz factors of the electron energies. 3To find the best fits, we varied the parameters µ and σ 2 s with an imposed constraint that the average Lorenz factors, ⟨γ⟩ = exp µ + σ 2 s /2 , have the fixed values presented in the Table .Table 2 also shows the ratio of the field-parallel and field-perpendicular components of the electron pressure tensor.Since the bulk plasma motion is only mildly relativistic, the ratio of the electron pressure tensor components is numerically calculated as the ratio of the components of the electron energy-momentum tensor, P ⊥ /P ∥ ≈ T ⊥ /T zz .
We also note that the log-normal particle energy distributions generated by turbulence with a strong guide field are different from power-law distributions numerically observed in turbulence with a weak guide field (e.g., Zhdankin et al. 2017Zhdankin et al. , 2018;;Comisso & Sironi 2019;Vega et al. 2022b).This may indicate that different particle acceleration mechanisms operate in these two regimes.In the case of a strong guide field, a particle gyroradius is much smaller than the plasma inertial scale.In this limit, the acceleration may be provided by the electric field parallel to the background magnetic field or by the particle curvature drifts in the vicinity of strong magnetic structures.As a result, particles are accelerated in the field-parallel direction.Particles with sufficiently high energies can however have gyroradii larger than the plasma inertial scale (e.g., Vega et al. 2022b), which happens when γc/Ω ce > d rel , that is, Such particles interact with Alfvénic turbulent fluctuations more efficiently.Their magnetic moments may be not conserved, and they can be accelerated in the field-perpendicular direction.As a result, their energy distribution functions develop power-law tails.Since in our simulations, w 0 ≈ 1.88, σ 0 = 4000, and w ≈ 2⟨γ⟩ (see Eq. ( A18)), we obtain a rather large value of the critical energy, γ crit ≈ 274.This may explain why we observe only the log-normal part of the distributions.
As the initial field-perpendicular magnetic perturbations relax, they drive turbulence at large scales.In a magnetically dominated plasma, the excited large-scale fluctuations can be a combination of the two modes whose magnetic polarizations are normal to the background field: the shear-Alfvén mode and the ordinary mode.The frequency of the ordinary mode is given by where ω 2 p = 2ω 2 pe ⟨1/γ 3 ⟩, see Eqs. (A12) and (A15) in Appendix.This mode is excited in our setup with a relatively low amplitude, contributing only a small fraction of the total turbulent energy.Such a mode is not important at the hydrodynamic scales.
The energies of electric and magnetic fluctuations in the large-box run are shown in Fig. 5.The total energy spectrum is slightly steeper than k −3 , however, it is consistent with the Kraichnan spectrum of turbulence including a logarithmic intermittency correction, in agreement with Eq. ( 12).We find that k −1 0 = 1.2 d e provides a good match.Also, as we mentioned before, Landau damping may play a role in the steepening of the spectrum.
In a the small-box run, the imposed magnetic fluctuations at t = 0 are relatively strong at small scales; The measurements of the coefficients of the log-normal fit in Eq. ( 16), the average electron Lorentz factor ⟨γ⟩, as well as the components of the electron pressure tensor, are made at ct/l = 76 for run I and ct/l = 33 for run II.
their decay leads to the production of a weak ordinary mode.Such a mode is most strongly generated at the smallest scale available for Aflvénic fluctuations since the current is largest there.This is the inertial scale of a non-relativistic pair plasma, expressed as We will see below that the energy of the ordinary modes is indeed concentrated at these kinetic scales that are the focus of our study.Since the phase velocity of the ordinary mode exceeds the speed of light, see Eq. ( 18), such fluctuations are not significantly damped.Therefore, we need to make sure that in our statistical analysis, the fluctuations associated with such a mode are separated from the Alfvénic fluctuations.
The ordinary mode can be detected in simulations if we observe that its electric field is polarized in the z direction, while the electric polarization of the shear-Alfvén mode is normal to the z direction.The frequency of the E z fluctuations can be numerically found from the Maxwell-Ampere law by calculating the ratio: Figure 6 shows this frequency calculated for the smallbox run.The frequency indeed agrees with the analytic dispersion relation given by Eq. ( 18).Since k ⊥ ≫ k z , the frequency of the ordinary mode is much larger than the frequency of shear-Alfvén fluctuations.The oblique shear-Alfvén fluctuations are characterized by a mostly potential electric field.The frequency of the potential electric fluctuations can be similarly measured from the Maxwell-Ampere law: Fig. 6 illustrates that it is indeed much smaller than the frequency of the ordinary mode.We also note that the frequency calculated in Eq. ( 21) coincides with the frequency of electric charge fluctuations since the divergence of the electric field coincides with the charge density.21).Here, we denote dp ≡ c/ωp, where the relativistic plasma frequency ωp is defined in Eq. ( 18), and Ω ce,rel ≡ Ωce/⟨γ⟩, where Ωce is the non-relativistic electron cyclotron frequency and ⟨γ⟩ is given in Table 2.
Run U e z,rms Uz,rms U ⊥,rms ⟨γ e ⟩ I 0.47c 0.070c 0.32c 0.053c 1.25 II 0.61c 0.079c 0.39c 0.056c 1.73 Table 3. Parameters of the velocity fluctuations for the electron fluid, and for the bulk plasma motion, U = (niU i + neU e )/(ni + ne).Here, γe is the Lorenz factor associated with the electron fluid velocity.The measurements are made at ct/l = 76 for run I and ct/l = 33 for run II.
The nearly linear scaling of the measured frequency of potential fluctuations with the wavenumber can be a consequence of two independent effects.First, it may reflect the almost linear dispersion relation of the Alfvén mode given by Eq. ( 8).Indeed, in our 2D runs, the magnetic-field direction deviates from the z-axis by a small angle, sin(θ) ∼ δB ⊥ /B ∼ 0.1.Alfvén fluctuations with the wavenumber k ⊥ then correspond to the field-parallel wavenumber k ∥ ∼ k ⊥ sin(θ).According to Eq. ( 8), this gives the frequency of linear Alfvén mode, ω ∼ k ⊥ c sin(θ) ∼ 0.1 k ⊥ c, which is not far away from the measurements in Fig. 6.Second, it may correspond to the linear "Doppler shift" of the frequency provided by the passive advection of the small-scale plasma structures by large-scale Alfvén fluctuations, U ⊥,rms .Since U ⊥,rms ∼ 0.1 c in our runs (see Table 3), the corresponding angle-averaged frequency shift, ω ∼ k ⊥ U ⊥,rms / √ 2, is also consistent with the measurements in Fig. 6.
The magnetic-field spectrum associated with the ordinary mode is found from the Faraday law: where the frequency should be substituted from Eq. ( 20).Since the frequency of the ordinary fluctuations is much larger than the frequency of the Alfvén mode, we may average high-frequency ordinary fluctuations independently of the low-frequency Alfvén fluctuations.We may then obtain the electromagnetic spectrum associated with the Alfvén modes by subtracting the spectrum of the ordinary mode (22) from the total spectrum of magnetic fluctuations.Fig. 7, left panel, shows the spectra of electric and magnetic fields in the small-box run II.The right panel shows the spectra of the electric and magnetic fields where the fluctuations associated with the ordinary mode have been removed.The observed spectrum is close to k −3 , which is consistent with the Kraichnan spectrum expected for the kinetic-scale turbulence; see Eqs. ( 11) and ( 12).
Figure 7.The electric and magnetic spectra for the small-box run II.The left panel shows the spectra of the perpendicular components of electric and magnetic fields, as well as the spectrum of Ez.The Ez-spectrum is concentrated at the non-relativistic inertial scale defined in Eq. ( 19) and corresponds to the ordinary mode produced by the initial magnetic perturbations.The right panel shows the kinetic-Alfvén spectra where the fluctuations associated with the ordinary mode have been removed.The spectra indicated by the solid lines are given for the reader's orientation.

CONCLUSIONS
We presented an analytical and numerical study of kinetic-scale Alfvén turbulence in a magnetically dominated ultrarelativistic pair plasma.We assumed that a relatively strong background magnetic field is imposed on a plasma, which is a situation relevant, e.g., for magnetospheres of neutron stars (Arons & Barnard 1986;Gedalin et al. 1998;Nättilä & Beloborodov 2022).This leads to the separation of the inertial and gyro scales and the existence of the kinetic interval of turbulence at scales 1/d rel ≪ k ≪ 1/ρ e .Such kinetic-scale turbulence may be relevant for energy dissipation and particle energization in a turbulent plasma, and it is somewhat analogous to the kinetic-Alfvén or inertial-Alfvén turbulence previously studied in non-relativistic cases.We however demonstrated that the kinetic-scale energy cascade in the ultrarelativistic case is qualitatively different from the non-relativistic counterparts.
First, contrary to non-relativistic kinetic-Alfvén or inertial-Alfvén turbulence, the thermal and inertial effects are necessarily on the same order in the ultrarelativistic kinetic-scale turbulence.Second, the scaling of the energy spectrum is slightly steeper than k −3 , which is different from the kinetic-Alfvén case (the energy spectrum ∼ k −8/3 (e.g., Alexandrova et al. 2009;Boldyrev & Perez 2012;TenBarge & Howes 2012;Boldyrev et al. 2015;Zhou et al. 2023)) and the inertial-Alfvén case (the spectrum is ∼ k −11/3 (e.g., Loureiro & Boldyrev 2018;Milanese et al. 2020)).We have proposed that in the ultrarelativistic case, the energy spectrum is consistent with the Kraichnan spectrum of incompressible 2D turbulence corresponding to the enstrophy cascade, k −3 or k −3 ln −1/3 (k/k 0 ) if the intermittency corrections are taken into account.We noted that the kinetic-scale cascade may be affected by Landau damping which, in general, is not weak in the relativistic case.The spectrum, however, exhibits a near power-law behavior despite being affected by the damping.We conjectured that this might be the consequence of the non-locality of kinetic-scale turbulence.
We performed numerical simulations for two cases, with the box size much larger than the electron inertial scale d rel and the box size moderately larger than the inertial scale.In the first case, the hydrodynamic Alfvénic cascade (kd rel ≪ 1) was well established, while the kinetic-scale turbulence was resolved with about 15 cells per non-relativistic inertial scale d e .In the small-box case, the hydrodynamic interval was shorter, but the kinetic-scale turbulence was resolved significantly better, with 80 cells per d e .The small box also had stronger initial magnetic fluctuations at the kinetic scales.As a result, the "contamination" of the spectrum by the O-mode was stronger, and the measurement of the power-law exponent of the kinetic energy spectrum was less straightforward in the small-box run.In both cases, however, the spectrum of kinetic scale turbulence was consistent with the Kraichnan scaling, from k −3 to k −3 ln −1/3 (k/k 0 ).
It is known that particles get strongly energized in collisionless magnetically-dominated relativistic turbulence (e.g., Zhdankin et al. 2017;Comisso & Sironi 2019;Nättilä & Beloborodov 2021;Vega et al. 2022b).In the case of a strong guide field considered here, the particles are accelerated predominantly in the direction par-allel to the background magnetic field (e.g., Nättilä & Beloborodov 2022), leading to a quasi-one-dimensional particle distribution function.We found that in both cases of the large and small boxes, the resulting particle distribution functions are qualitatively the same; they are well approximated by a universal log-normal distribution.This is in contrast with particle acceleration in weak guide-field turbulence, where energetic particles develop power-law energy distribution functions (e.g., Zhdankin et al. 2019;Wong et al. 2020;Comisso & Sironi 2019;Vega et al. 2022bVega et al. , 2023)) We consider a strongly magnetized, magnetically dominated pair plasma, where the cyclotron frequency is much larger than the plasma frequency and the frequency of the considered wave modes.Here, we imply the relativistic versions of the cyclotron and plasma frequencies that depend on the details of the particle distribution function, see the discussion below.In what follows, we denote by ω pe = 4πn 0 e 2 /m e the non-relativistic electron plasma frequency.We also assume that the electron gyroscale is negligibly small, k 2 ⊥ ρ 2 e ≪ 1.It is convenient to choose the coordinate frame such that k = (k ⊥ , 0, k z ), where z is the direction along the background magnetic field.Under these conditions, the plasma dielectric tensor simplifies to: where the function P (ω, k) depends on the particle distribution function and will be discussed later.In order to find the frequencies and polarizations of the plasma modes, we need to solve the wave equation: which in the matrix form reads Equating the determinant of the matrix to zero, we obtain: Setting the first multiplicative term to zero, one gets the dispersion relation of the electromagnetic extraordinary mode, whose electric-field polarization is normal to both the background magnetic field and the wave vector, E X = (0, E y , 0). (A6) By equating the second term to zero, we obtain In order to specify the function P in this expression, we need to know the particle distribution function.Following (Godfrey et al. 1975;Arons & Barnard 1986;Gedalin et al. 1998), we assume that the particle velocity distribution function is one-dimensional, f (u) = f (u z )δ(u ⊥ ), with the normalization In this expression, we denote u z = v z / 1 − v 2 z /c 2 , where v z is the particle velocity, so that γ 2 = 1 + u 2 z /c 2 .This simplifying assumption is motivated by two independent considerations.First is the fact that in a very strong guide magnetic field such as that relevant, e.g., for pulsar magnetospheres and winds, (e.g., Arons & Barnard 1986;Gedalin et al. 1998), the field-perpendicular components of particle momenta are significantly reduced with respect to their field-parallel ones due to strong synchrotron cooling.Second is the numerical observation that magnetically dominated Alfvénic turbulence (σ ≫ 1) strongly heats plasma particles in the field-parallel rather than field-perpendicular direction (Nättilä & Beloborodov 2022).We also note that in astrophysical applications, plasma can stream along the background magnetic field with relativistic velocity.Our consideration will apply to the plasma rest frame, where we assume that the particle momentum distribution is symmetric with respect to the ±z-directions.
In the considered limit of a very strong large-scale magnetic field and a one-dimensional particle velocity distribution function, one obtains for a pair plasma (Gedalin et al. 1998): The function W in this expression is given by

Figure 1 .
Figure 1.Time evolution of the electromagnetic energy.

Figure 2 .
Figure 2. Visualization of magnetic fluctuations shows similar structures formed in large-and small-box runs.

Figure 3 .
Figure 3. Electron energy distribution functions.The black dashed lines show the log-normal fits with the parameters given in Table2.The thin vertical lines show the ranges of γ where the fits are made.

Figure 4 .
Figure 4.The electron distribution functions.Here, p ∥ and p ⊥ are particle momenta in the directions parallel and perpendicular to the background magnetic field B0.The functions are strongly anisotropic.The anisotropy at very large energies is slightly less pronounced in run II, which is likely a consequence of the smaller box size.

Figure 5 .
Figure5.Spectra of electric and magnetic fluctuations in the large-box run I.The total energy spectrum is slightly steeper than k −3 , but is consistent with the spectrum including a logarithmic intermittency correction, cf.Eqs (11) and (12).

Figure 6 .
Figure 6.Frequencies of the ordinary fluctuations and the potential fluctuations measured in Run II according to Eqs. (A15) and (21).Here, we denote dp ≡ c/ωp, where the relativistic plasma frequency ωp is defined in Eq. (18), and Ω ce,rel ≡ Ωce/⟨γ⟩, where Ωce is the non-relativistic electron cyclotron frequency and ⟨γ⟩ is given in Table2.

Table 1 .
The parameters of the runs.

Table 2 .
Parameters of the particle distribution functions.
, which may indicate different acceleration mechanisms in these two cases.We thank the referee for the valuable comments.This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under award number DE-SC0024362.The work of CV and SB was also partly supported by the NSF grant PHY-2010098 and the Wisconsin Plasma Physics Laboratory (US Department of Energy Grant DE-SC0018266).VR was also partly supported by NASA grant 80NSSC21K1692.Computational resources were provided by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin and by the NASA High-End Computing (HEC)