Weak Alfvénic turbulence in relativistic plasmas. Part 1. Dynamical equations and basic dynamics of interacting resonant triads

Abstract

Alfvén wave collisions are the primary building blocks of the non-relativistic turbulence that permeates the heliosphere and low- to moderate-energy astrophysical systems. However, many astrophysical systems such as gamma-ray bursts, pulsar and magnetar magnetospheres and active galactic nuclei have relativistic flows or energy densities. To better understand these high-energy systems, we derive reduced relativistic magnetohydrodynamics equations and employ them to examine weak Alfvénic turbulence, dominated by three-wave interactions, in reduced relativistic magnetohydrodynamics, including the force-free, infinitely magnetized limit. We compare both numerical and analytical solutions to demonstrate that many of the findings from non-relativistic weak turbulence are retained in relativistic systems. But, an important distinction in the relativistic limit is the inapplicability of a formally incompressible limit, i.e. there exists finite coupling to the compressible fast mode regardless of the strength of the magnetic field. Since fast modes can propagate across field lines, this mechanism provides a route for energy to escape strongly magnetized systems, e.g. magnetar magnetospheres. However, we find that the fast-Alfvén coupling is diminished in the limit of oblique propagation.

1. Introduction

Turbulence provides the transport of mass, momentum and energy in a wide range of plasmas throughout the universe, from the intracluster medium and magnetar magnetospheres to the solar wind and laboratory fusion confinement experiments. In space and astrophysical plasmas, turbulence plays the important role of transferring large-scale motions, often driven by violent processes and instabilities, to small scales at which damping, dissipation and plasma heating can occur. This cascade of turbulent energy is governed by nonlinear interactions that occur within the plasma, and the process is well studied in the non-relativistic (Newtonian) limit. However, many astrophysical plasmas of interest are relativistic and often magnetically dominated ($b^2 \gg h$, where $b^2 = b_\mu b^\mu = B^2/\gamma ^2 + (\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {v})^2$ is the magnetic energy density, $h$ is the enthalpy density and $\gamma$ is the Lorentz factor). How these energetic plasmas are heated is fundamental for interpreting the electromagnetic radiation we observe at Earth.

We will begin our study of relativistic, magnetically dominated turbulence by reviewing and building upon results from Newtonian turbulence. The magnetic fields that universally permeate plasmas imply that Alfvénic fluctuations (Alfvén 1942) will govern the hierarchy of turbulent fluctuations rather than the eddies that compose hydrodynamic turbulence. The shear (or transverse) Alfvén wave has the property that the fluid motions corresponding to it are entirely transverse to the background magnetic field, with no compressional component. Based on these ideas, Iroshnikov (1963) and Kraichnan (1965) (IK) employ incompressible magnetohydrodynamics (MHD) to propose that the nonlinear interactions in plasma turbulence are composed of counter-propagating and overlapping Alfvén waves or wave packets.

The fundamental assumption underlying the IK picture of turbulence is that the so-called $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity is the dominant nonlinear term in the plasma. In terms of fluid, MHD equations, this term appears as $\boldsymbol {v}_\perp \boldsymbol {\cdot } \boldsymbol {\nabla } \star$, where $\boldsymbol {v}_\perp \simeq \boldsymbol {E} \times \boldsymbol {B}_0/B_0^2$ is the dominant drift velocity in the plasma, $\boldsymbol {B}_0$ is a mean magnetic field, $\star$ indicates either $\boldsymbol {v}$ or $\boldsymbol {B}$ and perpendicular, $\perp$, indicates perpendicular to $\boldsymbol {B}_0$. The form of this term makes clear immediately that the nonlinearity requires fluctuations in the plane perpendicular to the mean magnetic field, i.e. $k_\perp \ne 0$. This term also requires that the perpendicular wavevectors of interacting Alfvén waves are not collinear, but this point is not obvious from the simplified expression above and is explored further in § 2.1. Although there are cases in which the $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity is not dominant, e.g. parametric instabilities such as the decay, modulation and beat instabilities, we will adopt the assumption that the $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity is the dominant nonlinear term. Additionally, the above discussion is based on a fluid description of plasma; however, most high-energy astrophysical systems are in the weakly collisional limit ($l \ll \lambda _{mfp}$, where $l$ is an intermediate turbulence scale and $\lambda _{mfp}$ is the collisional mean free path), formally requiring a kinetic description. Fortunately, kinetic Newtonian turbulence retains many of the same basic properties as MHD turbulence for scales larger than ion kinetic scales, e.g. the ion inertial length and gyroradius. Importantly, the $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity is the dominant nonlinearity in kinetic Newtonian turbulence, even at turbulence scales below ion kinetic scales (Schekochihin et al. 2009).

The IK theory of Alfvénic turbulence assumes that the turbulent cascade proceeds isotropically across scales, i.e. $k_\parallel \sim k_\perp$ at all scales. Montgomery & Turner (1981) and Shebalin, Matthaeus & Montgomery (1983) revisit the IK theory under the assumption that the turbulence is weak, $\omega _{NL} \ll \omega$, where $\omega _{NL} \sim k_\perp \delta v_\perp$ is the nonlinear frequency and $\omega \sim k_\parallel v_A$ is the linear, Alfvén frequency. Under this weak assumption, they demonstrate that the three-wave interaction is the dominant nonlinear interaction of the $\boldsymbol {E} \times \boldsymbol {B}$ nonlinearity. Importantly, this three-wave interaction involves a frequency and wavevector matching condition for two counter-propagating Alfvén waves colliding to produce a $k_\parallel = 0$ mode, leading to an anisotropic cascade of energy which fixes $k_\parallel$ at the outer scale. Based on this interaction, Montgomery & Turner (1981), Montgomery (1982) and Higdon (1984) develop a theory of highly anisotropic incompressible MHD turbulence consisting of two-dimensional velocity and magnetic field fluctuations in the plane perpendicular to the mean magnetic field¹. To describe this system, they employ the reduced MHD (RMHD) equations introduced by Kadomtsev & Pogutse (1974) and Strauss (1976) for systems in which there is a strong mean magnetic field, leading to the ordering assumptions of RMHD: $\delta B_\perp / B_0 \sim \delta v_\perp / v_A \sim k_\parallel / k_\perp \sim \epsilon \ll 1$. Note that the assumptions of RMHD necessarily imply that the Alfvén and fast modes are well separated in frequency since $\omega _{A} \propto k_\parallel \ll \omega _F \propto k \sim k_\perp$. This ordering implies that the fast mode is ordered out of the RMHD system. All of the following discussion regarding the development of MHD turbulence is predicated on this assumption that the fast wave is well separated from the Alfvén mode and is therefore ignorable. To develop a theory of relativistic, magnetically dominated astrophysical plasmas, we revisit the concept of RMHD in the relativistic limit in § 3.2.

Sridhar & Goldreich (1994) take a different approach in expanding upon the work of IK by first pointing out that IK is a theory of weak turbulence, and then loosening the isotropy assumption of IK by building anisotropy into weak turbulence theory. However, Sridhar & Goldreich (1994) argue that the three-wave interaction is empty because it involves a mode with $\omega = 0$, which cannot be a linear wave mode with finite power. Therefore, they invoke a four-wave interaction ($\mathrm {A}+\mathrm {A} \rightarrow \mathrm {A} + \mathrm {A}$) which maintains $k_\parallel = 0$, and as the weak cascade proceeds to smaller scales, the strength of nonlinear interactions increases, eventually reaching a state of strong turbulence. Goldreich & Sridhar (1995) carry the weak turbulence theory developed in Sridhar & Goldreich (1994) into the strong limit. In the strong limit, they argue that the resonance condition for interaction is broadened, and due to the broadening, a parallel cascade develops. Further, they argue that the parallel cascade leads the turbulence towards a state of critical balance in which $\chi = \omega _{NL} / \omega \simeq 1$, where $\chi$ is the nonlinearity parameter and $\chi \ll 1$ corresponds to weak turbulence. Critical balance is a condition in which the nonlinear frequency or cascade time is balanced by the linear time in the system. The critically balanced turbulence cascade is predicted to have an energy spectrum $E \sim k_\perp ^{-5/3}$ and a spectral anisotropy of the form $k_\parallel \sim k_\perp ^{2/3}$, thus developing a scale-dependent anisotropy. Note that Alfvénic turbulence, whether weak or strong, always leads to a case at small scales in which $k_\parallel \ll k_\perp$ and $\delta B \ll B_0$ regardless of the isotropy or amplitude of fluctuations at the outer scale².

Montgomery & Matthaeus (1995) and Ng & Bhattacharjee (1996) concurrently claim that Sridhar & Goldreich (1994) are incorrect to claim that three-wave interactions are empty, because $k_\parallel = 0$ modes are valid nonlinear fluctuations. Ng & Bhattacharjee (1997) further employ perturbation theory to demonstrate explicitly that the three-wave interaction is non-empty, $k_\parallel = 0$ modes do develop and the three-wave interaction dominates over the four-wave one. They also present the energy spectrum of weak turbulence, $E \sim k_\perp ^{-2}$. Admitting their mistake in omitting the three-wave interaction, Goldreich & Sridhar (1997) reformulate their weak turbulence theory to include three-wave interactions. The basic prediction for the spectral anisotropy (no parallel cascade), energy spectrum $E \sim k_\perp ^{-2}$ and the strengthening of the cascade as it proceeds to smaller scales remain unchanged. Galtier et al. (2000) provide a rigorous derivation of the $k_\perp ^{-2}$ scaling result based on the wave kinetic approach (Zakharov, L'Vov & Falkovich 1992), and Perez & Boldyrev (2008) numerically verify the derivation and examine the weak to strong turbulence transition. Lithwick & Goldreich (2003) present further extensions of weak turbulence theory.

More recently, a series of papers were written examining the building blocks of weak turbulence through the interaction of waves analytically (Howes & Nielson 2013), numerically (Nielson, Howes & Dorland 2013) and in an experiment (Drake et al. 2013; Howes et al. 2013) conducted at the Large Plasma Device (Gekelman et al. 1991). This series of papers focuses on a heuristic, analytical solution beginning with a collision of counter-propagating Alfvén waves at first order in the fluctuation amplitude and following their evolution order-by-order. At second order, the three-wave interaction involving the nonlinear, $k_\parallel = 0$ mode is found to be dominant. This mode does not involve a secular exchange of energy and oscillates with frequency $\omega = 2\omega _0$, where $\omega _0 = k_\parallel v_A$ is the frequency of the incident Alfvén waves. At third order, the incident Alfvén waves interact with the $k_\parallel = 0$ mode, and the $k_\parallel = 0$ mode shears the Alfvén waves in the perpendicular plane, providing a secular exchange of energy to smaller perpendicular scales but fixed $k_\parallel$ scale.

The theory of critically balanced strong turbulence has also been further refined since Goldreich & Sridhar (1995). Boldyrev (2005,  2006) note that the vector nature of the nonlinearity leads to a state called dynamic alignment wherein the velocity and magnetic field fluctuations align themselves as the cascade proceeds to smaller scales. This alignment leads to the formation of thin current sheets at small scales, and these current sheets can eventually disrupt (Loureiro & Boldyrev 2017; Mallet, Schekochihin & Chandran 2017; Comisso et al. 2018; Dong et al. 2018). Intermittency has also been built into the theory of critically balanced and dynamically aligned turbulence (Chandran, Schekochihin & Mallet 2015), leading to the theory of refined critical balance (Mallet, Schekochihin & Chandran 2015). For a more complete, contemporary (but biased, according to the authors) review of the current state of MHD turbulence including weak, strong and imbalanced turbulence, see Schekochihin (2020).

Although fundamentally important for understanding energy dissipation and astrophysical observations, relativistic turbulence has received much less attention. Thompson & Blaes (1998) and Troischt & Thompson (2004) examine weak turbulence in the magnetically dominated, relativistic regime. Thompson & Blaes (1998) follow Sridhar & Goldreich (1994) and Goldreich & Sridhar (1995) in arguing that the four-wave interaction is the dominant Alfvénic interaction because the $k_\parallel = 0$ mode is not a linear mode of the system. However, they note that in the extreme relativistic regime, one cannot assume that the fast and Alfvén modes are well separated, i.e. the intuition gained from incompressible MHD no longer holds. Therefore, they argue that the dominant three-wave interactions are of the form $\mathrm {A}+\mathrm {A} \rightarrow \mathrm {F}$, $\mathrm {A}+\mathrm {F} \rightarrow \mathrm {A}$ or $\mathrm {F}+\mathrm {F} \rightarrow \mathrm {F}$. Heyl & Hernquist (1999) extend the work of Thompson & Blaes (1998) to include quantum electrodynamic effects making the same assumptions regarding three- versus four-wave interactions. Li, Zrake & Beloborodov (2019) return to the problem of weak, magnetically dominated turbulence following the work of Thompson & Blaes (1998), again assuming that the dominant Alfvénic interactions are the four-wave coupling or the three-wave couplings involving fast modes, as listed above. Li et al. (2019) follow the analytical discussion with relativistic MHD simulations in the force-free limit considering both weak and strong turbulence limits. They focus on the nonlinear, turbulent conversion of Alfvén to fast mode energy as a possible mechanism to release energy from magnetar magnetospheres, since Alfvén waves propagate along field lines and remain trapped in the magnetosphere, while fast modes can propagate across fields lines, thereby escaping confinement and releasing energy. A variety of other recent papers have explored various aspects of relativistic turbulence theory and simulations, finding broadly similar results to those of Newtonian MHD turbulence, except they consistently find a small portion ($\lesssim$10 %–15 %) of the initial Alfvénic energy leaks into the fast mode branch (Cho 2005; Takamoto & Lazarian 2016,  2017; Li et al. 2019).

This work is organized as a sequence of papers. In this paper (Part 1), we derive a set of relativistic RMHD (RRMHD) equations which have the same form and properties as their Newtonian counterpart, including the wave kinetic equation governing spectral evolution. We then employ an approach similar to that of Howes & Nielson (2013) to obtain an analytical solution for three-wave interactions in relativistic systems. We emphasize that although we make a connection with the wave kinetic equation for the relativistic system, the primary analysis is intended to be heuristic to highlight the role of the three-wave interactions and other similarities to Newtonian turbulence rather than a formal weak turbulence theory. Our approach, wherein we obtain the solution order-by-order, is well suited to comparison with numerical simulations (Part 2) and complements the variational approach of Thompson & Blaes (1998). We build upon the wisdom gained from Newtonian turbulence and begin by outlining and reviewing some of the salient properties for an astrophysical audience of both incompressible and compressible MHD turbulence in § 2. In § 3, we derive the relativistic Elsasser equations in the reduced, relativistic limit and discuss their connection with weak Newtonian turbulence. We then derive through third order the weak Alfvénic turbulence solutions in § 4. In § 5, we compare our solutions with direct numerical simulations and consider the role of fast waves in both relativistically hot and magnetically dominated turbulence. Finally, we summarize the results in § 6. More in-depth numerical simulation studies of weak, relativistic turbulence are presented in Part 2 (Ripperda et al. 2021).

2. General properties of non-relativistic turbulence

2.1. Incompressible MHD

Before exploring relativistic Alfvénic turbulence, it is important to first review some of the fundamental knowledge learned from Newtonian turbulence to provide a framework for discussing the relativistic limit. This discussion is far from exhaustive, because incompressible MHD turbulence is an exceptionally broad and deep field.

We will begin our discussion by presenting some of the basic properties of incompressible turbulence. Incompressible MHD is the basis from which Newtonian turbulence theories have been derived. Incompressible MHD assumes $v \ll c_s$, where $c_s$ is the sound speed. In other words, compressive fluctuations are carried away from the source at essentially infinite speed, an assumption that is not applicable for relativistic systems. This assumption implies $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} = 0$, leading to the incompressible MHD equations

(2.1)\begin{gather} \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} ={-}\boldsymbol{\nabla} P/\rho_0 + \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B}/\rho_0, \end{gather}

(2.2)\begin{gather}\frac{\partial \boldsymbol{B}}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B} = \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}, \end{gather}

(2.3)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} = 0, \end{gather}

where $P$ is the total (thermal plus magnetic) pressure and $\rho _0$ is the equilibrium mass density.

For turbulence analysis, the incompressible MHD equations are typically cast in Elsasser form (Elsasser 1950) by adding and subtracting the momentum and induction equations to arrive at

(2.4)\begin{gather} \frac{\partial \boldsymbol{z}^\pm}{\partial t} \mp \boldsymbol{v}_A \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{z}^\pm{=}{-}\boldsymbol{z}^\mp \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{z}^\pm{-}\boldsymbol{\nabla} P / \rho_0, \end{gather}

(2.5)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{z}^\pm{=} 0, \end{gather}

where the magnetic field has been separated into equilibrium and fluctuating parts, $\boldsymbol {B} = B_0 \hat {\boldsymbol {z}} + \delta \boldsymbol {B}$, $\boldsymbol {v}_A = \boldsymbol {B}_0 / \sqrt {4 {\rm \pi}\rho _0}$ and $\boldsymbol {z}^\pm = \boldsymbol {v} \pm \delta \boldsymbol {B} / \sqrt {4 {\rm \pi}\rho _0}$ are the Elsasser fields. Because the Elsasser fields are divergenceless, this is a closed system of equations (Montgomery 1982). To obtain an equation for the pressure, one can simply take the divergence of (2.4) to find

(2.6)\begin{equation} \nabla^2 P/\rho_0 ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \left( \boldsymbol{z}^\mp \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{z}^\pm\right). \end{equation}

The Elsasser equations, (2.4), are the progenitor equations to describe Alfvénic turbulence, and one can discern many important points about turbulence simply from the form of the equations. First, the terms on the left-hand side are linear, while those on the right are nonlinear. Therefore, linearizing the equations is as simple as setting the right-hand side to zero. By linearizing, one can immediately see that the system supports two propagating linear wave modes with $\omega = \pm k_\parallel v_A$: (i) Alfvén waves with $\boldsymbol {z}^\pm$ polarized in the $\hat {\boldsymbol {z}} \times \widehat {\boldsymbol {k}}^\pm$ direction and (ii) pseudo-Alfvén waves, the incompressible limit of magnetosonic slow modes, with polarization $\widehat {\boldsymbol {k}}^\pm \times (\hat {\boldsymbol {z}} \times \widehat {\boldsymbol {k}}^\pm )$, where $\widehat {\boldsymbol {k}}^\pm$ corresponds to the unit vector along the wavevector of $\boldsymbol {z}^\pm$. We adopt the convention that $\omega \ge 0$, and the sign of $k_\parallel$ determines the propagation direction. Note that the incompressible assumption orders the magnetosonic fast mode out of the system by imposing $c_s \rightarrow \infty$. Thus, one can interpret the Elsasser field $\boldsymbol {z}^+$ ($\boldsymbol {z}^-$) as describing the evolution of Alfvén or pseudo-Alfvén waves propagating down (up) the equilibrium magnetic field.

Next, we turn our focus to the primary nonlinear term, $\boldsymbol {z}^\mp \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {z}^\pm$. The most immediate point one can see from the form of this term is that the nonlinearity only survives if both $\boldsymbol {z}^+$ and $\boldsymbol {z}^-$ are non-zero, i.e. the nonlinearity is one that does not involve self-interaction of an Alfvén wave with itself but its interaction with an oppositely propagating Alfvén wave. If either Elsasser field is zero, then the opposite Elsasser variable is an exact solution. For instance, if $\boldsymbol {z}^- = 0$, then $\boldsymbol {z}^+(x,y,z+v_A t)$ is an exact, nonlinear solution representing an arbitrary-amplitude Alfvén or pseudo-Alfvén wave propagating in the $-\hat {\boldsymbol {z}}$ direction.

Physically, the counter-propagating waves shear one another when they interact through this nonlinear term, and the shearing leads to a transfer of energy to smaller scales. For simplicity, let us focus on the case of counter-propagating Alfvén waves and examine the nonlinear term in more detail. In Fourier space, the nonlinear interaction term for a $\boldsymbol {z}^+$ wave distorted by a $\boldsymbol {z}^-$ wave is

(2.7)\begin{equation} \boldsymbol{z}^- \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{z}^+{\propto} \left[\left(\hat{\boldsymbol{z}} \times \widehat{\boldsymbol{k}}^- \right) \boldsymbol{\cdot}\widehat{\boldsymbol{k}}^+\right] \left(\hat{\boldsymbol{z}} \times \widehat{\boldsymbol{k}}^+ \right) = \left[\hat{\boldsymbol{z}} \boldsymbol{\cdot} \left( \widehat{\boldsymbol{k}}^-{\times} \widehat{\boldsymbol{k}}^+ \right)\right] \left(\hat{\boldsymbol{z}} \times \widehat{\boldsymbol{k}}^+ \right), \end{equation}

where we have used the polarization properties of the Alfvén waves to write $\boldsymbol {z}^\pm \propto \hat {\boldsymbol {z}} \times \widehat {\boldsymbol {k}}^\pm$. Therefore, for the nonlinear interaction of counter-propagating Alfvén waves to be non-zero, $\hat {\boldsymbol {z}} \boldsymbol {\cdot } (\widehat {\boldsymbol {k}}^- \times \widehat {\boldsymbol {k}}^+ ) \ne 0$, i.e. variation is required in both directions perpendicular to the equilibrium field. An alternative way to state this requirement is that the waves must be polarized with respect to one another in the plane perpendicular to the equilibrium field so that the perpendicular wavevector components are not collinear.

From examination of the linear and nonlinear terms of the Elsasser equation, we have gleaned three crucial facts for turbulence: (i) the system supports two linear waves modes, both of which require $k_\parallel \ne 0$ to propagate; (ii) the nonlinearity requires counter-propagating fluctuations; and (iii) the nonlinearity requires that the fluctuations be relatively polarized with one another in the plane perpendicular to the equilibrium field. Thus, to fully capture the physics of the turbulent cascade one requires all three dimensions (Tronko, Nazarenko & Galtier 2013; Howes 2015). This requirement to retain all three dimensions to capture Alfvénic turbulence persists in both full MHD and kinetic plasmas (Schekochihin et al. 2009; Howes 2015).

Finally, we highlight one other important fact about turbulence, which can be seen by examining the Elsasser energy equation, obtained by taking the dot product of $\boldsymbol {z}^\pm$ with (2.4) and integrating over all of space. Assuming periodic boundary conditions or that the fields vanish at infinity, one obtains

(2.8)\begin{equation} \int \textrm{d}^3 \boldsymbol{r} \frac{\textrm{d} |\boldsymbol{z}^\pm|^2}{\textrm{d} t} = 0. \end{equation}

Equation (2.8) implies that there is no exchange of energy between the upward- and downward-propagating fluctuations, even during nonlinear interactions (Maron & Goldreich 2001; Schekochihin et al. 2009). The collisions of counter-propagating fluctuations are therefore elastic: one wave packet can scatter another, but the individual energies of the $z^+$ and $z^-$ fluctuations do not change.

2.2. The connection to RMHD

As noted in § 1, Alfvénic turbulence proceeds in an anistropic fashion as it cascades to smaller scales, eventually leading to a state in which $k_\parallel \ll k_\perp$ (and $\delta B / B_0 \ll 1$) regardless of the initial isotropy of the plasma at the outer scale. This scale-by-scale anisotropic turbulence cascade has been well observed in simulations (Cho & Vishniac 2000; Maron & Goldreich 2001; Chen et al. 2012) and in situ solar wind observations (Horbury, Forman & Oughton 2008; Wicks et al. 2010; Chen et al. 2012; Chen 2016). Thus, it would be advantageous to consider an ordering framework that leverages this fact. Fortunately, the minimal ordering assumptions for RMHD are anisotropic fluctuations ($k_\parallel \sim \epsilon k_\perp$), small-amplitude fluctuations relative to the background (e.g. $\delta B \sim \epsilon B_0$) and characteristic time scales $\omega \sim k_\parallel v_A$, where $\epsilon \ll 1$. Note that early derivations of RMHD (Kadomtsev & Pogutse 1974; Strauss 1976) further assumed strong magnetization, implying plasma $\beta \ll 1$. However, more recent derivations (Schekochihin & Cowley 2007; Schekochihin et al. 2009) have demonstrated that the RMHD equations are valid for arbitrary plasma $\beta$, assuming a homogeneous background. Conveniently, the equations of RMHD are essentially identical to (2.4), with the exception that the gradients in the nonlinear terms reduce to gradients perpendicular to the equilibrium magnetic field.

Much like the incompressible MHD equations, the fast wave is ordered out of RMHD. In the solar wind, this is well justified and supported by observations, which show that fast modes generally compose a small fraction of the solar wind (Tu & Marsch 1994; Howes et al. 2012; Klein et al. 2012). The RMHD equations have gained prominence as the preferred set of equations for describing Alfvénic turbulence because they have a few favourable properties compared with incompressible MHD, despite their close similarity. First, the RMHD equations are a rigorous set of equations for describing collisional or collisionless Alfvénic physics at scales large compared to the ion gyroradius, $k_\perp \rho _i \ll 1$ (Schekochihin et al. 2009). Second, in the anisotropic limit of RMHD, the plus and minus Alfvén mode, plus and minus slow (or pseudo-Alfvén) mode and the lone entropy mode cascades decouple. In other words, there are five independent cascade channels in RMHD, and the energy in each channel is independently conserved. Thus, the five channels do not exchange energy with one another, analogous to the two independent channels described above for incompressible MHD. The Alfvénic cascade is described by the perpendicular vector components of (2.4), the slow mode cascade by the parallel component and the entropy mode cascade by the pressure balance condition (Schekochihin et al. 2009). Formally, the RMHD equations are only valid for the slow and entropy modes in the collisional limit, because these modes are subject to collisionless wave–particle interactions via Landau (Landau 1946) or Barnes (Barnes 1966) damping, even for scales $k_\perp \rho _i \ll 1$. To describe these modes in the collisionless limit, one can instead use kinetic RMHD, which evolves the perpendicular dynamics using conventional RMHD and the parallel dynamics using the reduced ion kinetic equation (Schekochihin et al. 2009). Third, in the anisotropic limit of RMHD, the Alfvén waves are not affected by the slow or entropy modes, and the slow and entropy modes do not cascade on their own. The slow and entropy modes behave as passive scalars, which can only be cascaded to small scales by interacting with Alfvén waves. This fact also implies that slow and entropy fluctuations will not be generated in situ in a purely Alfvénic RMHD turbulence cascade³.

Finally, before moving forward to discuss compressible MHD turbulence, we note that RMHD is equally valid for describing weak and strong (critically balanced) turbulence (Galtier & Chandran 2006). Comparing the strength of the nonlinear and linear terms in RMHD, we find

(2.9)\begin{equation} \chi = \frac{\boldsymbol{z}^\mp\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{z}^\pm}{\boldsymbol{v}_A \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{z}^\pm} \sim \frac{z^\mp k_\perp}{v_A k_\parallel} \sim \frac{k_\perp \delta B }{k_\parallel B_0}. \end{equation}

The final expression is the ratio of asymptotically ordered quantities in RMHD and is therefore unordered. In other words, RMHD can equally well describe weak ($\chi \ll 1$) and critically balanced ($\chi \sim 1$) turbulence.

As a brief aside, it is worthwhile at this point to provide a physical interpretation of critical balance, and justify why we have chosen to ignore the case $\chi > 1$. Physically, critical balance amounts to equating the linear, propagation time, $\tau _A = l_\parallel / v_A$, with the nonlinear (or ‘eddy turnover’) time, $\tau _{NL} \sim l_\perp / u_\perp$. The case $\tau _A \ll \tau _{NL}$ ($\chi \ll 1$) is the weak turbulence case, which will eventually, at sufficiently small scales, transition to $\tau _A \sim \tau _{NL}$ because the weak turbulence cascade does not produce a cascade in the direction parallel to the equilibrium field: $l_\parallel$ is fixed at the outer scale. Thus, in weak turbulence $\tau _A$ remains fixed while $\tau _{NL}$ will decrease with scale. The case $\tau _A \gg \tau _{NL}$ corresponds essentially to creating two-dimensional structures perpendicular $\boldsymbol {B}_0$ and is not sustainable at small scales separated from any external forcing. A fluctuation can only remain two-dimensional if it is causally connected, but for a system with an equilibrium $B_0$, the maximum parallel length over which a fluctuation can be coherent is $l_\parallel \sim \tau _{NL} v_A$, i.e. $l_\parallel / v_A \sim \tau _A \sim \tau _{NL}$. Thus, if a system is driven such that $\tau _A \gg \tau _{NL}$, it will rapidly relax back to critical balance, $\tau _A \sim \tau _{NL}$, at small scales. It is also worth noting that the critical balance condition, or predictions that follow from critical balance, has been observed in simulations (Cho & Vishniac 2000; Maron & Goldreich 2001; Perez & Boldyrev 2008; TenBarge & Howes 2012; Mallet et al. 2015; Mallet & Schekochihin 2017), in situ solar wind observations (Horbury et al. 2008; Chen et al. 2010; Wicks et al. 2010; Chen et al. 2012; Chen 2016) and laboratory experiments (Ghim et al. 2013)⁴.

Although the current discussion is focused on RMHD, which requires a strong equilibrium magnetic field, the preceding discussion concerning critical balance can equally well apply to small scales in a system without an equilibrium magnetic field. Alfvénic fluctuations at a scale $l$ propagate along the total, local magnetic field at position $\boldsymbol {x}_0$, $\boldsymbol {B}_0^{local} (\boldsymbol {x}_0) = \boldsymbol {B}_0 + \sum _{l' \gtrsim l} \delta \boldsymbol {B}_{l'}(\boldsymbol {x}_0)$. Since the propagation and nonlinear times both decrease with scale, the fluctuations with $l' \gtrsim l$ are approximately static relative to small-scale turbulence fluctuations, and those with $l' \lesssim l$ are more rapid and will average to zero. Therefore, small-scale Alfvénic turbulence always sees an approximately constant, local mean magnetic field.

2.3. Compressible MHD

Given the complexity of incompressible turbulence, it is somewhat unsurprising that compressible MHD turbulence has received less attention. The first important point about compressible turbulence is the nature of the fast mode cascade. The fast mode propagates isotropically ($\omega _F \propto k$), and therefore the fast mode turbulence cascade is also isotropic, as confirmed in numerical simulations of compressible MHD (Cho & Lazarian 2002,  2003). In the weak limit, Chandran (2005), Luo & Melrose (2006) and Chandran (2008) have examined the wave kinetic equation in detail to explore the couplings between the Alfvén and fast modes. For quasi-parallel fluctuations ($k_\parallel > k_\perp$), the frequency and wavevector matching conditions permit any of the following three-wave interactions on approximately equal footing: $\mathrm {A}+\mathrm {A}\rightarrow \mathrm {A}$, $\mathrm {A}+\mathrm {A}\rightarrow \mathrm {F}$, $\mathrm {A}+\mathrm {F}\rightarrow \mathrm {F}$ and $\mathrm {F}+\mathrm {F}\rightarrow \mathrm {F}$. However, in the obliquely propagating limit ($k_\parallel \lesssim k_\perp$), the cascades decouple, leaving only $\mathrm {A}+\mathrm {A}\rightarrow A$ and $\mathrm {F}+\mathrm {F}\rightarrow \mathrm {F}$, because in this limit, the frequency of the fast modes exceeds significantly that of the Alfvén modes, making the frequency matching conditions involving mixed wave modes impossible. Although the wavevector and frequency matching conditions are broadened as the turbulence becomes stronger, one still expects the turbulence to be dominated by interactions that are nearby in scale (Kolmogorov 1941; Frisch 1995; Howes, TenBarge & Dorland 2011), and thus also nearby in wavevector and frequency space. Considering these facts combined with the isotropic cascade of fast modes and the anisotropic cascade of Alfvénic modes, it is expected that the cascades decouple at small scales, regardless of the initial wavevector distribution. In other words, regardless of the way a system is driven or initialized, the Alfvénic portion of the cascade will eventually decouple and behave the same as the RMHD cascade described in the preceding section. Further, moderate-amplitude fast modes rapidly form shocks as they propagate, and in weakly collisional plasmas fast modes are moderately to strongly damped via resonant wave–particle interactions (Landau 1946; Barnes 1966) for a wide range of plasma parameters (Klein et al. 2012). Therefore, in the Newtonian limit, focusing only on the Alfvénic turbulence cascade is, generally, well justified.

3. The equations of relativistic weak turbulence

3.1. Relativistic Elsasser equations

To begin our exploration of weak turbulence in relativistic, magnetically dominated plasmas, we start from an equation set that: (i) makes direct contact with the Newtonian Elsasser equations, (2.4), and (ii) highlights the fundamental role of Alfvén wave collisions in establishing the turbulence cascade. To this end, we turn to Chandran, Foucart & Tchekhovskoy (2018), who present a set of Elsasser-like equations for general relativistic MHD, including an inhomogeneous background. Here, we outline the derivation and consider only the special relativistic form with a homogeneous background. From this point forward, we employ Lorentz–Heaviside units and all speeds are normalized to the speed of light, $c=1$, and we assume a fixed Minkowski metric with signature $\eta ^{\mu \nu } = \text {diag}\{-1,1,1,1\}$. The ideal relativistic MHD equations are mass conservation

(3.1)\begin{equation} \partial_\nu(\rho u^\nu) = 0, \end{equation}

the stress energy equation

(3.2)\begin{equation} \partial_\nu T^{\mu \nu} = 0 \end{equation}

and the induction equation

(3.3)\begin{equation} \partial_\nu (b^\mu u^\nu - b^\nu u^\mu) = 0, \end{equation}

where $\rho$ is the rest-mass density, $u^\mu = (\gamma,\gamma \boldsymbol {v})$ is the four-velocity and $\gamma = 1 / \sqrt {1 - v^2}$ is the Lorentz factor. Furthermore,

(3.4)\begin{equation} T^{\mu \nu} = \mathcal{E} u^\mu u^\nu + \left(p + \frac{b^2}{2}\right) \eta^{\mu \nu} - b^\mu b^\nu \end{equation}

is the stress-energy tensor,

(3.5)\begin{equation} b^\mu = \frac{1}{2}\varepsilon^{\mu \nu \kappa \lambda} u_\nu F_{\kappa \lambda} = \left(\gamma (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{B}), \frac{1}{\gamma} \left[B^i + \gamma^2 (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{B}) v^i \right]\right) \end{equation}

is the magnetic field four-vector, $b^2 = b^\mu b_\mu$, $F_{\kappa \lambda }$ is the Faraday tensor, $\mathcal {E} = h + b^2$, $h = \rho (1 + \epsilon ) + p$ is the enthalpy density, $\varepsilon$ is the specific internal energy and $p$ is the gas pressure. Repeated indices indicate summation, Greek indices span the four-dimensional spacetime (0 to 3), while Latin indices (1 to 3) correspond to the spatial directions in a suitably chosen $3+1$ foliation of spacetime. We adopt an ideal equation of state $p = \rho \epsilon (\varGamma -1)$, with an adiabatic index $\varGamma = 4/3$ appropriate for relativistic plasmas. For this equation of state and adiabatic index, $h = \rho + 4 p$.

An Elsasser-like formulation of the relativistic MHD equations can be achieved by simply multiplying (3.3) by $\pm \sqrt {\mathcal {E}}$, adding the result to (3.2) and dividing the sum by $\mathcal {E}$, yielding

(3.6)\begin{equation} \partial_\nu (z^\mu_\pm z^\nu_\mp{+} \varPi \eta^{\mu \nu}) + \left(\frac{3}{4} z^\mu_\pm z^\nu_\mp{+} \frac{1}{4} z^\mu_\mp z^\nu_\pm{+} \varPi \eta^{\mu \nu}\right) \frac{\partial_\nu \mathcal{E}}{\mathcal{E}} = 0, \end{equation}

where

(3.7)\begin{equation} z^\mu_\pm{=} u^\mu \pm \frac{b^\mu}{\sqrt{\mathcal{E}}} \end{equation}

are the relativistic Elsasser four-vector fields and

(3.8)\begin{equation} \varPi = \frac{2p + b^2}{2\mathcal{E}} = \frac{1}{2} - \frac{2 p + \rho}{2\mathcal{E}}. \end{equation}

Much like in non-relativistic MHD, (3.6) represents the evolution of upward- and downward-propagating Elsasser fields; however, (3.6) represents the fully compressible system and is thus not completely analogous to (2.4). In the relativistic limit, one cannot simply go to the incompressible limit typically used to derive incompressible MHD, because the maximum speed is the speed of light, and thus compressible fluctuations cannot be ordered in such a way as to instantaneously carry information away from a source. We can, however, consider a limit similar to that of RMHD to isolate the Alfvénic fluctuations described by (3.6). Such a limit, in which fluctuations are highly elongated relative to a mean magnetic field, is reasonable to consider for many astrophysical plasmas which often have strong mean fields, e.g. black hole accretion disks and jets (Narayan, Igumenshchev & Abramowicz 2003), coronae (Chandran et al. 2018; Yuan, Blandford & Wilkins 2019) and magnetar magnetospheres (Parfrey, Beloborodov & Hui 2012). Further, as we show in § 4, relativistic Alfvénic turbulence proceeds in the same anisotropic sense as Newtonian turbulence, eventually leading to highly anisotropic, small-amplitude fluctuations at small scales, regardless of the outer-scale conditions. Additionally, many astrophysical systems are characterized by exceptionally large inertial ranges⁵, vastly exceeding the three- to four-order-of-magnitude-wide inertial range observed in the solar wind near Earth (Howes et al. 2008; Bruno & Carbone 2013; Kiyani, Osman & Chapman 2015; Chen 2016).

3.2. Relativistic RMHD

To derive a set of Elsasser equations for RRMHD, we begin by separating the fluid into mean (background) quantities and fluctuating quantities of the form

(3.9)\begin{equation} b^\mu = \langle b^\mu \rangle + \delta b^\mu. \end{equation}

We also construct an average fluid rest frame in which $\langle u^i\rangle$ vanishes⁶. We define $\lambda$ to be the correlation length transverse to $\langle B^i \rangle$, and $L$ to be the characteristic scale of the background, i.e. the length scale parallel to $\langle B^i \rangle$. As in RMHD, we assume $\lambda / L \sim \mathcal {O}\left (\epsilon \right )$, where $\epsilon \ll 1$. In addition to the anisotropy assumption, we also assume that the characteristic frequency is of the order of the Alfvén frequency, $\partial _t \sim \langle B^i\rangle \partial _i \sim k_\parallel v_A$, and that the fluctuations are ordered small:

(3.10)\begin{gather} \frac{\sqrt{\delta u^2}}{v_A} \equiv \frac{\sqrt{\delta u^\mu \delta u_\mu}}{v_A} \sim \sqrt{\delta b^2} \sim \mathcal{O}\left(\epsilon\right), \end{gather}

(3.11)\begin{gather}\frac{\delta \rho}{\langle\rho\rangle} \sim \frac{\delta p}{\langle p\rangle} \sim \mathcal{O}\left(\epsilon^2\right), \end{gather}

where

(3.12)\begin{equation} v_A^i = \frac{\langle b^i \rangle}{\sqrt{\langle\mathcal{E}\rangle}} \end{equation}

and $v_A = \sqrt {v_A^i v_{Ai}}$. We further assume that all background quantities, e.g. $\langle B^i\rangle = \boldsymbol {B}_0, \langle \rho \rangle \equiv \rho _0,\langle p\rangle \equiv p_0$, are constant with no background inhomogeneities.

As in Newtonian RMHD, in RRMHD, the fast wave will be ordered out of the system, since $\omega _F \sim k \gg \omega \sim k_\parallel$. Relativistic RMHD will describe Alfvén and pseudo-Alfvén fluctuations; however, we are not interested in the pseudo-Alfvén waves, which have fluctuations parallel to the background magnetic field and are a separate, passive, cascade channel. Therefore, we also assume $B_{0i} \delta z^i_\pm = \mathcal {O}\left (\epsilon ^2\right )$ to remove the pseudo-Alfvén modes. With all of the above assumptions, we note that $\gamma \sim 1 + \mathcal {O}\left (\epsilon ^2\right )$, $\partial _i \delta u^i \sim \partial _\mu \delta z^\mu _\pm \sim \langle b^t\rangle \sim \delta b^t \sim \delta u^t / v_A \sim \mathcal {O}\left (\epsilon ^2\right )$.

Applying the above reduced assumptions to the relativistic Elsasser equation (3.6), we arrive at

(3.13)\begin{equation} v^\nu_{A\mp} \partial_\nu \delta z^\mu_\pm{=}{-}\delta z^\nu_\mp \partial_\nu \delta z^\mu_\pm{-} \partial_\nu (\delta \varPi \eta^{\mu \nu}), \end{equation}

where to lowest non-zero order

(3.14)\begin{gather} v_{A\pm}^\mu = \left(1,\pm\frac{\boldsymbol{B}_0}{\sqrt{\mathcal{E}_0}}\right) \equiv (1,\pm \boldsymbol{v}_A), \end{gather}

(3.15)\begin{gather}\delta z^\mu_\pm{=} \delta u^\mu \pm \frac{\delta b^\mu}{\sqrt{\mathcal{E}_0}} = \left(0,\delta \boldsymbol{v}_\perp{\pm} \frac{\delta \boldsymbol{B}_\perp}{\sqrt{\mathcal{E}_0}}\right), \end{gather}

(3.16)\begin{gather}\delta \varPi ={-}\frac{2 \delta p + \delta \rho}{2 \mathcal{E}_0} + \frac{2 p_0 + \rho_0}{2 \mathcal{E}_0}\frac{\delta \mathcal{E}}{\mathcal{E}_0}, \end{gather}

$\mathcal {E}_0 = 4 p_0 + \rho _0 + B_0^2$ and $\delta \mathcal {E} = 4\delta p + \delta \rho + 2 B_0 \delta B_\parallel$. In three-vector form, the equation set is particularly simple and similar to the Newtonian Elsasser equations:

(3.17)\begin{equation} \frac{\partial \delta \boldsymbol{z}_{{\perp}{\pm}}}{\partial t} \mp \boldsymbol{v}_A \boldsymbol{\cdot} \boldsymbol{\nabla} \delta \boldsymbol{z}_{{\perp}{\pm}} ={-} \delta \boldsymbol{z}_{{\perp}{\mp}} \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \delta \boldsymbol{z}_{{\perp}{\pm}} - \boldsymbol{\nabla}_\perp \delta \varPi. \end{equation}

The terms of the left-hand side of (3.13) are linear, while those on the right-hand side are nonlinear. Due to the assumptions regarding pressure and parallel fluctuations, the only relevant components of (3.13) are the two components transverse to the mean magnetic field, as expressed in (3.17). The other two components, time-like and parallel, are one order higher in the ordering parameter, $\epsilon$. Importantly, the parallel fluctuations do not appear at lowest order in the nonlinear $\boldsymbol {E} \times \boldsymbol {B}$ term in the perpendicular, Alfvén wave equations, (3.13) and (3.17), and they therefore do not cascade the Alfvén waves. However, the primary nonlinear term for the parallel fluctuations is of the form $\delta \boldsymbol {z}_{\perp _\pm } \boldsymbol {\cdot } \boldsymbol {\nabla }_\perp \delta z_\parallel$. Thus, the parallel fluctuations are passively scattered/mixed by the Alfvén waves, much like in Newtonian RMHD turbulence (Schekochihin et al. 2009). Note that as in the Newtonian RMHD equations, the system is closed by taking the divergence of (3.17), or four-divergence of (3.13), to find an equation for $\delta \varPi$.

The reduced relativistic Elsasser equations, (3.17), for RRMHD are identical in form to the Newtonian RMHD equations, and at this point the standard approach to solve for the Alfvén dynamics would be to take the curl of (3.17) to eliminate the pressure term and solve for the Elsasser potentials rather than the Elsasser fields. Indeed, this is the approach we will also take; however, we note that it is possible to simplify the system even further, because there is a straightforward limit to consider in relativistic plasmas. This final simplification to consider is to remove the pressure fluctuations by assuming that we are in the magnetically dominated regime, $\sigma = b^2 / h \gg 1$. In this limit, $\langle \varPi \rangle = 1/2$, and $\delta \varPi \sim \mathcal {O}\left (\epsilon ^2\right ) / \sigma \ll \mathcal {O}\left (\epsilon ^2\right )$. Therefore, we arrive at the magnetically dominated RRMHD equations:

(3.18)\begin{equation} v^\nu_{A\mp} \partial_\nu \delta z^\mu_\pm{=}{-}\delta z^\nu_\mp \partial_\nu \delta z^\mu_\pm. \end{equation}

In three-vector form, the equation set is particularly simple:

(3.19)\begin{equation} \frac{\partial \delta \boldsymbol{z}_{{\perp}{\pm}}}{\partial t} \mp \boldsymbol{\nabla}_\parallel \delta \boldsymbol{z}_{{\perp}{\pm}} ={-} \delta \boldsymbol{z}_{{\perp}{\mp}} \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp \delta \boldsymbol{z}_{{\perp}{\pm}}. \end{equation}

Being in the $\sigma \gg 1$ limit, (3.18) and (3.19) are essentially the reduced analogue of the force-free limit of the ideal MHD equations (Gruzinov 1999; Komissarov 2002). In this limit, $v_A = c = 1$.

To arrive at this simple form, it may seem that we have rather seriously brutalized the relativistic MHD equations by applying a series of restrictive asymptotic orderings: (i) $k_\parallel \ll k_\perp$; (ii) $\omega \sim k_\parallel v_A$; (iii) small-amplitude fluctuations with constant (mean) backgrounds; (iv) second-order parallel and pressure fluctuations; and (v) magnetically dominated, $\sigma \gg 1$. The reduced assumptions (i)–(iii) and (v) are consistent with one another, i.e. they can all be obtained by assuming the system is strongly magnetized. In the Newtonian limit, the fast mode can be eliminated by assuming the system is incompressible; however, the finite speed of light prevents easy elimination of the fast mode (Takamoto & Lazarian 2017) without assuming the fluctuations are anisotropic ($k_\parallel \lesssim k_\perp$). To remove the slow mode, we assumed that parallel and pressure fluctuations are second-order quantities when deriving the reduced relativistic Elsasser equations. Moving to the magnetically dominated regime produces the same result, since in the $\sigma \rightarrow \infty$ limit, the slow wave is also ordered out of the system. Thus, despite the myriad assumptions to achieve a simple set of equations for describing relativistic Alfvénic turbulence, they are all self-consistent. This set of assumptions is also relevant to many astrophysical systems, such as magnetars (Li & Beloborodov 2015; Li et al. 2019; Yuan et al. 2020b), glitches affecting radio emission from pulsar magnetospheres (Bransgrove, Beloborodov & Levin 2020; Yuan et al. 2020a) and X-ray-emitting coronae around black hole accretion disks (Thompson & Blaes 1998; Chandran et al. 2018). We also note that the most important assumption is wavevector anisotropy, $k_\parallel \ll k_\perp$, which naturally arises as the Alfvénic turbulence cascades to small scales.

3.3. Connection to and comparison with Newtonian RMHD solutions

Equation (3.17) will form the basis for the following analysis in § 4, because it represents the simplest set of equations for describing Alfvénic turbulence in magnetized, relativistic environments and has a form that is nearly identical to the Newtonian RMHD Elsasser equations. Thus, the system shares many properties with the Newtonian RMHD system of equations, some of which can be seen immediately: (i) the system supports linear Alfvén modes, which require $k_\parallel \ne 0$ to propagate; (ii) the nonlinearity requires counter-propagating fluctuations; and (iii) the nonlinearity requires that the fluctuations be polarized with respect to each other in the perpendicular plane. Further, the wave kinetic equations for the system coincide with the kinetic equations for shear Alfvén waves derived in Galtier et al. (2000) in the incompressible MHD limit and Boldyrev & Perez (2009) in the RMHD limit neglecting imbalanced interactions, i.e. non-zero cross helicity. As demonstrated in those papers, the wave kinetic equation evolves the spectral energy $\textrm {e}^{\pm } = \langle |\delta \boldsymbol {z}_\perp ^\pm (\boldsymbol {k})|^2 \rangle$ and can be expressed as

(3.20)\begin{equation} \frac{\partial \,\textrm{e}^\pm(\boldsymbol{k})}{\partial t} = \int M_{k,pq} \textrm{e}^\mp(\boldsymbol{q}) [\textrm{e}^\pm(\boldsymbol{p}) - \textrm{e}^\pm(\boldsymbol{k})] \delta(q_\parallel)d_{k,pq}, \end{equation}

where $M_{k,pq} = ({\rm \pi} / v_A) (\boldsymbol {k}_\perp \times \boldsymbol {q}_\perp )^2 (\boldsymbol {k}_\perp \boldsymbol {\cdot } \boldsymbol {p}_\perp )^2 / (k_\perp q_\perp p_\perp )^2$ and $d_{k,pq} = \delta (\boldsymbol {k} - \boldsymbol {p} - \boldsymbol {q})d^3p \ d^3q$. Since the kinetic equation is unchanged from RMHD, we also know that the weak turbulence solutions: (i) are dominated by three-wave interactions and (ii) lead to an energy spectrum of the form $f(k_\parallel ) k_\perp ^{-2}$, where $f(k_\parallel )$ is set by external forcing or initial conditions. In the following sections, we explore in detail the three-wave interaction of Alfvén waves first analytically via heuristic solutions through third order, and then numerically to demonstrate the decoupling of the fast mode from the Alfvén waves in the obliquely propagating (reduced) limit. Note that although we employ the RRMHD equations (3.17) to derive the turbulence solutions in the following section, the Alfvén solutions are identical for the $\sigma \rightarrow \infty$ limit of the RRMHD equations (3.19).

4. Three-wave weak Alfvénic interactions

To construct analytical solutions, we now consider a subsidiary expansion of the form $\zeta _\pm = \varepsilon \zeta _{1\pm } + \varepsilon ^2 \zeta _{2\pm } + \cdots$, where $\delta \boldsymbol {z}_\pm = \hat {\boldsymbol {z}} \times \boldsymbol {\nabla }_\perp \zeta _\pm$ defines the Elasser potential, $\zeta$, and $\zeta _{i \pm }(t=0) = 0$ for $i > 1$. Note that since the relativistic equations are identical in form to the Newtonian limit, the solution to the RRMHD equations is also the same, and the full solutions appear in Appendix A. Therefore, we primarily summarize the solutions found in Howes & Nielson (2013); however, we note that we include corrections to equations (22), (28)–(29) of Howes & Nielson (2013), and all equations involving $\zeta _{3-}$, including the equations for $\boldsymbol {B}_{\perp 3}$ and $\boldsymbol {E}_{\perp 3}$. These errors have been corrected in Appendix A.

For specificity, we will consider a periodic domain of size $L_x \times L_y \times L_z$ with $\boldsymbol {B}_0 = B_0 \hat {\boldsymbol {z}}$. At $t=0$, we will initiate two counter-propagating, perpendicularly polarized fluctuations:

(4.1)\begin{equation} \left.\begin{gathered} \boldsymbol{z}_{1+} = z_+ \cos{(k_\perp x - k_\parallel z - \omega_0 t)}\hat{\boldsymbol{y}},\\ \boldsymbol{z}_{1-} = z_- \cos{(k_\perp y + k_\parallel z - \omega_0 t)}\hat{\boldsymbol{x}}, \end{gathered}\right\} \end{equation}

where $z_\pm$ are the initial amplitudes, $k_\perp = 2{\rm \pi} / L_x = 2 {\rm \pi}/ L_y$ and $k_\parallel = 2{\rm \pi} / L_z$. We maintain the convention that $k_\perp$, $k_\parallel$ and $\omega _0$ are positive constants, and the direction of propagation is supplied by the explicit sign of $k_\parallel$. To describe the wave modes that arise from the nonlinear evolution, we use the notation $(k_x/k_\perp, k_y/k_\perp, k_z / k_\parallel )$. For instance, the plus and minus initial wave modes in this notation are $(1,0,-1)$ and $(0,1,1)$, respectively.

The solutions through third order are lengthy, and as such can be found in Appendix A. Here, we summarize the important findings at each order. The initial conditions provided in (4.1) satisfy the lowest-order equations if $\omega _0$ is the linear Alfvén wave frequency, $\omega _0 = k_\parallel v_A$. Thus, at lowest order the solution describes counter-propagating, linear Alfvén waves.

The second-order solutions for the electric and magnetic field are given by (A 14) and (A 15), from which we can immediately see that at second order, all components are purely oscillatory, i.e. there is no secularly growing mode. We can also see that two Fourier modes are generated at this order, $(1,1,0)$ and $(-1,1,2)$. Both of these modes satisfy the wavevector matching conditions, which require $\boldsymbol {k}_2 = \boldsymbol {k}_{1-} \pm \boldsymbol {k}_{1+} = (0,1,1) \pm (1,0,-1)$. The $(-1,1,2)$ mode satisfies the conditions for a linear Alfvén wave. Specifically, the frequency of this mode is $\omega = \pm 2\omega _0 = \pm 2 k_\parallel v_A = \pm k_z v_A$, and the mode obeys the Alfvén wave eigenrelation

(4.2)\begin{equation} \frac{\delta \boldsymbol{B}_\perp}{\sqrt{\mathcal{E}_0}} ={\pm} \delta \boldsymbol{u}_\perp{=} \pm \frac{\boldsymbol{E}_\perp}{B_0} \times \hat{\boldsymbol{z}}. \end{equation}

These two linear modes are counter-propagating along the background magnetic field; however, they are not perpendicularly polarized. Therefore, their interaction is simply a linear superposition that forms a standing wave for this particular symmetric initial condition. The $(1,1,0)$ mode is a purely magnetic mode that has no structure along the background field ($k_z = 0$) but oscillates with frequency $\omega = 2 \omega _0$. This mode is not a linear Alfvén mode and corresponds to a nonlinear magnetic shear. The interaction of this shear mode with the initial wave modes is the nonlinear interaction that will provide secular growth at the next order.

As with the second-order equations for $\boldsymbol {B}_\perp$ and $\boldsymbol {E}_\perp$, we can once again straightforwardly interpret the third-order solutions, (A 23) and (A 24). First, we note that there are now secularly (proportional to $t$) growing modes which are boxed for clarity: $(2,1,-1)$ and $(1,2,1)$. Both of these modes are linear Alfvén waves. They have frequency $\omega = \pm \omega _0 = \pm k_\parallel z = \pm k_z z$ and obey the Alfvén wave eigenrelations (4.2). Therefore, the modes correspond to linear Alfvén waves propagating up, $(1,2,1)$, and down, $(2,1,-1)$, the background magnetic field. These secularly growing modes are phase-shifted relative to the initial modes by $-{\rm \pi} / 2$, have perpendicular wavevectors $\sqrt {5} k_\perp$, but $|k_z| = k_\parallel$. Therefore, the weak turbulence cascade proceeds to smaller scales in the perpendicular plane, but the parallel scales are conserved, confirming the prediction from simple three-wave matching conditions. Further, as alluded to in the previous section, these secular modes follow from the interaction of the second-order $(1,1,0)$ mode with the initial Alfvén waves: $\boldsymbol {k}_3 = \boldsymbol {k}_{2} + \boldsymbol {k}_{1\pm } = (1,1,0) + (1,0,-1) = (2,1,-1)$ and $(1,1,0) + (0,1,1) = (1,2,1)$, i.e. each third-order mode conserves not just the magnitude but also the sign of $k_z$. This fact confirms that there is no exchange of energy between upward- and downward-propagating fluctuations.

The other six purely oscillatory components are a mixture of linear Alfvén waves and nonlinear structures. Unlike the second-order solutions, there is no purely magnetic mode at third order. Four of the components, those with wavevectors $(2,1,-1), (1,2,1), (-2,1,3)$ and $(-1,2,3)$, have $\sqrt {5} k_\perp$, but these all have both linear and nonlinear components. Similarly, the $(0,1,1)$ and $(1,0,-1)$ components appear as both linear and nonlinear terms. Note that these final two components have the same wavevector as the initial Alfvén waves, but these third-order modes have different phases and will serve to cancel the initial modes.

5. Numerical comparison

5.1. Simulation description

To confirm the analytical results of §§  3 and 4, we consider the nonlinear interaction between two perpendicularly polarized Alfvén waves that counter-propagate in a three-dimensional, periodic domain along a uniform guide field $\boldsymbol {B_0} = B_0 \boldsymbol {\hat {z}}$ using the general relativistic MHD code BHAC (Porth et al. 2017; Olivares et al. 2019; Ripperda et al. 2019a, ). The recently added force-free limit for the resistive MHD code BHAC employs the numerical scheme of Ripperda et al. (2019a, ) and damps force-free violations $E>B$ and $\boldsymbol {E} \boldsymbol {\cdot } \boldsymbol {B} \neq 0$ on resistive time scales⁷. We employ both a periodic cubic domain with $L_{\perp }=L_x=L_y=L_{\parallel } = L_z = 2{\rm \pi}$ and a periodic elongated domain with $L_{\parallel } = 10 L_\perp = 20{\rm \pi}$, with resolution $N_x = N_y = N_z = 256$ and $N_z = 2560$ cells for the elongated case. Initially, we set a gas-to-magnetic-pressure ratio of $\beta =2p/B_0^2=0.02$ and magnetization $\sigma _{cold} \equiv b^2 / \rho \in [0.01;0.1;1;10;100]$ (corresponding to $\sigma \in [0.01;0.1;1;7;20]$), where we vary the density, $\rho$, but maintain a constant guide field, $\boldsymbol {B_0} = \boldsymbol {\hat {z}}$, and pressure $p=0.01$ with adiabatic index $\varGamma =4/3$ for an ideal relativistic gas. In the force-free case, $\beta \rightarrow 0$, $\sigma \rightarrow \infty$ and $v_A \rightarrow 1$.

As in § 4, we prescribe a scale-free definition of characteristic wavelengths $(k_x/k_\perp, k_y/k_\perp,$ $k_z / k_\parallel )$, and initialize our Alfvén wave simulations with two overlapping, counter-propagating and perpendicularly polarized Alfvén waves described by the wavevectors $\boldsymbol {k}_{+}=(1,0,-1)$ and $\boldsymbol {k}_{-}= (0,1,1)$. The magnetic field is initialized through a vector potential $\boldsymbol {A} = (-B_0 y, 0, \delta B_\perp [\sin (k_\perp x+ k_\parallel z) + \sin (k_\perp y-k_\parallel z)])$, representing the initial counter-propagating Alfvén waves with frequency $\omega _0 = k_{\parallel } v_A$. The electric field is initialized as $\boldsymbol {E} = (v_A B_y, v_A B_x,0)$ such that the velocity is equal to the drift velocity, $\delta \boldsymbol {u}_\perp = \boldsymbol {E} \times \boldsymbol {B}/B^2$. Note that the overlapping Alfvén waves, in contrast to a single Alfvén wave, are not an exact force-free equilibrium due to a small second-order violation $\boldsymbol {E}_{\mp } \boldsymbol {\cdot } \boldsymbol {B}_{\pm } \neq 0$ between the fields of the two waves, which is damped on a short time scale in BHAC.

The strength of the nonlinearity is characterized by $\chi = k_{\perp } \delta B_{\perp } / k_{\parallel } B_0$. To maintain weak turbulence, we fix $\chi = 0.01$ in all of the following simulations. By fixing $\chi$, we expect to maintain self-similar behaviour as we explore elongated domains that approximate the reduced limit. Thus, in the cubic domain with $k_{\perp } = k_{\parallel }$, $\delta u_{\perp } / v_A = \delta B_\perp / B_0 = 0.01$, and in the elongated case with $k_{\perp } = 10 k_{\parallel }$, $\delta u_{\perp } / v_A = \delta B_\perp / B_0 = 0.001$.

5.2. Alfvén wave–Alfvén wave collisions

5.2.1. Cubic domain

In figure 1, we present the evolution of the mode amplitudes of the $B_\perp$ (Alfvénic) fluctuations in a cubic domain scanning $\sigma _{cold} \equiv b^2 / \rho \in [0.01;0.1;1;10;100;\infty ]$, presented in descending order of $\sigma$ from top-left to bottom-right. The initial, counter-propagating and perpendicularly polarized Alfvén waves are represented by the red lines in each panel, the $(0,1,-1)$ mode by dashed lines and the $(1,0,1)$ mode by dotted lines. These modes interact to produce a secondary, nonlinear, magnetic shear mode indicated by the green line, representing the $(1,1,0)$ mode. Note that this mode is purely oscillatory in time with $\omega = 2 \omega _0$, as found in § A.3. Finally, the secondary mode interacts with each of the primary Alfvén waves to produce at third order two higher $k_\perp$ Alfvén waves, but with the same $k_\parallel$ as the primary waves, where the $(1,2,-1)$ and $(2,1,1)$ modes are represented by the blue dashed and dotted curves. These dynamics are in agreement with the results of § A.4. These modes grow secularly in time, $B_{\perp 3} \propto t$, as indicated by the black line in each panel, which is a curve proportional to time and scaled by the final amplitude of $B_{\perp 3}$. The net result of this Alfvén wave–Alfvén wave collision is the anisotropic transfer of energy from large to small scales, which is governed by three-wave interactions.

Figure 1. Mode evolution of $B_{\perp }$ for a cubic domain with $\chi =0.01$ and $\sigma _{cold} \equiv b^2 / \rho \in [0.01;0.1;1;10;100;\infty ]$, presented in descending order of $\sigma$ from (a) to (f). The red lines correspond to the initial counter-propagating Alfvén modes: the $(0,1,-1)$ and $(1,0,1)$ modes are shown by dashed and dotted lines, respectively. Dashed green lines correspond to the second-order nonlinear shear mode $(1,1,0)$. The $(1,2,-1)$ and $(2,1,1)$ third-order, secularly growing modes are shown by dashed and dotted blue lines, respectively. The black line corresponds to a secular growth directly proportional to $t$ and scaled by the final amplitude of the third-order mode.

5.2.2. What about fast waves? Exploring the Alfvén wave–fast wave coupling

Our analytical analysis in the preceding sections purposefully neglected the fast wave by choosing the reduced limit. However, the fast wave may play an important role in releasing energy in strongly magnetized, relativistic, astrophysical plasmas, because the fast wave can travel across field lines (Li et al. 2019; Yuan et al. 2020b). For instance, Alfvén waves travel along magnetic fields, and on closed magnetic field lines in pulsar and magnetar magnetospheres their energy is confined to the magnetosphere, but fast waves can release their energy into the surrounding medium by propagating across field lines. For this reason, the coupling between the fast and Alfvén branches in relativistic MHD and its force-free limit has been explored in detail in recent years (Cho 2005; Takamoto & Lazarian 2016,  2017; Li et al. 2019). Through a numerical simulation study of relativistic turbulence, Takamoto & Lazarian (2016,  2017) find that the fast-to-Alfvén mode power scales as $(\delta u_f / \delta u_A)^2 \propto \sqrt {(1+\sigma )} \delta u_A / c_{f\perp }$ in the $\sigma \gtrsim 1$ limit, where $c_{f\perp }$ is the fast mode speed in the perpendicular plane. In the $\sigma \ll 1$ limit, Cho & Lazarian (2002,  2003) find that $(\delta u_F / \delta u_A)^2 \propto \delta u_A / c_{f\perp }$. Thus, for $\sigma \ll 1$, the fast and Alfvén branches decouple as the background magnetic field strength is increased; however, for $\sigma \gtrsim 1$, the coupling between the modes increases with $\sqrt {1 + \sigma }$, asymptotically approaching unity for large $\sigma$.

In figure 2, we present the evolution of the highest-amplitude $B_\parallel$ modes for the same initial configuration as presented in figure 1. Here $B_\parallel$ serves as a proxy for the amplitude of compressible, fast mode fluctuations, since the Alfvén modes have negligible $B_\parallel$ components. As an additional test (not shown), we have projected the fluctuations onto the fast and Alfvén wave eigenfunctions and recovered quantitatively similar results.

Figure 2. Mode evolution of $B_{\parallel }$ for a cubic domain with $\chi =0.01$ and $\sigma _{cold} \equiv b^2 / \rho \in [0.01;0.1;1;10;100;\infty ]$, presented in descending order of $\sigma$ from (a) to (f). Line styles and colours are as defined in figure 1.

We first note that the $(0,1,-1)$ and $(1,0,1)$ $B_\parallel$ modes have zero amplitude at $t=0$ but do develop finite amplitude that is approximately independent of $\sigma$. These modes appear at a mixture of two or more frequencies, the most dominant of which is at the fast mode frequency, $\omega \sim \omega _f \sim k c_f \sim \sqrt {2} \omega _0$, and remain at the same or lower amplitude as the tertiary (blue) modes. The secondary mode, $(1,1,0)$ (green), is the highest-amplitude mode and is again composed of two or more frequencies, the most dominant of which are $\omega \sim \omega _f \sim k c_f \sim \sqrt {2} \omega _0$ and $\omega \sim \omega _0 / 2$. Unlike the Alfvén branch, the $(1,1,0)$ mode can be a linear fast mode, which is consistent with a component frequency being at the fast mode frequency. Finally, the tertiary modes (blue) again display secular growth with time, with a primary frequency given by the fast mode frequency $\omega \sim \omega _f \sim k c_f \sim \sqrt {6} \omega _0$. The tertiary mode also displays the most significant dependence on $\sigma$. For $\sigma \ll 1$, the amplitude is independent of $\sigma$.

5.2.3. Elongated domain

The scaling of the Alfvén–fast mode coupling discussed above focuses purely on the strength of the background magnetic field while neglecting any wavevector anisotropy that may naturally arise from a strong magnetic field. As noted in § 2.3, in the Newtonian, obliquely propagating limit ($k_\parallel \lesssim k_\perp$), the fast and Alfvén branch turbulent cascades decouple, and we expect an analogous decoupling to occur in the relativistic limit, since the fast and Alfvén wave frequencies remain well separated in the oblique limit. Thus, we now examine the results of an elongated box simulation with $L_\parallel = 10L_\perp$.

In figure 3(a,b) are presented the mode analysis results of a force-free simulation of obliquely propagating Alfvén waves, with $k_\perp = 10 k_\parallel$. The colour and line styles are as in the previous figures. Comparing the relative amplitudes between each of the $B_\perp$ modes with the relative amplitudes for the cubic domain case presented in figure 1(a), it is clear that the elongated domain does not change the Alfvénic $B_\perp$, cascade. However, the three parallel modes’ growth, amplitude and frequency in the elongated case compared with the cubic case in figure 2(a) are markedly different. First, no secular growth of parallel modes is apparent. If there is a secularly growing mode, it is sufficiently low in amplitude that it does not modulate or interfere with the non-growing, purely oscillatory modes. Second, the mode amplitudes relative to the Alfvén wave primaries are an order of magnitude smaller in the elongated case relative to the cubic domain, which is consistent with the compressible mode amplitudes scaling with the elongation factor. Third, in the elongated case, the first-order (red) and third-order (blue) parallel modes are dominated by low-frequencyoscillations that match the frequency of the Alfvénic secondary mode. The second-order (green) mode is much higher in frequency, well above even the fast mode frequency, suggesting nonlinear modes with multiple frequencies contribute to this component.

Figure 3. Mode evolution of $B_{\perp }$ and $B_{\parallel }$ with $\chi =0.01$. (a,b) Results for an Alfvén wave–Alfvén wave collision in an elongated domain, $L_z = 10 L_x$. (c–f) Results for an Alfvén wave–fast wave collision in a cubic domain (c,d) and elongated domain (e,f). The line style and colour are as in previous figures.

These results support our intuition that the Alfvén and fast mode cascades decouple as the wavevector anisotropy increases and extend the findings of Chandran (2005) to the relativistic limit. In the case examined here, $\omega _A = k_\parallel \simeq 0.1 \omega _F = 0.1 k$, leading to a poor frequency matching condition required for the resonant weak turbulence interaction. Although the wave matching condition broadens in strong turbulence, the Alfvén and fast mode cascades will remain well separated in the limit $k_\parallel \ll k_\perp$.

5.3. Alfvén wave–fast wave collisions

In the Alfvén wave–Alfvén wave collision case, we have shown that the amplitude of the product fast modes decreases as the initial modes become anisotropic, indirectly demonstrating that the modes decouple. Here, we directly diagnose the coupling of the Alfvén and fast mode cascades by examining the collision between an Alfvén wave and a fast wave in both cubic and elongated domains. Unlike the Alfvén wave–Alfvén wave collision explored in the previous sections, there is not a canonical choice for the counter-propagating fast mode. Therefore, we choose a fiducial fast mode with $\boldsymbol {k} = (0,1,-1)$, and an Alfvén wave with $\boldsymbol {k} = (1,0,1)$. This choice is largely for consistency with the Alfvén wave–Alfvén wave collision; however, it also obeys reasonable wavevector and frequency matching conditions, and it has the same expected product modes as the Alfvén wave–Alfvén wave collision. Further, the secondary modes will still be $(1,1,0)$ and $(-1,1,-2)$, and for these initial modes in a cubic domain, the secondary $(1,1,0)$ mode can satisfy the linear fast mode dispersion relation. For these simulations, we continue the convention of fixing $\chi = 0.01$ to satisfy weak turbulence.

5.3.1. Cubic domain

In figure 3(c,d) are plotted the results of the Alfvén wave–fast wave collision outlined above for a cubic domain. Note that in this case, the dotted and dashed red lines in the $B_\perp$ panel correspond to the initial Alfvén and fast waves respectively. Focusing on the behaviour in figure 3(c), the secondary, green, mode is the lowest-amplitude mode, displays no apparent secular growth and contains a mixture of linear and nonlinear frequencies. Similarly, the minus tertiary mode is oscillatory and dominated by high-frequency nonlinear modes. However, the plus tertiary mode does exhibit secular growth, and is dominated by the linear Alfvén wave frequency. Relative to the Alfvén wave–Alfvén wave collision in figure 1(a), this tertiary mode is decreased in amplitude by approximately a factor of three.

Turning to the Alfvén wave–fast wave in figure 3(d), we see different behaviour compared with figure 2(a). The largest-amplitude component remains the secondary mode, which is dominated by the fast mode linear frequency and an additional low-frequency component. Both the plus mode ‘primary’ and tertiary exhibit secular growth with comparable amplitude. Note that the plus mode $B_\parallel$ ‘primary’ is not initialized since the initial plus mode is an Alfvén wave. This ‘primary’ mode growth is due to the interaction of the $(-1,1,2)$ secondary with the fast mode primary with $(0,1,-1)$ summing to give $(1,0,1)$, i.e. the ‘primary’ is here a tertiary mode.

5.3.2. Elongated domain

In figure 3(e,f) are plotted the results of the Alfvén wave–fast wave collision outlined above for the elongated domain with $L_\parallel = 10L_\perp$. The overall behaviour is similar to that of the Alfvén wave–fast wave collision in the cubic domain with a few notable differences. Focusing first on figure 3(e), we see that the Alfvénic tertiary is further suppressed in amplitude, now reduced by nearly two orders of magnitude relative to the cubic and elongated domain Alfvén wave–Alfvén wave collisions in figures 1(a) and 3(a). The secondary and minus tertiary are also further suppressed in amplitude, and the secondary is dominated by higher-frequency modes relative to the cubic domain. Turning to the Alfvén wave–fast wave in figure 3(f), we again see that all product modes are suppressed by approximately an order of magnitude, consistent with the elongation factor weakening the interaction. Notably, both blue tertiary modes display no secular growth. The secondary mode is of higher frequency, and the ‘primary’ mode is the only $B_\parallel$ component with clear secular growth; however, the ‘primary’ is now dominated by a high-frequency, nonlinear component. These results further support our intuition from Newtonian weak turbulence (Chandran 2005) that the Alfvén and fast mode cascades decouple as the wavevector anisotropy increases by directly comparing the Alfvén wave–fast wave collision in both cubic and elongated domains extend the findings of Chandran 2005 to the relativistic limit.

6. Summary and conclusions

In this paper, we present an analytical and numerical study of relativistic, weak Alfvénic turbulence focusing on the three-wave interaction. We begin by reviewing the knowledge gained from non-relativistic (Newtonian) turbulence and use that framework to guide our study of relativistic turbulence. From Newtonian turbulence theories, we know that both weak and strong Alfvénic turbulence lead to anisotropic cascades. Therefore, regardless of the strength of the turbulence and isotropy at large scales, turbulence will become anisotropic and eventually satisfy $k_\parallel \ll k_\perp$. Leveraging the anisotropic cascade of energy and the fact that $\delta B/B_0 \ll 1$ is satisfied either at the outer scale or once the cascade reaches sufficiently small scales, we derive a set of RRMHD equations and cast them in Elsasser form. This set of equations has the same properties and basic form as RMHD and is appropriate for describing strongly magnetized relativistic plasmas for which Alfvénic fluctuations dominate. We note that RRMHD is not equivalent to the magnetically dominated equations of the force-free limit of relativistic MHD, which require $\sigma = b^2 / h \rightarrow \infty$. However, we do also derive a set of magnetically dominated RRMHD equations which are valid in the $\sigma \rightarrow \infty$ limit but retain the favourable properties inherent to RMHD for the analysis of Alfvénic turbulence.

We note that the similarity between RMHD and RRMHD extends to the wave kinetic equation, which describes the spectral evolution of a weakly turbulent plasma. As such, we conclude that RRMHD: (i) is dominated by three-wave interactions of Alfvén waves and (ii) leads to an energy spectrum of the form $f(k_\parallel ) k_\perp ^{-2}$, where $f(k_\parallel )$ is set by external forcing or initial conditions. We then employ the RRMHD equations in Elsasser form to analytically compute the building blocks of weak turbulence through third order to heuristically demonstrate the primacy of the three-wave interaction. The primary interaction is the collision of perpendicularly polarized Alfvén waves with $\omega = \omega _0$ and wavevectors $\boldsymbol {k}_{1\pm } = (1,0,1)$ and $(0,1,-1)$ that interact through a three-wave interaction to produce a nonlinear, magnetic shear mode, $\boldsymbol {k}_2 = (1,1,0)$, with $\omega = 2\omega _0$ at second order. The interaction of this secondary mode with the primaries through a subsequent three-wave interaction produces at third order two secularly growing linear Alfvén waves with $\boldsymbol {k}_3 = (1,2,1)$ and $(2,1,-1)$, confirming the prediction that there is no parallel cascade of energy. The analytical results also confirm the constraint that upward- and downward-propagating fluctuations do not exchange energy. In other words, the upward-propagating primary mode transfers energy to the upward-propagating tertiary mode, and the downward-propagating primary mode transfers energy to the downward-propagating tertiary mode. These analytical solutions are fundamentally identical to the Newtonian results of Howes & Nielson (2013) and highlight the fundamental role that three-wave interactions of Alfvén waves play in transferring energy from large to small scales in relativistic plasmas.

Since there is not a formally incompressible limit in relativistic systems wherein there is a finite speed of propagation, $c$, we turn to numerical simulations to confirm these analytical results and test the coupling of the Alfvén and fast modes in relativistic weak turbulence. We perform a set of numerical simulations for $\sigma \in [0.01;0.1;1;7;20;\infty ]$ while keeping $\chi = 0.01$ fixed to maintain weak turbulence. We examine both Alfvén wave–Alfvén wave collisions as well as Alfvén wave–fast wave collisions in both cubic and elongated domains, where the elongated domain is employed to approximate the reduced limit of $k_\parallel \ll k_\perp$. The numerical results confirm the analytical findings for the case of Alfvén wave–Alfvén wave collisions in both cubic and elongated domains, and the Alfvénic cascade is unaltered by the elongated domain. We find that in the cubic domain, both the Alfvén wave–Alfvén wave and Alfvén wave–fast wave collisions produce secondary and tertiary fluctuations consistent with fast waves. However, in the elongated domains, the Alfvén wave–fast wave interaction is suppressed, and the suppression is proportional to the elongation factor, thus confirming our theoretical expectation that the two modes decouple in the anisotropic limit. More detailed numerical analysis of relativistic Alfvén wave and Alfvén wave packet collisions, including the formation and evolution of current sheets, can be found in Part 2 (Ripperda et al. 2021).

The results of this analytical and numerical study of weak, relativistic, Alfvénic turbulence present a simple picture of nonlinear energy transfer through Alfvén wave collisions by highlighting the importance of the three-wave interaction and the decoupling of the fast and Alfvén modes in the anisotropic, reduced limit. Further, they demonstrate the fundamental importance of Alfvénic turbulence to high-energy astrophysical systems, extending the work of Ng & Bhattacharjee (1996), Galtier et al. (2000), Chandran (2005) and Howes & Nielson (2013). With these insights into the continued importance of Alfvénic interactions even in the relativistic limit, we may examine turbulence in these extreme astrophysical systems with newfound intuition.