Weak Alfv\'{e}nic turbulence in relativistic plasmas. Part 1. Dynamical equations and basic dynamics of interacting resonant triads

Alfv\'{e}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 MHD equations and employ them to examine weak Alfv\'{e}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 the relativistic system. 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\'{e}n coupling is diminished in the limit of oblique propagation.


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 h, where b 2 = b µ b µ = B 2 /γ 2 + (B · v) 2 is the magnetic energy density, h is the enthalpy density, and γ 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 socalled E × B nonlinearity is the dominant nonlinear term in the plasma. In terms of fluid, MHD, equations, this term appears as v ⊥ · ∇ , where v ⊥ E × B 0 /B 2 0 is the dominant drift velocity in the plasma, B 0 is a mean magnetic field, indicates either v or B, and perpendicular, ⊥, indicates perpendicular to 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 ⊥ = 0. This term also requires that the perpendicular wave vectors of interacting Alfvén waves are not collinear, but this point is not obvious from the simplified expression above and will be explored further in §2.1. Although there are cases in which the E × B nonlinearity is not dominant, e.g., parametric instabilities such as the decay, modulation, and beat instabilities, we will adopt the assumption that the E × 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 λ mf p , where l is an intermediate turbulence scale and λ mf p 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 E × 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 ∼ k ⊥ at all scales. Montgomery & Turner (1981) and Shebalin et al. (1983) revisit the IK theory under the assumption that the turbulence is weak, ω N L ω, where ω N L ∼ k ⊥ δv ⊥ is the nonlinear frequency and ω ∼ k 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 E × 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 = 0 mode, leading to an anisotropic cascade of energy which fixes k at the outer-scale. Based on this interaction, Montgomery & Turner (1981); Montgomery (1982); Higdon (1984) develop a theory of highly anisotropic incompressible MHD turbulence consisting of twodimensional 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); Strauss (1976) for systems in which there is a strong mean magnetic field, leading to the ordering assumptions of RMHD: δB ⊥ /B 0 ∼ δv ⊥ /v A ∼ k /k ⊥ ∼ 1. Note that the assumptions of RMHD necessarily imply that the Alfvén and fast modes are well separated in frequency since ω A ∝ k ω F ∝ k ∼ k ⊥ . 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 will revisit the concept of reduced MHD in the relativistic limit in section §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 ω = 0, which cannot be a linear wave mode with finite power. Therefore, they invoke a four-wave interaction (A + A → A + A) which maintains k = 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 χ = ω N L /ω 1, where χ is the nonlinearity parameter and χ 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 ∼ k −5/3 ⊥ and a spectral anisotropy of the form k ∼ k 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 k ⊥ and δB 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 = 0 modes are valid nonlinear fluctuations. Ng & Bhattacharjee (1997) further employ perturbation theory to demonstrate explicitly that the threewave interaction is non-empty, k = 0 modes do develop, and the three-wave interaction dominates over the four-wave. They also present the energy spectrum of weak turbulence, E ∼ k −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 ∼ k −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 −2 ⊥ scaling result based on the wave-kinetic approach (Zakharov et al. 1992), and Perez & Boldyrev (2008) numerically verify the derivation and examine the weak to strong turbulence † Note that in the purely 2D limit, i.e., the 2D plane not containing the mean field, the distinction between weak and strong Alfvénic turbulence vanishes, since the linear physics of Alfvén waves is mostly eliminated. Therefore, 2D simulations in this geometry are always strong, analogous to hydrodynamic turbulence that does not support linear waves.
‡ Built into this statement and all discussions of turbulence herein is the assumption that the viscosity and resistivity are sufficiently small that a turbulent cascade is able to fully develop. 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 , numerically , and in an experiment Drake et al. 2013) conducted at the Large Plasma Device (LaPD) (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 = 0 mode is found to be dominant. This mode does not involve a secular exchange of energy and oscillates with frequency ω = 2ω 0 , where ω 0 = k v A is the frequency of the incident Alfvén waves. At third order, the incident Alfvén waves interact with the k = 0 mode, and the k = 0 mode shears the Alfvén waves in the perpendicular plane, providing a secular exchange of energy to smaller perpendicular scales but fixed k scale.
The theory of critically balanced strong turbulence has also been further refined since Goldreich & Sridhar (1995). Boldyrev (2005Boldyrev ( , 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;Comisso et al. 2018;Dong et al. 2018). Intermittency has also been built into the theory of critically balanced and dynamically aligned turbulence , leading to the theory of refined critical balance . 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) ;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 = 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 A + A → F, A + F → A, or F + F → 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 et al. (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 fields 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 Newtonian MHD turbulence, except they consistently find a small portion ( 10 − 15%) of the initial Alfvénic energy leaks into the fast mode branch (Cho 2005;Takamoto & Lazarian 2016Li et al. 2019).
This work is organized as a sequence of papers. In this manuscript (Paper I), we derive a set of relativistic reduced MHD (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  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 with 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 (Paper II) 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 to 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 the second paper of the sequence (Ripperda et al. 2021), henceforth Paper II.

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 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 ∇ · v = 0, leading to the incompressible MHD equations where P is the total (thermal plus magnetic) pressure, and ρ 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 where the magnetic field has been separated into equilibrium and fluctuating parts, B = B 0ẑ + δB, v A = B 0 / √ 4πρ 0 , and z ± = v ± δB/ √ 4πρ 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 Eq. (2.4) to find (2.6) The Elsasser equations, Eq. (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 ω = ±k v A : (i) Alfvén waves with z ± polarized in theẑ ×k ± direction, and (ii) pseudo-Alfvén waves, the incompressible limit of magnetosonic slow modes, with polarizationk ± × ẑ ×k ± , wherek ± corresponds to the unit vector along the wavevector of z ± . We adopt the convention that ω 0, and the sign of k determines the propagation direction. Note that the incompressible assumption orders the magnetosonic fast mode out of the system by imposing c s → ∞. Thus, one can interpret the Elsasser field z + (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, z ∓ · ∇z ± . The most immediate point one can see from the form of this term is that the nonlinearity only survives if both z + and z − are nonzero, 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 z − = 0, then 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 −ẑ 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 z + wave distorted by a z − wave is where we have used the polarization properties of the Alfvén waves to write z ± ∝ẑ × k ± . Therefore, for the nonlinear interaction of counter-propagating Alfvén waves to be nonzero,ẑ · k− ×k + = 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 wave vector 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 = 0 to propagate; (ii) The nonlinearity requires counter-propagating fluctuations; (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 et al. 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 z ± with Eq. (2.4) and integrating over all of space. Assuming periodic boundary conditions or that the fields vanish at infinity, one obtains Eq. (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.

The Connection to Reduced MHD
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 k ⊥ (and δB/B 0 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 et al. 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 reduced MHD (RMHD) are anisotropic fluctuations (k ∼ k ⊥ ), small amplitude fluctuations relative to the background (e.g., δB ∼ B 0 ), and characteristic timescales ω ∼ k v A , where 1. Note that early derivations of RMHD (Kadomtsev & Pogutse 1974;Strauss 1976) further assumed strong magnetization, implying plasma β 1. However, more recent derivations (Schekochihin & Cowley 2007;Schekochihin et al. 2009) have demonstrated that the RMHD equations are valid for arbitrary plasma β, assuming a homogeneous background. Conveniently, the equations of RMHD are essentially identical to Eq. (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 favorable properties compared to 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 ⊥ ρ i 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 Eq. (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 ⊥ ρ i 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 to linear terms in RMHD, we find 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 (χ 1) and critically balanced (χ ∼ 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 χ > 1. Physically, critical balance amounts to equating the linear, propagation time, τ A = l /v A , with the nonlinear (or "eddy turnover") time, τ N L ∼ l ⊥ /u ⊥ . The case τ A τ N L (χ 1) is the weak turbulence case, which will eventually, at sufficiently small scales, transition to τ A ∼ τ N L because the weak turbulence cascade does not produce a cascade in the direction parallel to the equilibrium field: l is fixed at the outer-scale. Thus, in weak turbulence τ A remains fixed while τ N L will decrease with scale. The case τ A τ N L corresponds essentially to creating two dimensional structures perpendicular 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 Thus, if a system is driven such that τ A τ N L , it will rapidly relax back to critical balance, τ A ∼ τ N L at small scales. It is also worth noting that the critical balance condition, or predictions that follow from critical balance, have been observed in simulations (Cho & Vishniac 2000;Maron & Goldreich 2001;Perez & Boldyrev 2008;Mallet et al. 2015;, 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 x 0 , B local 0 (x 0 ) = B 0 + l l δB l (x 0 ). Since the propagation and nonlinear times both decrease with scale, the fluctuations with l l are approximately static relative to small scale turbulence fluctuations, and those with l l are more rapid and will average to zero. Therefore, small scale Alfvénic turbulence always sees an approximately constant, local mean magnetic field.

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 (ω F ∝ k), and therefore the fast mode turbulence cascade is also isotropic, as confirmed in numerical simulations of compressible MHD (Cho & Lazarian 2002. In the weak limit, Chandran (2005); Luo & Melrose (2006); 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 > k ⊥ ), the frequency and wavevector matching conditions permit any of the following three-wave interactions on approximately equal footing: A + A → A, A + A → F, A + F → F, and F + F → F. However, in the obliquely propagating limit (k k ⊥ ), the cascades decouple, leaving only A + A → A and F+F → 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 et al. 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 waveparticle interactions (Landau 1946;Barnes 1966) for a wide range of plasma parameters . Therefore, in the Newtonian limit, focusing only on the Alfvénic turbulence cascade is, generally, well justified.

Relativistic Elsasser Equations
To begin our exploration of weak turbulence in relativistic, magnetically dominated plasmas, we would like to start from an equation set that: (i) makes direct contact with the Newtonian Elsasser equations, Eq. (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 et al. (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 will employ Lorentz-Heaviside units and all speeds are normalized to the speed of light, c = 1, and we will assume a fixed Minkowski metric with signature η µν = diag {−1, 1, 1, 1}. The ideal relativistic MHD equations are mass conservation, the stress energy equation, and the induction equation, where ρ is the rest-mass density, u µ = (γ, γv) is the four-velocity, γ = 1/ √ 1 − v 2 is the Lorentz factor. Furthermore, is the stress-energy tensor, are the relativistic Elsasser four-vector fields, and Much like in non-relativistic MHD, Eq. (3.6) represents the evolution of upward and downward propagating Elsasser fields; however, Eq. (3.6) represents the fully compressible system and is thus not completely analogous to Eq. (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 reduced MHD to isolate the Alfvénic fluctuations described by Eq. (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 et al. 2003), coronae (Chandran et al. 2018;Yuan et al. 2019), and magnetar magnetospheres (Parfrey et al. 2012). Further, as we will 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 (Bruno & Carbone 2013;Howes et al. 2008;Kiyani et al. 2015;Chen 2016).

Relativistic Reduced MHD
To derive a set of Elsasser equations for relativistic reduced MHD (RRMHD), we begin by separating the fluid into mean (background) quantities and fluctuating quantities of the form, (3.9) We also construct an average fluid rest frame in which u i vanishes ‡. We define λ to be the correlation length transverse to B i , and L to be the characteristic scale of the background, i.e., the length scale parallel to B i . As in RMHD, we assume λ/L ∼ O ( ), where 1. In addition to the anisotropy assumption, we also assume that the characteristic frequency is of order the Alfvén frequency, ∂ t ∼ B i ∂ i ∼ k v A , and that the fluctuations are ordered small: and v A = v i A v Ai . We will further assume that all background quantities, e.g., B i = B 0 , ρ ≡ ρ 0 , p ≡ p 0 , are constant with no background inhomogeneities.
As in Newtonian RMHD, in relativistic RMHD, the fast wave will be ordered out of the system, since ω F ∼ k ω ∼ k . RRMHD 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 will also assume B 0i δz i ± = O 2 to remove the pseudo-Alfvén modes. With all of the above assumptions, we note that γ Applying the above reduced assumptions to the relativistic Elsasser equation, Eq. (3.6), we arrive at where to lowest non-zero order 14) 16) † The range of scales well separated from driving, dissipation, and kinetic effects. ‡ Note that u µ is not formally a four-velocity, since u µ uµ = −(1 + δv 2 ). E 0 = 4p 0 + ρ 0 + B 2 0 , and δE = 4δp + δρ + 2B 0 δB . In three-vector form, the equation set is particularly simple and similar to the Newtonian Elsasser equations, The terms of the left-hand side of Eq. (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 Eq. (3.13) are the two components transverse to the mean magnetic field, as expressed in Eq. (3.17). The other two components, time-like and parallel, are one order higher in the ordering parameter, . Importantly, the parallel fluctuations do not appear at lowest order in the nonlinear E × B term in the perpendicular, Alfvén wave equations, Eq. (3.13) and Eq. (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 δz ⊥± · ∇ ⊥ δz . 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 Eq. (3.17), or four-divergence of Eq. (3.13), to find an equation for δΠ.
The reduced relativistic Elsasser equations, Eq. (3.17), for relativistic RMHD 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 Eq. (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, σ = b 2 /h 1. In this limit, Π = 1/2, and δΠ ∼ O 2 /σ O 2 . Therefore, we arrive at the magnetically dominated, relativistic, reduced MHD equations v ν A∓ ∂ ν δz µ ± = −δz ν ∓ ∂ ν δz µ ± . (3.18) In three-vector form, the equation set is particularly simple, Being in the σ 1 limit, Eq. (3.18) and Eq. (3.19), are essentially the reduced analog 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: (iii) Small amplitude fluctuations with constant (mean) backgrounds; (iv) Second order parallel and pressure fluctuations; (v) Magnetically dominated, σ 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 k ⊥ ). 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 σ → ∞ 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 et al. 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 k ⊥ , which naturally arises as the Alfvénic turbulence cascades to small scales.

Connection and Comparison to Newtonian Reduced MHD Solutions
Eq. (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 = 0 to propagate; (ii) The nonlinearity requires counterpropagating fluctuations; (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 these papers, the wave kinetic equation evolves the spectral energy e ± = |δz ± ⊥ (k)| 2 and can be expressed as 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 )k −2 ⊥ , where f (k ) 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, Eq. (3.17), to derive the turbulence solutions in the following section, the Alfvén solutions are identical for the σ → ∞ limit of the RRMHD equations, Eq. (3.19).

Three-wave Weak Alfvénic Interactions
To construct analytical solutions, we now consider a subsidiary expansion of the form ζ ± = εζ 1± + ε 2 ζ 2± + · · · , where δz ± =ẑ × ∇ ⊥ ζ ± defines the Elasser potential, ζ, and ζ i± (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 §A. Therefore, we primarily summarize the solutions found in ; however, we note that we include corrections to  Eqs. (22), (28)(29), and all equations involving ζ 3− , including the equations for B ⊥3 and E ⊥3 . These errors have been corrected in the §A, and the changes are highlighted in red.
For specificity, we will consider a periodic domain of size L At t = 0, we will initiate two counter-propagating, perpendicularly polarized fluctuations where z ± are the initial amplitudes, k ⊥ = 2π/L x = 2π/L y , and k = 2π/L z . We maintain the convention that k ⊥ , k , and ω 0 are positive constants, and the direction of propagation is supplied by the explicit sign of k . To describe the wave modes that arise from the nonlinear evolution, we will use the notation (k x /k ⊥ , k y /k ⊥ , k z /k ). 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 §A. Here, we summarize the important findings at each order. The initial conditions provided in Eq. (4.1) satisfy the lowest order equations if ω 0 is the linear Alfvén wave frequency, ω 0 = k v A . Thus, at lowest order the solution described counter-propagating, linear Alfvén waves.
The second order solutions for the electric and magnetic field are given by Eqs. (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 k 2 = k 1− ± k 1+ = (0, 1, 1) ± (1, 0, −1). The (−1, 1, 2) mode satisfies the conditions for a linear Alfvén wave. Specifically, the frequency of this mode is ω = ±2ω 0 = ±2k v A = ±k z v A , and the mode obeys the Alfvén wave eigenrelation 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 ω = 2ω 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 B ⊥ and E ⊥ , we can once again straightforwardly interpret the third order solutions, Eqs. (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 ω = ±ω 0 = ±k z = ±k z z and obey the Alfvén wave eigenrelations, Eq. (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 −π/2, have perpendicular wavevectors √ 5k ⊥ , but |k z | = k . Therefore, the weak turbulence cascade proceeds to smaller scales in the perpendicular plane, but the parallel scales are conserved, confirming the prediction from simple threewave 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: k 3 = k 2 + k 1± = (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 non-linear structures. Unlike the second-order solutions, there is not a 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 √ 5k ⊥ , but these all have both linear and non-linear components. Similarly, the (0, 1, 1) and (1, 0, −1) components appear as both linear and non-linear 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.

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 3D, periodic domain along a uniform guide field B 0 = B 0ẑ using the general relativistic MHD code BHAC (Porth et al. 2017;Olivares et al. 2019;Ripperda et al. 2019a,b). The recently added force-free limit for the resistive MHD code BHAC employs the numerical scheme of Ripperda et al. (2019a,b) and damps force-free violations E > B and E·B = 0 on resistive time scales †. We employ both a period cubic domain with L ⊥ = L x = L y = L = L z = 2π and a periodic elongated domain with L = 10L ⊥ = 20π, with resolution N x = N y = N z = 256 and N z = 2560 cells for the elongated case. Initially, we set a gas-to-magneticpressure ratio of β = 2p/B 2 0 = 0.02 and magnetization σ cold ≡ b 2 /ρ ∈ [0.01; 0.1; 1; 10; 100] (corresponding to σ ∈ [0.01; 0.1; 1; 7; 20]), where we vary the density, ρ, but maintain a constant guide field, B 0 =ẑ, and pressure p = 0.01 with adiabatic index Γ = 4/3 for an ideal relativistic gas. In the force-free case, β → 0, σ → ∞, and v A → 1.
As in §4, we prescribe a scale-free definition of characteristic wavelengths (k x /k ⊥ , k y /k ⊥ , k z /k ), and initialize our Alfvén wave simulations with two overlapping, counterpropagating, and perpendicularly polarized Alfvén waves described by the wavevectors k + = (1, 0, −1), and k − = (0, 1, 1). The magnetic field is initialized through a vector potential A = (−B 0 y, 0, δB ⊥ [sin(k ⊥ x + k z) + sin(k ⊥ y − k z)]), representing the initial counter-propagating Alfvén waves with frequency ω 0 = k v A . The electric field is initialized as E = (v A B y , v A B x , 0) such that the velocity is equal to the drift velocity, δu ⊥ = E × 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 E ∓ · B ± = 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 χ = k ⊥ δB ⊥ /k B 0 . To maintain weak turbulence, we fix χ = 0.01 in all of the following simulations. By fixing χ, we expect to maintain self-similar behavior as we explore elongated domains that approximate the reduced limit. Thus, in the cubic domain with k ⊥ = k , δu ⊥ /v A = δB ⊥ /B 0 = 0.01, and in the elongated case with k ⊥ = 10k , δu ⊥ /v A = δB ⊥ /B 0 = 0.001.

Cubic Domain
In Fig. 1, we present the evolution of the mode amplitudes of the B ⊥ (Alfvénic) fluctuations in a cubic domain scanning σ cold ≡ b 2 /ρ ∈ [0.01; 0.1; 1; 10; 100; ∞], presented in descending order of σ from top-left to bottom-right. The initial, counter-propagating, and perpendicularly polarized Alfvén waves are represented by the red lines in each figure, † Full details of the simulation code, including convergence, numerical diffusion, and numerical dispersion studies, can be found in Paper II (Ripperda et al. 2021).
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 ω = 2ω 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 ⊥ Alfvén waves, but with the same k as the primary waves, where the (1, 2, −1) and (2, 1, 1) modes are represented by the blue dashed and dotted curves. These dynamic are in agreement with the results of §A.4. These modes grow secularly in time, B ⊥3 ∝ 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 ⊥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.

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 2016Li et al. 2019). Through a numerical simulation study of relativistic turbulence, Takamoto & Lazarian (2016 find that the fast-to-Alfvén mode power scales as (δu f /δu A ) 2 ∝ (1 + σ)δu A /c f ⊥ in the σ 1 limit, where c f ⊥ is the fast mode speed in the perpendicular plane. In the σ 1 limit, Cho & Lazarian (2002 find that (δu F /δu A ) 2 ∝ δu A /c f ⊥ . Thus, for σ 1, the fast and Alfvén branches decouple as the background magnetic field strength is increased; however, for σ 1, the coupling between the modes increases with √ 1 + σ, asymptotically approaching unity for large σ. In Fig. 2, we present the evolution of the highest amplitude B modes for the same initial configuration as presented in Fig. 1. B serves as a proxy for the amplitude of compressible, fast mode fluctuations, since the Alfvén modes have negligible B 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.
We first note that the (0, 1, −1) and (1, 0, 1) B modes have zero amplitude at t = 0 but do develop finite amplitude that is approximately independent of σ. These modes appear at a mixture of two or more frequencies, the most dominant of which is at the fast mode frequency, ω ∼ ω f ∼ kc f ∼ √ 2ω 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 ω ∼ ω f ∼ kc f ∼ √ 2ω 0 and ω ∼ ω 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 ω ∼ ω f ∼ kc f ∼ √ 6ω 0 . The tertiary mode also displays the most significant dependence on σ. For σ 1, the amplitude is independent of σ.

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  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. for a cubic domain with χ = 0.01 and σ cold ≡ b 2 /ρ ∈ [0.01; 0.1; 1; 10; 100; ∞], presented in descending order of σ from top-left to bottom-right. Line styles and colors are as defined in Fig. 1. that may naturally arise from a strong magnetic field. As noted in §2.3, in the Newtonian, obliquely propagating limit (k k ⊥ ), 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 = 10L ⊥ .
In the top two panels of Fig. 3 are presented the mode analysis results of a force-free simulation of obliquely propagating Alfvén waves, with k ⊥ = 10k . The color and line styles are as in the previous figures. Comparing the relative amplitudes between each of the B ⊥ modes to the relative amplitudes for the cubic domain case presented in the upper left panel of Fig. 1, it is clear that the elongated domain does not change the Alfvénic B ⊥ , cascade. However, the three parallel modes' growth, amplitude, and frequency in the elongated case compared to the cubic case in the upper left panel of Fig. 2 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 (red) and third (blue) order parallel modes are dominated by low frequency oscillations 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.
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, ω A = k 0.1ω F = 0.1k, 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 k ⊥ .

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 counterpropagating fast mode. Therefore, we choose a fiducial fast mode with k = (0, 1, −1), and an Alfvén wave with 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 will continue the convention of fixing χ = 0.01 to satisfy weak turbulence.

Cubic Domain
In the middle two panels of Fig. 3 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 ⊥ panel correspond to the initial Alfvén and fast waves respectively. Focusing on the behavior in the middle left B ⊥ panel, 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 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 the upper left panel of Fig. 1, this tertiary mode is decreased in amplitude by approximately a factor of three.
Turning to the middle right B panel for the Alfvén wave-fast wave in Fig. 3, we see different behavior compared to the upper left panel of Fig. 2. 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 "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.

Elongated Domain
In the bottom two panels of Fig. 3 are plotted the results of the Alfvén wave-fast wave collision outlined above for the elongated domain with L = 10L ⊥ . The overall behavior is similar to the Alfvén wave-fast wave collision in the cubic domain with a few notable differences. Focusing first on the B ⊥ panel in the lower left of Fig. 3, 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 the upper left panels of Fig. 1 and Fig. 3. 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 lower right B panel for the Alfvén wave-fast wave in Fig. 3, 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 higher frequency, and the "primary" mode is the only B 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 wavefast wave collision in both cubic and elongated domains extend the findings of to the relativistic limit.

Summary & 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 k ⊥ . Leveraging the anisotropic cascade of energy and the fact that δB/B 0 1 is satisfied either at the outerscale or once the cascade reaches sufficiently small scales, we derive a set of relativistic reduced MHD (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 σ = b 2 /h → ∞. However, we do also derive a set of magnetically dominated RRMHD equations which are valid in the σ → ∞ limit but retain the favorable 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 turbulence 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 )k −2 ⊥ , where f (k ) 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 ω = ω 0 and wavevectors k 1± = (1, 0, 1) and (0, 1, −1) that interact through a three-wave interaction to produce a non-linear, magnetic shear mode, k 2 = (1, 1, 0), with ω = 2ω 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 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. These analytical solutions are fundamentally identical to the Newtonian results of  and highlight the fundamental role 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 σ ∈ [0.01; 0.1; 1; 7; 20; ∞] while keeping χ = 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 k ⊥ . 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 Paper II, (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); . With these insights on the continued importance of Alfvénic interactions even in the relativistic limit, we may examine turbulence in these extreme astrophysical systems with newfound intuition.

A.1. Conversion to Elsasser Potentials and Characteristic Variables
To prepare the system of equations for analysis of their asymptotic solutions, we will convert to the Elsasser potential formulation. This approach eliminates the nonlinear pressure term, and reduces the equation system to two scalar equations for the Elsasser potentials, ζ ± rather than two vector equations for the Elsasser fields, z ± . † The Elsasser potentials are defined by the relation δz ± =ẑ × ∇ ⊥ ζ ± , which follows from the fact that the Elsasser fields are solenoidal. Note that δv ⊥ and δB ⊥ can be reconstructed from the potentials Also, using the fact that in general in relativistic MHD we use E = −v × B, to lowest order Taking the curl of Eq. (3.17) and substituting the expression for the Elsasser potentials yields the Elasser potential equations which are identical in form to those derived by (Schekochihin et al. 2009). The Poisson bracket is defined by Eq. (A 3) retains the form that the left-hand side describes the linear evolution, and we can further simplify the task of solving the equation set by converting to characteristic variables, φ ± = z ± v A t ‡. In terms of these characteristic variables, Eq. (A 3) becomes Finally, the Elsasser potential form for the initial conditions provided in Eq. (4.1) is where the final equality follows from conversion to the characteristic variables φ ± .
A.2. Linear, O (ε), Solutions At lowest, linear, order, Eq. (A 5) reduces to The initial conditions above satisfy Eq. (A 7) if ω 0 is the linear Alfvén wave frequency, ω 0 = k v A . Thus, at lowest order the solution describes counter-propagating, linear Alfvén waves, as expected of RRMHD.