Spectra of magnetic turbulence in a relativistic plasma

We present a phenomenological and numerical study of strong Alfv\'enic turbulence in a magnetically dominated collisionless relativistic plasma with a strong background magnetic field. In contrast with the non-relativistic case, the energy in such turbulence is contained in magnetic and electric fluctuations. We argue that such turbulence is analogous to turbulence in a strongly magnetized non-relativistic plasma in the regime of broken quasi-neutrality. Our 2D particle-in-cell numerical simulations of turbulence in a relativistic pair plasma find that the spectrum of the total energy has the scaling $k^{-3/2}$, while the difference between the magnetic and electric energies, the so-called residual energy, has the scaling $k^{-2.4}$. The electric and magnetic fluctuations at scale $\ell$ exhibit dynamic alignment with the alignment-angle scaling close to $\cos\phi_\ell\propto \ell^{1/4}$. At scales smaller than the (relativistic) plasma inertial scale, the energy spectrum of relativistic inertial Alfv\'en turbulence steepens to $k^{-3.5}$.


INTRODUCTION
Relativistic plasma turbulence may be present in astrophysical objects, such as jets from active galactic nuclei, pulsar-wind nebulae, magnetospheres of stars and accretion discs. In particular, interest has been attracted to the magnetically-dominated regime, where the magnetic energy exceeds the rest-mass energy of the plasma. It has been discovered that such turbulence leads to efficient energization of plasma particles. This leads not only to thermal plasma heating but also to particles accelerating to ultrarelativistic energies, a process manifested in power-law tails of the particle energy distribution functions (Zhdankin et al. 2017(Zhdankin et al. , 2018a(Zhdankin et al. ,b, 2019(Zhdankin et al. , 2021Comisso & Sironi 2018, 2019Wong et al. 2020;Nättilä & Beloborodov 2021;Demidem et al. 2020;Vega et al. 2021;Chernoglazov et al. 2021;Nättilä & Beloborodov 2022). In this respect, magnetically dominated relativistic turbulence can be considered as a mechanism of particle acceleration complementary to previously studied particle acceleration by collisionless shocks or magnetic-reconnection events (e.g., Marcowith et al. 2016;Guo et al. 2020).
Numerical studies indicate that particle energization crucially depends on the properties of relativistic plasma turbulence itself. Recent numerical and phenomenological consideration suggests that initially ultrarelativistic large-scale (outer-scale) turbulent fluctuations tend to become mildly relativistic in a few turnover times, whereas most of the energy of large-scale fluctuations is converted into thermal energy of plasma particles so that plasma temperature becomes ultrarelativistic. Moreover, the magnetic fluctuations exhibit Kolmogorov-like power-law spectra, bearing some similarity with the nonrelativistic case (e.g., Vega et al. 2021;Zhdankin et al. 2018b;Comisso & Sironi 2019).
Let us assume that a plasma is immersed in a uniform background magnetic field B 0 and denote the typical (rms) strength of magnetic fluctuations as δB 0 . The magnetization parameter related to the magnetic guide field, characterizes the ratio of the magnetic energy to the kinetic energy of plasma particles. In the relativistic case, it may be estimated as σ ∼ B 2 0 /(4πnγmc 2 ), where n is the average number density of plasma particles, m is their mass, and γ is the typical gamma factor of their thermal distribution. In this work, we concentrate on the regime σ ≫ 1 and B 0 ≫ δB 0 . This may correspond to the case when a strong guide field is imposed by external sources (e.g, magnetoshperes), or when a small subregion of a turbulent domain is considered where magnetic fluctuations are much smaller than the mean magnetic field.
Recent numerical and analytical studies of relativistic plasma turbulence using magnetohydrodynamic (MHD) approximation suggested analogies with the non-relativistic case, such as the description in terms of the Elsasser fields, the similarities in the power-law exponents of the spectra of magnetic fluctuations, and the existence of dynamic alignment between magnetic and velocity fields (TenBarge et al. 2021;Chernoglazov et al. 2021). It was also noted that the relative magnitude of the electric field fluctuations in this case exceeds that of the kinetic fluctuations (Vega et al. 2021;Chernoglazov et al. 2021).
In this work, we study turbulence in a magnetically dominated relativistic collisionless plasma. We derive two-fluid equations governing dynamics of Alfvénic fluctuations and compare their predictions to first-principle kinetic particle-in-cell (PIC) simulations. We demonstrate that there is a remarkable mathematical correspondence of the derived equations with those describing dynamics of a non-relativistic plasma, which suggests that a similarity exists between properties of turbulence in relativistic in non-relativistic plasmas. Such a similarity may look rather counter-intuitive from the physical point of view though, since in classical non-relativistic magnetohydrodynamic turbulence, electric fluctuations are negligible while magnetic fluctuations are approximately equipartitioned with fluid motions. In stark contrast with non-relativistic cases, the turbulence discussed here is energetically dominated by fluctuations of electromagnetic field alone.
To resolve this seeming controversy, we note that non-relativistic plasma dynamics may, in fact, also contain a regime where electric field fluctuations are relatively strong as compared to kinetic fluctuations. Such a regime corresponds to broken plasma quaisneutrality, and it is rarely discussed in non-relativistic literature due to its limited applicability to practical systems. However, it turns out to be valuable for understanding magnetically dominated relativistic turbulence. The revealed physical similarity between relativistic and non-relativistic cases allows us to address the total energy spectrum of relativistic turbulence (the sum of electric and magnetic spectra) in both Alfvén and inertial Alfvén intervals, the scale-dependent angular alignment between electric and magnetic fluctuations, and the spectrum of the residual energy (the difference between the magnetic and electric spectra).

ALFVÉN DYNAMICS OF MAGNETICALLY DOMINATED PLASMA
Consider non-relativistic motion of a relativistically hot collisionless plasma with temperarure T α ≫ m α c 2 , where α = {e, i} denotes the types of particles. The momentum equation takes the form (e.g., Mihalas & Mihalas 1999): where ϵ α = n α m α c 2 + u α + p α is the enthalpy density, n α is the particle number density, u α is the internal energy density, and p α is the pressure of the particles. For simplicity, we assume isotropic pressure. For untrarelativistic temperatures of the species, we can approximate ϵ α ≈ u α + p α ≈ 4p α . In addition, one can assume some equation of state, for instance, the adiabatic law of relativistically hot plasma, p α ∝ n 4/3 α . We now assume that the fluctuations of magnetic and electric fields are much weaker than the large-scale uniform magnetic field B 0 = B 0ẑ . We will be interested in Alfvénic plasma fluctuations that are relatively low frequency as compared to cyclotron frequencies, ω ≪ Ω α /γ α , where γ α is the typical gamma factor of particle thermal motion and Ω α is the non-relativistic gyrofrequency. We also impose a self-consistent assumption that the Fourier spectra of the fields are anisotropic in the Fourier space with respect to the background magnetic field B 0 , k ∥ ≪ k ⊥ , and the fluctuations obey the critical balance condition δB/B 0 ∼ k ∥ /k ⊥ ≪ 1 (e.g., Goldreich & Sridhar 1995).
The following derivation is analogous to the procedure developed for the non-relativistic case (e.g., Chen & Boldyrev 2017;Loureiro & Boldyrev 2018;Boldyrev et al. 2021;Milanese et al. 2020). Here, we outline its main steps, the details can be found in the above references. As follows from the momentum equation, to the leading order in the small parameter ωγ α /Ω α the particle motion is the E-cross-B drift, while to the next order it is the polarization drift, where The field-parallel component of the velocity field is expressed through the parallel electric current, n α v α,∥ = J α,∥ /q α . These velocities are substituted into the continuity equation, To the leading order in magnetic, electric, and density fluctuations, one then obtains: In this equation, n 0 is the unperturbed density of each species, δn α = n α − n 0 is the corresponding density perturbation, and ϕ is the electric potential, E = −∇ϕ − 1 c ∂A/∂t. The parallel gradient is taken along the direction of the magnetic field, ∇ ∥ = ∂/∂z − 1 B0 (ẑ ×∇A z )·∇, where the field-perpendicular magnetic perturbation is expressed as δB ⊥ = −ẑ × ∇A z . In order to find the electric current, we turn to the magnetic-field-parallel component of the momentum equation: We now multiply each of the Equations (4) and (5) by n 0 q α and sum over particle species. As a result, Equation (4) turns into the charge conservation law: where ρ = q i δn i + q e δn e is the electric charge density, J ∥ = J e,∥ + J i,∥ is the parallel current, and ϵ 0 = (ϵ i,0 + ϵ e,0 ) /2 is the mean unperturbed enthalpy. To simplify the formulae, we have assumed without loss of generality that q i = |q e | ≡ q. We may also assume that in the relativistic case, the unperturbed enthalpy is the same for both species, ϵ i,0 = ϵ e,0 = ϵ 0 . In the ultrarelativistic limit, ϵ 0 = 4p 0 . The parallel momentum equation (5) will then lead to where δp = δp i − δp e denotes the pressure imbalance. We note (see also Vega et al. 2021) that the terms containing the electric charge density ρ and pressure imbalance δp in Eqs. (6, 7) reflect deviation of plasma fluctuations from quasineutrality, that is, from the condition δn i = δn e . It is easy to see by using the Gauss law, −∇ 2 ϕ = 4πρ, that these terms are relatively small in the case of weak magnetization, that is, when σ ∼ B 2 0 /(8πϵ 0 ) ≪ 1. In the opposite case of strong magnetization that we consider in this work, the deviation from quasineutrality is essential and as a result, both the electric charge and pressure imbalance become dynamically significant. In the case of a non-relativistic electron -proton plasma, we need to replace ϵ 0 → n 0 m i c 2 /2, and the magnetization parameter turns into the so-called quasineutrality parameter, Ω 2 i /ω 2 pi = λ 2 i /ρ 2 i . Here, λ i is the ion Debye scale, ρ i is the gyroscale, and ω pi is the ion plasma frequency. Therefore, Ω 2 i /ω 2 pi is the nonrelativistic analog of relativistic magnetization σ.
In order to close system (6)- (7), we replace the charge density by using the Gauss law, ρ = −(1/4π)∇ 2 ⊥ ϕ, express the parallel electric current as J ∥ = −(c/4π)∇ 2 ⊥ A z , and use the adiabatic law for each particle species to evaluate the pressure gradients, so that Here, we assume that the unperturbed pressure is the same for both species, p 0 = p i,0 = p e,0 . Finally, we introduce the new variables for the scalar and vector Below, we will use only these variables and omit the overtilde sign. Substituting these expressions into Eqs. (6) and (7), we finally obtain the system of equations describing Alfvén dynamics of a relativistic plasma in both magnetohydrodynamic and inertial regimes: where v A = cB 0 / (8πϵ 0 ) is the Alfvén speed, 2 . d = ϵ 0 /(8πn 2 0 q 2 ) is the relativistic inertial length, and the field-parallel gradient has the form The dispersion relation of the linear waves described by Eqs. (9, 10) is 2 In relativistic studies, it is common to use the relativistic Alfvén speed defined asṽ A = v A / 1 + v 2 A /c 2 , so that it is always bounded by the speed of light. In our discussion, we use the quantity v A to make the analogy with the non-relativistic case more transparent. The magnetically dominated case considered in this work corresponds toṽ A ≈ c.
which can be termed the relativistic inertial Alfvén waves.
Except for the very last term in Eq. (10) describing the relativistic pressure contribution, the system of equations (9) and (10) is analogous to the non-relativistic case. The nonrelativistic electron-ion case is recovered by replacing the Alfvén speed and the inertial length by their non-relativistic counterparts by means of the substitutions ϵ 0 → n 0 m i c 2 /2 in the Alfvén velocity and ϵ 0 → 2n 0 m e c 2 in the inertial length. It may also be instructive to compare the dispersion relation (12) with the dispersion relation of non-relativistic inertial kinetic Alfvén waves (e.g., Streltsov & Lotko 1995;Lysak & Lotko 1996), (Boldyrev et al. 2021, Eq. 19) where, similarly, the kinetic correction coming from thermal particle motion enters the numerator, while the inertial correction appears in the denominator. In our relativisic case, these two corrections are necessarily on the same order. We also note that similarly to the previous study (Ten- Barge et al. 2021), at large scales k 2 ⊥ d 2 ≪ 1, Eqs. (9, 10) describe shear Alfvén waves in a relativistically hot plasma and they are mathematically analogous to the equations of reduced magnetohydrodynamics.
Finally, as can be checked directly, Eqs. (9, 10) conserve two quadratic integrals, the energy: and the generalized cross-helicity: 3. NUMERICAL RESULTS In the hydrodynamic range of scales, k 2 ⊥ d 2 ≪ 1, the energy integral becomes: where in the second expression, we have restored the corresponding physical fields. We see that the ratio of the electric to kinetic energy is given by the parameter v 2 A /c 2 = B 2 0 /(8πϵ 0 ) ∼ σ. In both relativistic and non-relativistic cases, when σ ≫ 1, the charge fluctuations are significant and the electric energy dominates the kinetic energy. Therefore, the energy of fluctuations is mostly contained in magnetic and electric fields.
In this limit, the system of equations (9, 10) is mathematically analogous to the equations of non-relativistic reduced magnetohydrodynamics, with the only difference that in the magnetically dominated case, the field ϕ in these equations should be associated with the intensity of electric rather than kinetic fluctuations. Based on this analogy, we may conjecture that the spectrum of the total energy of relativistic magnetically dominated plasma turbulence, the spectrum of its residual energy, and the alignment angle of turbulent fluctuations, should be similar to their reduced MHD counterparts when reinterpreted in terms of magnetic and electric fields as discussed above. Here, we analyze these spectra using kinetic particle-in-cell simulations of a relativistic collisionless plasma.
We numerically study decaying 2D turbulence in a pair plasma where the imposed uniform magnetic field B 0 = B 0ẑ is much stronger than the initial magnetic perturbations. We use the fully relativistic particle-incell code VPIC (Bowers et al. 2008). The fluctuating fields (magnetic, electric, particle velocities) have three vector components, but depend on only two spatial coordinates x and y. Similarly to our previous study (Vega et al. 2021), we denote the root-mean-square of initial magnetic perturbations as δB 0 = ⟨δB 2 (x, t = 0)⟩ 1/2 , and define the two magnetization parameters related to the guide magnetic field and the fluctuations, respectively: In these formulae, n 0 denotes the mean density of each plasma species, w 0 mc 2 is the initial particle enthalpy, with w 0 = K 3 (1/θ 0 )/K 2 (1/θ 0 ). Here K ν is the modified Bessel function of the second kind, and we use the dimensionless initial temperature parameter θ 0 = kT 0 /mc 2 . The particle distribution function is initialized with an isotropic Maxwell-Jüttner distribution, with the temperature parameter θ 0 = 0.3. For such an initialization, w 0 ≈ 1.88. We use a double-periodic square simulation domain with dimensions L x = L y ≈ 2010 d e , where d e is the (nonrelativistic) inertial scale of the plasma particles. The adopted numerical resolution is N x = N y = 16640. The simulations have 100 particles per cell per species. The simulation plane is normal to the uniform magnetic field. The runs are initiated by randomly phased magnetic fluctuations of the Alfvénic type where the wave numbers are chosen in the interval k = {2πn x /L y , 2πn y /L y }, with n x , n y = 1, ..., 8, and χ k are the random phases. The field polarizations correspond to the shear-Alfvén modes (e.g., Lemoine et al. 2016;Demidem et al. 2020),ξ k = k × B 0 /|k × B 0 |, and all the amplitudes δB k are chosen to be the same. Table 1 summarizes the simulations discussed in this paper. As was previously noted in (Vega et al. 2021), the bulk turbulent fluctuations are only mildly relativistic, with average Lorentz factor ⟨γ⟩ ≳ 1, while the particles are strongly relativistic, with average Lorentz factor ⟨γ⟩ ≈ 3.3 − 8.4. This is consistent with the phenomenological discussion of (Vega et al. 2021) that as the initially strong magnetic fluctuations relax, the magnetic energy is mostly converted into heat rather than kinetic energy of collective plasma motion. The Lorentz factors were measured at simulation time t = 16l/c, or about two light crossing times of the simulation box, where l = 2π/k x,y (n = 8) = L x,y /8. By this time, quasisteady states for the distributions of fields have been established. The light crossing time of the simulation box corresponds to a few large-scale dynamical times of turbulence, l/v, where v ≈ 0.3 − 0.4 c is the typical velocity of turbulent fluctuations.
In Figure 1, we present the spectra of magnetic and electric fluctuations as well as the total energy spectrum, W k ⊥ = |B k ⊥ | 2 + |E k ⊥ | 2 . The phase-volume compensated scaling of the energy spectrum in the magnetohydrodynamic interval of scales k ⊥ d e ≪ 1 is close to . Such a spectrum is expected in non-relativistic magnetohydrodynamic turbulence (e.g., Boldyrev 2006;Boldyrev et al. 2009;Mason et al. 2006Mason et al. , 2012Perez et al. 2012;Tobias et al. 2013;Chandran et al. 2015;Chen 2016;Kasper et al. 2021), where the energy is contained in magnetic and kinetic fluid fluctuations. We see that it also holds in relativistic collisionless plasma turbulence dominated by magnetic and electric fields. Our result is also consistent with the recent relativistic MHD studies (TenBarge et al. 2021;Chernoglazov et al. 2021).
At kinetic scales, 1/d 2 e ≪ k 2 ⊥ ≪ 1/ρ 2 e , the inertial Alfvén waves are transformed into ω 2 = 1 3 k 2 z v 2 A /(1 + v 2 A /c 2 ), while the energy integral (13) takes the form: Figure 1. The magnetic and electric spectra, and the total electromagnetic energy spectrum, for two different magnetizations. The total energy spectrum approaches a k −3/2 power-law at k ⊥ de ≪ 1, and a k −3.5 power-law at k ⊥ de ≫ 1.
In this asymptotic region, dimensional estimates give for turbulent fluctuations at field-perpendicular scale ℓ: A z,ℓ ∼ ϕ ℓ , and the nonlinear interaction time is estimated as τ ℓ ∼ ℓ 2 /ϕ ℓ . The condition of constant flux of the conserved quantity (18) then reads ϕ 2 ℓ /(ℓ 4 τ ℓ ) = const, which gives for the scaling of fluctuations ϕ ℓ ∝ ℓ 2 and for the electromagnetic energy spectrum of relativistic inertial Alfvén waves W k ⊥ 2πk ⊥ ∝ k −3 ⊥ . The nonlinear interaction time for such modes turns out to be independent of scale, which implies that for the critically balanced cascade, ω ∝ 1/τ ℓ , we have k z ∝ const, that is, the cascade proceeds in the field-perpendicular direction.
We note that such a turbulent cascade is somewhat analogous to the cascade of enstrophy in incompressible non-relativistic two-dimensional fluid turbulence. Such a cascade is however only marginally local, so it depends on the conditions at the low-k boundary of the inertial interval. Numerical simulations and experiments typically reveal the spectra steeper than the −3 in this case, unless the formation of large-scale structures is controlled (suppressed) by some forcing mechanism and/or the inertial interval is sufficiently large (e.g., Boffetta & Ecke 2012). Our numerical simulations found the spectral scaling close to k −3.5 ⊥ . In addition to the limited separation between the particle inertial and gyroradius scales and possible spatial intermittency, in our case the steepening may also be related to Landau damping at kinetic scales (e.g., Nättilä & Beloborodov 2022) as the phase velocity of the inertial-Alfvén waves, c/ √ 3, is smaller than the thermal speed of particles, ∼ c.
In non-relativistic magnetohydrodynamic turbulence, the magnetic energy is known to exceed the energy of kinetic fluctuations. Phenomenological and numerical considerations demonstrated that the difference between the magnetic and kinetic energies, the so-called residual energy, is positive and has the spectrum close to −2 (e.g., Boldyrev et al. 2011Boldyrev et al. , 2012Chen et al. 2013). In the relativistic magnetically dominated case, we may introduce the analog of residual energy as the difference between the magnetic and electric energies, This quantity is measured in Figure 2. In order to compensate for overall energy decline in decaying turbulence we have normalized the residual energy by the total energy of fluctuations and then averaged over several data cubes. A powerlaw spectrum is indeed observed, however, the scaling is slightly steeper than in its non-relativistic counterpart, more consistent with R k ⊥ 2πk ⊥ ∝ k −2.4 ⊥ . Finally, a characteristic feature of non-relativistic magnetic plasma turbulence in the presence of a strong guide field, is the dynamic alignment between the shear- Alfvén magnetic and velocity fluctuations. Such fluctuations become spontaneously aligned in a turbulent cascade, with the alignment angle scaling as sin θ ℓ ∝ ℓ 1/4 , (e.g., Boldyrev 2006;Podesta et al. 2009;Chen et al. 2011;Perez et al. 2012). Such an alignment progressively reduces the strength of nonlinear interaction at small scales, which arguably explains the shallowerthan-Kolmogorov spectrum of turbulent energy, k −3/2 ⊥ . In the case of magnetically-dominated turbulence, one may similarly expect a scale-dependent dynamic alignment between the electric and magnetic field fluctuations, where ℓ is the scale of the fluctuations in the fieldperpendicular plane, e.g., δB ≡ B(x + ℓ) − B(x) where ℓ ⊥ B 0 . In order to see whether a similar alignment exists in our numerical simulations of collisionless relativistic turbulence, we plot the cosine of the angle ϕ ℓ vs scale ℓ in Figure 3. We observe a scaling close to that of the non-relativistic case, and also to the relativistic MHD case . We however notice that the scaling varies slightly with the plasma magnetization parameter σ 0 ; it is slightly shallower in the case of stronger magnetization, possibly reflecting a shorter inertial range due to larger relativistic inertial scale.

DISCUSSION
We have presented a numerical and phenomenological study of relativistic turbulence in a collisionless magnetically dominated plasma. We proposed that dynamic equations (9, 10) provide a useful phenomenological model for magnetically dominated relativistic plasma turbulence. Using the 2D particle-in-cell simulations, we demonstrated that in such turbulence, the energy is contained mostly in collective magnetic and electric fluctuations. In the magnehydrodynamic range of scales, k 2 ⊥ d 2 e ≪ 1, the total energy (magnetic plus electric) exhibits the spectrum close to that of Alfvénic turbulence, W k ⊥ 2πk ⊥ ∝ k −3/2 ⊥ . In the kinetic range, 1/d 2 e ≪ k 2 ⊥ ≪ 1/ρ 2 e , the turbulence is governed by relativistic inertial Alfvén modes, and the spectrum steepens to approximately k −3.5 ⊥ .
We further established that in relativistic magnetically dominated turbulence, there is excess of magnetic energy over electric one, which becomes progressively smaller at smaller scales. We propose that this phenomenon is analogous to the generation of the so-called residual energy, that is, the excess of magnetic over kinetic energy known in the non-relativistic quasineutral case. The measured spectrum of the residual energy is close to R k ⊥ 2πk ⊥ ∝ k −2.4 ⊥ , which is slightly steeper than its non-relativistic counterpart, indicating an interesting difference with the non-relativistic case.
An additional intriguing similarity of relativistic turbulence with the nonrelativistic case is manifested in the presence of the so-called scale-dependent dynamic alignment between magnetic and electric fluctuations.
We have found that such fluctuations become progressively more orthogonal to each other at smaller scales. The cosine of the angle between the fluctuations was found to scale close to cos ϕ ℓ ∝ ℓ 0.25 . Based on the analogy with the nonrelativistic case, we proposed that such an alignment reduces the strength of nonlinear interactions in the relativistic dynamics, thus explaining the observed −3/2 scaling of the energy spectrum. We also note that the run with a stronger magnetization resulted in a slightly shallower scaling of the alignment angle. This may be related to a shorter inertial interval of turbulence due to a larger relativistic inertial scale of thermal particles.
Finally, we note that our consideration may provide a useful framework for phenomenological studies of relativistic magnetically dominated plasma turbulence in a presence of a strong background magnetic field. In particular, it may be relevant for the analysis of particle heating and acceleration mediated by such turbulence, as well as the effects of magnetic reconnection, struture formation, and intermittency generated by a turbulent cascade.