Nonlinear electromagnetic-wave interactions in pair plasma: (I) Non-relativistic regime

This paper is the first in a series devoted to the numerical study of nonlinear interactions of electromagnetic waves with plasma. We start with non-magnetized pair plasmas, where the primary processes are induced (Compton) scattering and the filamentation instability. In this paper, we consider waves in which electron oscillations are non-relativistic. Here, the numerical results can be compared to analytical theory, facilitating the development of appropriate numerical tools and framework. We distill the analytic theory, reconciling plasma and radiative transfer pictures of induced scattering and developing in detail the kinetic theory of modulation/filamentation instability. We carry out homogeneous numerical simulations using the particle-in-cell codes EPOCH and Tristan-MP, for both monochromatic waves and wave packets. We show that simulations of both processes are consistent with theoretical predictions, setting the stage for analyzing the highly nonlinear regime.


INTRODUCTION
The interaction of intense radiation with plasma is of great interest both for laboratory physics and for astrophysics. In the laboratory, the main interest is in studying the propagation of short, narrow laser pulses through a dense electron-ion plasma, in which case the main effects are self-focusing, induced Raman scattering, and the generation of wake-fields due to charge separation (see, e.g., the review by Mourou et al. 2006). In astrophysical systems, the spatial and time scales of the radiation are macroscopically large; even the millisecond pulses may be considered as macroscopic as compared with the wave period. In this case, self-focusing yields to filamentation of the radiation flux: the wide beam splits to narrow subbeams. In most astrophysical applications, the radiation frequency is well above the plasma frequency, in which case the non-linear effects become weaker. Nevertheless, even weak non-linear corrections could accumulate considerably over a long enough path, thus significantly affecting the propagation of the wave. In the high frequency, low density regime, induced Raman scattering becomes weaker than induced Compton scattering. Moreover, in many of the systems of interest, the radiation propagates through an electron-positron (henceforth pair) plasma, rendering charge separation effects, including Raman scattering, irrelevant. Therefore, for astrophysical applications, induced Compton scattering and the filamentation instability are the most important processes.
In the modern plasma literature, induced Compton scattering is sometimes called stimulated Brilluin scattering (e.g., Edwards et al. 2016;Schluck et al. 2017). The term Brilluin scattering implies a three-wave interaction of electromagnetic and sound waves. The latter do not exist in non-magnetized, collisionless plasmas (with the exception of electron-ion plasmas with electrons much hotter than ions), in which the electromagnetic waves are scattered off density fluctuations. This † Electronic address: arka@post.bgu.ac.il process is called the induced Compton scattering both in classical plasma works (e.g., Litvak & Trakhtengerts 1971;Galeev & Syunyaev 1973;Drake et al. 1974;Ott et al. 1974;Lin & Dawson 1975) and in the astrophysical literature. Therefore, we use this term (henceforth referred to simply as induced scattering), retaining the term Brilluin scattering for scattering off ion sound waves in two-temperature plasmas or off magnetosonic waves in magnetized plasmas.
The filamentation instability attracted much less interest in the astrophysical community, presumably because this effect is unable to prevent the escape of the waves. Nevertheless, the instability may affect the properties of the outgoing radiation (Sobacchi et al. 2020). Namely, the filamentation instability implies the formation of narrow subbeams, which are strongly diffracted when the radiation escapes the plasma slab. This could lead to smearing of short bursts in time and/or to frequency modulations of the observed intensity due to interference between the subbeams. Filamentation was observed in numerical simulations of relativistic, magnetized shocks, in which case a strong electromagnetic precursor is emitted upstream (Iwamoto et al. 2017;Sironi et al. 2021).
Our goal is to study the observational features of the above nonlinear processes, by simulating numerically the propagation of intense radiation through a pair plasma. In the case of induced scattering, this opens a possibility not only to estimate the optical depth, but also to find the characteristics of the outgoing radiation. In the case of the filamentation instability, we are interested in the quantitative parameters of the sub-beams, which can facilitate quantitative predictions for the timefrequency structure of the observed radiation. Among other sources, these effects are of special interest for FRBs. Their extremely high brightness temperature implies that close enough to the source, the electrons oscillate relativistically in the field of the FRB radiation (Luan & Goldreich 2014). A quantitative study of the above effects in such an extreme case is possible only numerically.
In this paper, we develop the setup for particle-in-cell (PIC) simulations of the propagation of intense radiation through pair plasmas, and demonstrate both induced scattering and the filamentation instability. Both processes may be considered as an instability of the pumping beam. Here, we consider only the case where the electron oscillations in the field of the wave are not relativistic. In this regime, the linear stage of the instability may be considered analytically, so we can verify the simulation results. The analysis of the highly nonlinear evolution, including the case of truly strong waves, in which electron oscillations become relativistic, is deferred for a future publication. We present the theoretical background in §2, the PIC simulation setup in §3, and the results in §4. We summarize and discuss our findings in §5. We use Gaussian cgs units with k B = c = 1.

INTERACTIONS WITH PAIR PLASMAS
The propagation of electromagnetic waves is conveniently studied using the wave equation for the vector potential A, In this paper, we consider a nonrelativistic pair plasma.
Here, the current is presented as where the summation is over the species, q = ±e, e is the electron charge, N is the density of the species, and v is the particle velocity. The particle oscillations in the field of the wave are found from the equation of motion, where m is the electron mass. In this paper, we assume that the oscillation velocity is nonrelativistic, v ≪ 1. Then to first order in the small parameters v and eA m, The unperturbed plasma is assumed to be neutral, with a particle number density N tot = 2N 0 , where N 0 is the unperturbed electron density. In the first approximation, Eq.
(2) then becomes where is the plasma frequency squared. Now the wave equation may be written as where j nl is the nonlinear current, which could arise from non-linear corrections to the electron velocity or/and from variations of the plasma density caused by the wave. In our approximation, the nonlinear current is small, so the nonlinear corrections are small at the wavelength scale. However, when the nonlinear current is in resonance with the wave, the corrections accumulate, and so could lead to a significant modification of the initial wave at longer time scales.
2.1. Induced scattering Even though the term "induced scattering" refers to the quantum picture, the effect is purely classical: the beating between the initial and the scattered waves produces a ponderomotive force, which gives rise to a plasma density modulation, which in turn leads to an enhanced scattering rate. With such an approach, the process can be considered as an instability: when the plasma is illuminated by a pumping radiation beam, the amplitudes of both the scattered wave and the density modulation grow. The derivation of the induced scattering of electromagnetic waves in plasma was given by, e.g., Litvak & Trakhtengerts (1971) and Drake et al. (1974). Here, we outline the derivation in the simplest case of a pair plasma, and find the condition under which the plasma theory results reduce to those of radiation transfer theory.
Having in mind 2D simulations, we consider scattering in the xy plane, with the initial wave propagating along the x-axis and the scattered wave being offset by an angle θ, assuming both waves being polarized in the xy plane (see illustration in Fig. 1; the generalization to a general setup is discussed farther below). Then the electromagnetic vector potential may be written as where indices 0 and 1 refer respectively to the initial and scattered waves, A 0 and A 1 are in the xy plane, c.c. is the complex conjugate, and r = {x, y, z} are Cartesian coordinates.
The beating wave, of angular frequency and wave vector given respectively by is in Landau resonance with particles moving with the phase velocity of the wave, v = ω d k d k 2 d . These particles "see" a force constant in time that produces a density modulation, δN e ik d ⋅r−iω d t + c.c. , with an amplitude growing in time. Therefore, the rate of scattering on these fluctuations grows, too. Thus the mechanism redistributes the energy of the initial wave into the energies of the scattered wave and the plasma x y (k 0 ) p particles. An important point is that inasmuch as v ≪ 1 and the scattering angle is typically large, k d ∼ k 0 and The density perturbations may be derived from the Vlasov equation, where f is the distribution function, p is the particle momentum, and F is the force exerted on the particles by the beating wave. Inasmuch as the beating frequency ω d is small as compared with ω 0 and ω 1 , one could express such a slowly varying force via the ponderomotive potential (e.g., Schmidt 1979), where the angular brackets denote averaging over a time period longer than ω −1 0 but short as compared with ω −1 d . Substituting Eq. (8) in Eq. (12) yields ϕ = const.
and µ ≡ cos θ. The ponderomotive potential is a quadratic function of the particle charge, therefore the electron and positron densities remain equal even in the presence of fluctuations. Looking for perturbations of the form (10), we may consider the 1D Vlasov equation for the distribution of particle velocities in the directionk d ≡ k d k d . One can conveniently normalize the 1D distribution function f (p) by unity, where p ≡ p ⋅k d , and integrals are henceforth carried out over their full parameter range. In the oscillating ponderomotive potential, the initial particle distribution function, f 0 , gets an oscillating perturbation of the form δf exp(ik d ⋅r−iω d t)+c.c. Linearising Eq. (11) in the small parameters δf f and ϕ , where ≡ p 2 2m is the particle energy, one gets The density perturbations are now found as where the Landau rule for bypassing the pole was used, with i0 designating an infinitesimal positive imaginary constant. Now we can find the evolution of the scattered wave.
The wave equation (7) for the scattered wave becomes where j nl 1 is the component of the non-linear current having the angular frequency ω 1 and the wave vector k 1 , and A 1 ≡ A 1 A 1 is the unit vector in the direction of A 1 . The nonlinear current is produced by beating between the low frequency density fluctuations (10) and the velocity oscillations in the field of the pumping wave. According to Eq. (4), these velocities are just v Substituting this current into Eq. (17), one obtains the equation for the scattered wave, Plugging in δN from Eq. (16), the wave equation turns into the dispersion equation Instead of the incident wave amplitude, one can conveniently use the dimensionless strength parameter where the factor 2 appears because our definition (8) implies that the wave amplitude is 2A 0 . According to Eq. (4), the condition that the oscillation velocity in the field of the wave is non-relativistic now becomes a ≪ 1. Using a, the dispersion equation becomes The scattering rate is determined by the imaginary part of ω 1 . The wave ω 0 is scattered into ω 1 , so that A 1 grows, if Im ω 1 > 0. In the opposite case, A 1 decreases as the wave ω 1 is scattered into ω 0 . The imaginary part of ω 1 is determined by the imaginary part of the integral in the RHS of the dispersion equation, which arises because of the pole in the integrand. Therefore, only particles in Landau resonance with the beating wave are responsible for the effect. Generally, the dispersion equation is solved numerically, as we show later. However, if we consider the growth rate to be small enough, where T is the plasma temperature in energy units, we can use the Sokhotski-Plemelj formula in the RHS of the dispersion equation. Then, approximating Im ω 1 ≪ ω 1 in the LHS, we obtain Im ω 1 = π 8 (ω p aµm) 2 ω 1 As df 0 d < 0, we see that if ω 1 < ω 0 , then Im ω 1 > 0, i.e., the pumping wave ω 0 is scattered into the wave ω 1 . Thus, the wave frequency decreases due to induced scattering. In astrophysical applications, one typically defines the rate of induced scattering as the growth rate of the scattered wave energy, κ = 2Im ω 1 . Substituting the Maxwell distribution for f 0 and taking into account that one finds This result coincides with the induced scattering rate used in radiation transfer theory (Zel'Dovich et al. 1972).
The term µ 2 in the numerator applies to scattering in the plane; it should be replaced by (1 + µ 2 ) 2 for the typical radiation-transfer case of 3D scattering of non-polarized radiation.
Equation (26) indicates that the growth rate is maximal for backward scattering, µ = −1; in fact, this applies for the general case of 3D scattering with arbitrary polarization. Induced scattering outside the emission beam occurs due to a weak background radiation (produced, e.g., by spontaneous scattering). In this case, the energy of the scattered radiation grows exponentially, ∝ exp (κt). Many e-folding times are necessary in order to scatter a significant fraction of the beam energy. Therefore, the beam is scattered predominantly into the states corresponding to the maximal scattering rate. The fastest growing mode has the angular frequency the corresponding wave vector being Here, we denote k 1 =k 0 ⋅ k 1 , so the minus sign indicates backscattering. Substituting Eq. (27) into Eq. (26), expanding around T m → 0 and retaining the leading order term yields the maximal scattering rate The validity condition for the asymptotic estimate (26), given by Eq. (23), now becomes Figure 2 presents the growth rate of the backscattered wave, derived both from the numerical solution to the dispersion equation (22) and from the analytic asymptotics (26). We find a good agreement between the two as long as Eq. (30) is satisfied, which is indeed the case in the figure. More examples are provided below, comparing PIC simulations of the scattering process to the theoretical estimates. The above results were derived for a monochromatic pumping wave. For astrophysical applications, the case of a broad spectrum is more interesting. Assuming that the initial beam is composed of waves with different frequencies and random phases, the electromagnetic potential becomes (cf. Eq. (8)) (31) Generally, the summations should be carried out over k modes, but as the directions of both initial and scattered beams are fixed, suffice to sum over absolute values k, equivalent to summations over frequencies. Now, Eq. (13) for the ponderomotive potential is replaced by which yields the density fluctuations (cf. Eq. (16)) The nonlinear current having angular frequency ω 1 and wave vector k 1 (cf. Eq. (18)) now becomes Averaging this current over random phases using Eq. replaced by ∑ ω A ω 2 . The resulting induced-scattering rate in the slow-growth limit of Eq. (23) is where the dimensionless amplitude, a ω , is related to A ω just as a is related to A 0 ; see Eq. (21). It will be shown below that the induced scattering rate of a wide spectrum beam is weaker than that of a monochromatic beams. Therefore the condition (23) is less restrictive in the wide spectrum case.
If the width of the spectrum is small compared to the width of the particle distribution, ∆ω ω ≪ v T ≃ (T m) 1 2 , scattering proceeds as for a monochromatic beam. In the opposite limit of a broad spectrum, the particle distribution function may be approximated by the delta-function, so the frequency shift in the scattering vanishes, ω 1 → ω. Then one can substitute in Eq. (35) the expression and integrate by parts, which yields where ω is both the incident and scattered frequency. If the total energy in the pumping beam, ∼ ω 2 a ω 2 ∆ω, remains fixed, one sees that the scattering rate decreases with increasing spectral width, as κ ∝ ∆ω −2 .
Astrophysical applications usually use the radiation transfer equation, for the radiation intensity I(ν,n) in frequency ν and directionn. For the induced scattering of non-polarized emission with a wide spectrum, this equation is written as (e.g., Wilson 1982) Here, the time derivative is taken along the ray, and the integral is over solid angle Ω. In order to compare Eqs. (37) and (39), note that in our case, the incident radiation is directed along the x-axis, such thatn ′ ⋅n = µ and the radiation flux is The radiation flux can also be expressed using the wave amplitudes as where the factor 2 in the first equality takes into account both polarizations. Comparing the last two equations shows that in our case, ∫ I(ν, n ′ )dΩ ′ should be replaced by m 2 ω 2 a 2 ω 2e 2 . Then, Eq. (39) reduces to our Eq. (37) if one also substitutes (1 + µ 2 ) 2 by µ 2 , in order to take into account that we considered the scattering of polarized radiation in the polarization plane.

Filamentation and modulation instabilities
Due to the nonlinear interactions of radiation and matter, the refraction index depends on the radiation power, which could lead, under certain conditions, to selfmodulation or/and self-focusing of the radiation beam (see, e.g., Karpman 1975). If the radiation beam is wide, which is typically the case in astrophysical applications, self-focusing leads to filamentation of the beam into a set of subbeams. Note that formally, the filamentation instability is a special case of the modulation instability, when modulation is perpendicular to the direction of the beam. However, the physics and the conditions for filamentation substantially differ from those for modulation in the longitudinal direction. Specifically in a pair plasma, only the filamentation instability is possible, as shown below.
The modulation of the amplitude of the pumping wave could be described by introducing two sidebands satisfying Here, ω 0 and k 0 are the angular frequency and wave vector of the pumping wave, respectively, and ω d and k d are the angular frequency and wave vector of modulation.
In order to study evolution of the modulation, one can consider the non-linear interaction of these waves (Drake et al. 1974;Max et al. 1974), as for induced scattering. Instead, we exploit the fact that in this case, both ω d ≪ ω and k d ≪ k 0 , so we deal with a slow, long-wave modulation of a plane, monochromatic wave; the generalization for a broad spectrum is discussed at the end of this section.
The propagation of the incident wave in plasma is described by the wave equation (7). The non-linear current for a monochromatic wave is found in the high-frequency limit as (Montgomery & Tidman 1964;Sluijter & Montgomery 1965; see also Appendix A) where a 2 ≡ 2⟨(eA m) 2 ⟩ is the dimensionless amplitude of the wave; this definition is equivalent to Eq. (21). Taking into account that the non-linear current is directed along A, Eq. (7) reduces to an equation for a scalar A. For a slowly varying wave, one can still use the monochromatic non-linear current (43), in which one has to substitute the local amplitude of the wave. Therefore, the wave equation takes the form where a is now the slowly varying dimensionless amplitude of the wave. Note that in the field of a modulated wave, the plasma density is expected to be modulated; therefore, one has to use the local density in the plasma frequency. Equation (44) has an exact solution describing a planar, monochromatic wave in a homogeneous medium, where a 0 is a dimensionless constant, with a nonlinear dispersion law Let us consider the stability of this wave with respect to a slow, long-wavelength modulation.
One can conveniently use the complex dimensionless amplitude, a(r, t), deviating from a 0 , such that Then the slowly varying, complex function a(r, t) contains information on both the amplitude and the phase of the wave. One must take into account that in Eqs. (43) and (44), a is the dimensionless real amplitude of the wave, so in the new notations, one has to replace a 2 by a 2 . An important point is that the nonlinear evolution of the wave occurs not only due to the non-linear current (43), but also because the plasma is expelled by the ponderomotive force from the regions of enhanced radiation intensity, which affects the refraction index (Drake et al. 1974). Therefore, we substitute ω 2 p by ω 2 p (1 + δN N 0 ) in the wave equation (44) (but not in the dispersion equation (46), which is written for an unperturbed density), so henceforth ω p denotes the plasma frequency for the unperturbed density. Now the wave equation takes the form where is the group velocity of the plane monochromatic wave, and we assumed, without loss of generality, that the constant a 0 is real and positive. In order to derive the evolution of a small modulation of the plane wave, we express the wave amplitude as a = a 0 + δa(r, t), where δa ≪ a 0 . Linearising the wave equation with respect to δa a 0 and δN N 0 and neglecting a 2 with respect to unity in the fourth term, we get The modulation of the incident wave (45) may be found by substituting into this equation the ansatz The density perturbation, δN , is found, as in §2.1, from the Vlasov equation (15) with the ponderomotive potential This yields Plugging Eqs. (51) and (53) into the wave equation (50) and collecting terms with phases exp , one obtains two coupled algebraic equations for a + and a − , The condition that a ± have non-trivial solutions gives the dispersion equation (55) Generally, this equation is solved numerically. We can however obtain analytical results in certain asymptotic regimes. Taking into account that ω 0 ≫ {ω d , k d , ω p }, one sees that the LHS becomes too large unless ω d ≈ v g ⋅ k d . Let us first consider the case of the filamentation instability, where v g ⋅ k d = 0. Then one can neglect ω 2 d with respect to k 2 d . In the integral, one can drop ω d only under a more restrictive condition, The condition (56) ensures that the characteristic instability time is small as compared with the filament crossing time by an electron. In this case, the filaments grow quasistatically: the ponderomotive force is balanced by the plasma pressure at each moment. In the opposite case, which will be considered below, the ponderomotive force is balanced by the electron inertia; then the growth rate of the instability is independent of the plasma temperature. We first assume that condition (56) is fulfilled (the results will be verified as consistent a posteriori), and later explore the regime where the opposite condition prevails.
Inasmuch as the plasma is nonrelativistic, T ≪ m, the first term in brackets in the RHS of Eq. (55) may be neglected. This term originates from the nonlinear current in the wave equation (the RHS of Eq. (44) and the last term in Eq. (50)), so the nonlinear current does not play a role in the filamentation instability of the pair plasma; the effect is determined only by the ponderomotive force. The dispersion equation now reduces to Hence, the instability develops if k d < (m 2T ) 1 2 ω p a 0 , as ω d becomes purely imaginary. The maximum growth rate, is achieved at One sees that the growth rate of the filamentary instability is of the order of the induced scattering rate (29). This result was obtained under the condition (56), i.e., for Γ = Im ω d ≪ k d v T . Making use of Eqs. (59) and (60), the result is seen to be applicable for In Fig. 3, we present the growth rate of the filamentation instability found both from the numerical solution to the dispersion equation (55) and from the asymptotic formula (58). For this purpose, we use parameters for which the condition (61) for the validity of the asymptotic solution is satisfied. Solutions for a wide range of parameters are presented farther below as we compare the results of PIC simulations of the filamentation instability with the theory. -The growth rate of the filamentation instability, normalized to the incident angular frequency, evaluated by numerical integration of the dispersion relation (55) (solid blue) and according to the asymptotic theoretical approximation (58) (dashed green), which is valid under the condition (61). Parameters: T m = a = 0.00137, ωp ω 0 = 0.107, and a 1D Maxwellian f 0 .
As k d ≫ ω d , one can again neglect the first term in brackets on the RHS of Eq. (55). Hence, the nonlinear current does not play a role, whether or not condition (56) is satisfied. In the present case, the dispersion equation takes the form This bi-quadratic equation always has an unstable solution. The growth rate increases with increasing k d , and One sees that at the growth rate approaches the maximal one, Eq. (64), to within 2.5%. On the other hand, the condition Γ This in fact means, that all perturbations with the wave vectors k d,1 < k d < k d,2 grow with practically the same growth rate (64). Now let us turn to the case of modulation at some angle θ d with respect to the pumping wave. We may again approximate ω d ≃ k d ⋅ v g in the dispersion equation (55), otherwise its LHS becomes too large. When θ is sufficiently far from π 2, such that µ d ≡ cos θ d ≫ T m, one can then expand the integral in the RHS of the equation in small v, Then the dispersion equation becomes The RHS of the dispersion equation is seen to be positive, implying that the instability does not develop in this case. One thus concludes that in a pair plasma, only the filamentation instability develops, so the pumping wave is modulated in the direction perpendicular to k 0 . The main effect in this case is density variations due to the ponderomotive force. In an electron-ion plasma, the density variations are suppressed because of a high ion inertia. In this case, the effect of the nonlinear current comes into play; then both the filamentation and the modulation instabilities could develop (Max et al. 1974;Sobacchi et al. 2020).
An important point is that the ponderomotive force depends only weakly on the radiation spectrum provided it is composed from waves with random phases. In this case, one substitutes A 2 in Eq. (12) by ∫ A ω 2 dω. The variations of the plasma density, which lead to the focusing of the radiation, are produced by the total ponderomotive force, therefore the filamentation instability is determined by the total radiation energy and the characteristic frequency, the growth rate of the instability being only weakly dependent on the shape of the radiation spectrum.
Formally, one has to substitute the ansatz (47) by e m A = Re where a ω0 2 is the spectral distribution of the initial beam, δa(r, t) ≪ 1 describes the common envelope. As we have already seen, the nonlinear current (43) may be neglected when considering the filamentation instability in the pair plasma. This means that one can neglect the a 2 0 term in the dispersion equation for the pumping wave, Eq. (45) and also neglect the term in the RHS of the wave equation (44). Substituting the new anzats in thus modified wave equation and linearizing with respect to small δa and δN N 0 , we get (cf. Eq. (50)) The ponderomotive potential is now written as which yields Eq. (53), in which a 2 0 should be substituted by ∫ a ω 2 dω. An important point is that in the case under consideration, ω d and k d are independent of ω because the envelope is common for the whole pumping beam. Substituting thus obtained δN into Eq. (71), one finds the dispersion equation similar to (55), but without the first term in the brackets in the RHS (which is not important anyway) and with a 2 0 substituted by ∫ a ω 2 dω. Therefore the same substitution should be done in all the obtained results for the filamentary instability.
As shown in §2.1, the induced scattering rate decreases with the width of the spectrum. Therefore, although the growth rates of induced scattering and the filamentation instability are comparable to each other for a monochromatic wave, in the case of a wide radiation spectrum, the filamentation instability develops significantly faster than induced scattering. This conclusion is valid only for a pair plasma, because in an electron-ion plasma, the growth rate of the filamentation instability is diminished (Max et al. 1974;Sobacchi et al. 2020).

PIC SIMULATION SETUP
We carry out a suite of PIC simulations to investigate the non-linear interactions of radiation with pair plasma.
The primary code used in this paper is EPOCH (Arber et al. 2015, version 4.17.10; https://github.com/Warwick-Plasma/epoch), while some simulations are also carried out using the code Tristan-MP (Spitkovsky 2008, version 1; https://github.com/ntoles/tristan-mp-pitp), for the purposes of comparison and verification. We focus on EPOCH because of the numerical considerations discussed in §3.2, and for the efficiency gains resulting from the 1D simulation option it provides. Here, we outline the setup of the simulations and their numerical considerations, methods, and convergence, postponing the results to §4.

Initial conditions
We are interested in the time evolution of the pair plasma due to its interaction with a strong electromagnetic wave. Therefore, we initialize the simulations with a homogeneous, neutral, thermal plasma distribution, and introduce a planar, monochromatic wave, called the pumping wave. The wave has the form as illustrated in Fig. 1, filling the simulation box. Here, k 0 = 2π λ 0 and λ 0 are the angular wavenumber and wavelength. The amplitude E 0 = amω 0 e of the electric and magnetic fields is defined in accordance with Eq. (21), where ω 0 = k 0 is the angular frequency of the wave in vacuum.
We explore 2D setups for the study of both induced scattering and the filamentation instability. Taking advantage of the preferential backscattering of the former, we utilize 1D setups to explore induced scattering in more detail. Periodic boundary conditions are applied in all directions of the simulation box, as applicable in each case. The particles are initialized in a relativistic Maxwellian, accounting for oscillations in the field of the wave, as discussed in §3.3.

Electromagnetic wave initialization
Initializing a strong wave in a PIC code involves several subtle considerations. First, the electric and magnetic fields are located in different locations within the Yee lattice, and are advanced at different times. This combination greatly reduces the discretization errors, but require that the initial electric and magnetic fields be phase-shifted relative to each other, by a distance of the order of the grid length, in order to correctly capture the physical behavior of the electromagnetic wave. The initial magnetic field (on its staggered grid) thus requires an added phase (with respect to the initial electric field) of ∆x + 2∆t in EPOCH and ∆x + ∆t in Tristan-MP, where ∆x is the grid spacing and ∆t is the time step.
Second, there are errors associated with the interpolation of the fields to the particle locations, which is needed to calculate the electromagnetic force in the particle mover. In the Tristan-MP version used in the analysis, only linear interpolation is available, which entails a significant fractional error in the electromagnetic field, of order We explicitly alter the interpolation routine to correct the results specifically for the case of a monochromatic pumping wave. With this correction, we are able to demonstrate good agreement between EPOCH and Tristan-MP in this case (see Appendix §C). For multiple modes, such a tailored correction is no longer possible with linear interpolation, so we focus on EPOCH, which uses a triangular (second order) interpolation. Third, the discrete nature of the PIC code changes the speed of light. This effect exists for both EPOCH and Tristan-MP. Fourth-order interpolation in the field solver, employed in both codes, greatly reduces this source of error. For a sufficiently long wavelength relative to the grid spacing, the exact speed of light does not have a strong effect on the dynamics. One should however take the corrected speed of light into account, for example when using test particles.
Convergence tests, which indirectly also gauge these three problems, are discussed in §3.5.

Particle initialization and orbit drifts
As our simulations involve strong waves, initializing the particles in a thermal distribution with no regard to the field would lead to a violent energy conversion at the very beginning of the simulation. To reduce this effect, we place each particle in an orbit in the field of the vacuum electromagnetic wave (where the angular frequency is ω 0 , and a precise solution for the particle orbits is known), as follows.
The field of the wave is given in this case by the vector potential where is the phase. The particle orbit is described by the first integrals of motion (e.g., Melrose 1980), written in the notations of Lyubarsky (2019) as and where u x and u y are the respective components of the particle 4-velocity u = γv, and γ = (1 − v 2 ) −1 2 is the Lorentz factor. The particle velocities are expressed via the integrals of motion as The velocity of the particle guiding center (denoted by tilde) is found by averaging the velocities over the wave period,ṽ = ⟨v⟩. The integrals of motion are expressed via the velocities of the particle guiding center, whereγ = (1−ṽ 2 ) −1 2 is the Lorentz factor of the particle guiding center. To initialize our particles in approximate orbits, we first calculate the guiding center velocity,ṽ, and the corresponding Lorentz factor from a relativistic Maxwellian at the initial temperature T . Then we calculate the phase, η = −ω 0 x 0 , at the initial position x = x 0 and time t = 0. Finally, we use Eqs. (79)-(82) to calculate the initial values of v x and v y for the particle. This procedure results in much milder numerical transients.

Measuring growth rates in simulations
Sufficiently large simulations in 2D can resolve both induced scattering and the filamentation instability, simultaneously. However, to better study each process separately, we focus on two types of simulations: 1D simulations along the pumping wave, which can only resolve induced scattering, and 2D simulations confined to only one wavelength along the beam but many wavelengths in the perpendicular direction, which can only resolve the filamentation instability.
Both induced scattering and the filamentation instability result in perturbations in both the density and the electromagnetic field. We measure these perturbations in our simulations, in each diagnosed time step (snapshot), as follows. For the density field, we compute the root mean square (rms) of fractional density perturbations δN rms N 0 across the entire simulation box. For the electromagnetic perturbations, we perform a Fourier decomposition in traveling modes of wavenumber k, as described in Appendix B, and compute the normalized spectral distribution of the Poynting flux, S k S 0 , where S 0 ≡ S(t = 0; k = k 0 ) is the initial beam flux. We also measure the normalized total electromagnetic energy outside the beam, where t = 0 corresponds to the beginning of the simulation. More generally, when injecting multiple modes, δS S 0 is defined as the flux in all other modes, normalized to the initial, injected flux. The snapshots of a simulation are labelled by an integer τ = 0, 1, 2, . . . , τ max . The growth rate of a given quantity is measured by fitting an optimal range of snapshots during an exponential growth phase. The time interval between snapshots is denoted ∆T and chosen such that 1 ≲ ∆T λ 0 ≲ 2πT (ω 2 p λ 2 0 ma 2 ), in order to resolve the theoretical growth rate (of both induced scattering and filamentation), but not necessarily the pumping wave frequency. Our choice of ∆T λ 0 varies from 1 to 100, depending on the simulation parameters.
To outline the general procedure to measure growth rates from simulations, we denote the relevant perturbed quantity as Q; namely, Q = δS S 0 for 1D (induced scattering) simulations, and Q = δN rms N 0 for 2D (filamentation) simulations. To minimize noise and transients effects in the growth-rate measurement, we exclude earlytime snapshots where Q is small. To avoid saturation effects, we exclude late-time snapshots where Q is large. Quantitatively, our simulations indicate that the growth of Q may be susceptible to these spurious effects outside the interval 10 −3 < δS S 0 < 10 −2 in 1D, and the interval 10 −2 < δN rms N 0 < 2.5 × 10 −2 in 2D.
Consider two snapshots, τ 1 and τ 2 , confined within the relevant interval above. Modeling Q as exponential growth during all snapshots between τ 1 and τ 2 yields a best fit d(ln Q) dt = µ ± σ. The values of τ 1 and τ 2 are selected as the cleanest exponential behavior, identified with the minimal σ. Finally, depending on Q, we assign the measurement µ ± σ to either κ or Γ.

Simulation parameters and convergence
The nominal physical parameters in our simulations are chosen, as in Figs. 2 and 3, to represent a cold plasma with T m = 0.00137 and a monochromatic pumping wave with a short wavelength ω p k 0 = 0.107. The pumping wave amplitude a = 0.034 is chosen larger than in the above theoretical figures, in order to accelerate perturbation growth. The nominal numerical parameters are chosen differently for 1D and for 2D simulations, otherwise the latter would require too many resources.
For 1D simulations, we take N ppc = 120 particles per cell (both species) and λ 0 ∆x = 160 cells per wavelength, i.e., (ω p ∆x) −1 = 238 cells per skin depth. The simulation box length is chosen as L x λ 0 = 200 wavelengths, sufficient to spectrally resolve induced scattering. To see this, note that the wavelength shift for the fastest growing mode, as given by Eq. (28), is very small (≲ 10% for our nominal value of T m = 0.00137), so the box must be larger than giving ∼ 14λ 0 for nominal physical parameters. Our Tristan-MP simulations are effectively 1D, using only one non-guard cell in the y direction. The nominal numerical parameters for the 2D simulations are chosen as N ppc = 20 particles per cell and λ 0 ∆x = 80 cells per wavelength, i.e., (ω p ∆x) −1 = 119 cells per skin depth. The simulation box consists of square, ∆x = ∆y cells, and has dimensions L x λ 0 = 1, x being the incident wave direction, and L y λ 0 = 100. The extent along the y direction is chosen to accommodate the typical wavelength of filamentation in Eq. (60), giving λ f ∼ 20λ 0 for nominal physical parameters. We carry out a suite of convergence tests, and find good convergence (to 15% in growth rate) in each of our numerical parameters n ppc , λ 0 ∆x and L λ 0 . The agreement between EPOCH and Tristan-MP and the parametric convergence are demonstrated in Appendix C.
Simulation times are normalized to the angular frequency ω 0 of the pumping wave in vacuum. In addition, for easier comparison with theory, we also show times normalized to the theoretical maximal growth rates. For our nominal physical parameters, used in all evolution figures below, Eq. (22) yields κ max = 1.94 × 10 −3 ω 0 for induced scattering, and Eq. (55) yields Γ max = 6.98 × 10 −4 ω 0 for filamentation.
For brevity, nominal 1D simulations henceforth refer to 1D (induced scattering) EPOCH simulations with a monochromatic pumping wave and nominal physical and numerical parameters. Similarly, nominal 2D simulations refer to 2D (filamentation) EPOCH simulations with a monochromatic wave and nominal parameters.

RESULTS
We present the results from our 1D simulations, which resolve only induced scattering, in §4.1. Results from our 2D simulations, which resolve only the filamentation instability, are presented in §4.2. Overall, the results demonstrate a good agreement with the theory throughout the relevant parameter and dynamical ranges.

Induced scattering
Here we study the physical manifestations of induced scattering, from onset to saturation, and compare them to the theory. The case of a monochromatic wave is discussed in §4.1.1. The more realistic scenario, of an incident wave with a broad spectrum, is studied in §4.1.2.

Monochromatic beam
The simulated evolution of the induced scattering of a monochromatic strong wave is presented in Figs. 4-8. All figures are based on nominal 1D (induced scattering, nominal parameter, EPOCH) simulations, with two exceptions: the simulation of Fig. 4 uses a very high N ppc = 1200 (albeit a smaller L x = 100λ 0 ) to resolve the density perturbations, and Fig. 8 uses simulations where one physical parameter is varied, to present results for a wide range of physical parameters. showing the normalized evolution of (see legend) the incident beam energy (dot-dashed cyan), scattered electromagnetic flux (solid red with right axis), density perturbations (teal dotted; after removing small-scale noise by applying a high, k > 4k 0 filter; see text), and average plasma temperature (orange dashed). The nominal (physical and numerical parameters, EPOCH) 1D simulation is modified with a high Nppc = 1200 (but smaller Lx = 100λ 0 ) to resolve the density perturbations. The measured growth rate of δS S 0 agrees with its theoretical fastest growth rate (thin magenta dot-dashed slope; obtained by solving Eq. (22) numerically) and our fitting procedure (thin black dashed slope).
The temporal evolution of the different plasma components is shown in Fig. 4. After an initial, oscillatory transient, the scattered electromagnetic modes are seen to grow exponentially, at a rate consistent with theory (nearly parallel measurement in solid red curve, theory in thin dot-dashed magenta, and automatic fitting result in thin dashed black). As anticipated above and shown numerically below, this growth is strongly dominated by the fastest growing mode, so a good agreement is obtained even when comparing the growth rate of all modes combined, δS S 0 , to the theoretical growth rate κ max of the fastest growing mode. The generation of scattered modes is accompanied by the growth of density perturbations (dotted teal). To better resolve the growth in δN rms N 0 , we lower the respective noise level from ∼ 20% to ∼ 0.7% by removing small-scale noise using a high, k > 4k 0 filter, as motivated by Fig. 5 below. The exponential growth of induced scattering modes is seen to saturate around t ≃ 8000ω −1 0 ≃ 15κ −1 max . Subsequently, the electromagnetic modes still grow, but no longer exponentially, and the built-up density perturbations decline nearly exponen-tially. The decay of the incident wave (dot-dashed cyan) and the heating of the plasma (dashed orange) accelerate after saturation, and are approximately exponential with a rate ∼ 3 × 10 −5 ω 0 .
max ; dot-dashed black) and late (t ≃ 7700ω −1 0 ≃ 15κ −1 max ; solid blue) exponential growth, and in saturation phase (t ≃ 15000ω −1 0 ≃ 30κ −1 max ; dashed red). Beating at an angular wavenumber k d is seen in the density field (inset; smoothed with a Gaussian of standard deviation λ 0 ) and as the main peak in the power spectrum (of data folded over 20λ 0 to reduce noise), consistent with the expected (magenta triangle and vertical line) k d k 0 ≃ 2 − 2(T m) 1 2 of Eq. (28). A shift to slightly lower k and the growth of the first harmonic are seen to develop at the later times. Fluctuations in electric charge (Np − Ne) Ntot (thick cyan, for the later time) remain small and consistent with small-scale white noise.
The density perturbations can be seen at earlier times if one takes a power spectrum of the density field, as demonstrated in Fig. 5 for three characteristic stages: early linear phase (κ max t = 7), late linear phase (κ max t = 15) and saturation phase (κ max t = 30). Even at early times, the figure shows pronounced beating at the expected angular wavenumber (magenta triangle), k d k 0 ≃ 2 − 2(T m) 1 2 ≃ 1.93 (for nominal temperature); see Eq. (28). At later times, a second harmonic emerges, and a shift of both primary and second harmonic to slightly lower k is seen to develop, consistent with plasma heating. At small scales with k ≳ 4k 0 we find only numerical noise even at late times, justifying the high k filter used in Fig. 4. We verify in Fig. 5 (cyan thick curves) that fluctuations in electric charge (N p − N e ) N tot remain small and are consistent with small-scale white noise.
To illustrate the scattered spectrum in more detail, Fig. 6 presents the electromagnetic spectrum derived from a Fourier decomposition in traveling modes (see Appendix B). The spectrum is shown for different times on the left panel, with an evolution plot (similar to Fig. 4) for reference in the right panel. At early times, the back-scattered beam is seen to grow exponentially, sharply peaked at k ≃ −0.9k 0 wavelengths: consistent with the fastest growing mode being reversed and with a diminished, ∆ν ν ≃ 2(T m) 1 2 frequency. After sufficient energy is transferred from the incident beam to the backscattered wave, the exponential growth saturates, and the reflected beam becomes broader in k-space. The reflected beam itself is then back-scattered, leading to the growth of a broad-k beam in the original direction, peaked around k ≃ 0.8k 0 . At early times, there is a good agreement between the energy in the back-scattered mode and in all k ≠ k 0 modes combined, but this is no longer true after saturation. The figure demonstrates the basic picture of the induced scattering process (e.g., Zel'Dovich et al. 1972): radiation energy is distributed to lower frequencies while conserving the total number of photons, so the process necessarily involves some plasma heating.
To quantitatively compare the simulated mode growth with the induced scattering theory, Fig. 7 shows the growth rate spectrum both as measured in the exponential phase of our nominal 1D simulation, and as derived asymptotically in Eq. (26). For simplicity, we estimate the growth rate here as max was chosen manually, and not by our nominal best-fit procedure. Although no smoothing is used in this figure, one can appreciate the good agreement between the simulated result (solid red curve) and the theoretical curve (dashed green).
Finally, we present in Fig. 8 a suite of simulations, showing how the induced scattering rate changes as a function of the plasma temperature (left panel) and the strength of the wave (right panel). The measured rates from the otherwise nominal 1D simulations (error bars inside symbols) are compared to theory, showing both the numerical solution to the dispersion equation (22), and the asymptotic relation Eq. (29). Qualitatively, we find a good agreement between simulation and theory, within their respective limits, although we are comparing slightly different quantities: all k ≠ k 0 modes outside the incident wave in the simulations, with the fastest growing mode in theory. The red error bars inside circles show the fitting uncertainties of our rate measurements. The black error bars inside diamonds show that even when combining fitting and convergence (from Richardson extrapolation) uncertainties, the resulting error estimate remains smaller than the approximation accuracy.

Broad-spectrum beam
Next, consider the general case of a nonmonochromatic incident beam, composed of many modes with random phases spanning a broad spectrum, as typically found in astrophysical applications. We simulate the evolution of different such beams, defined as in §2.1 by a superposition of monochromatic waves, analogous to the distribution a ω used e.g., in Eq. (35). To simplify the comparison to the monochromatic case, we normalize the fields such that the total energy of any given beam equals that of the nominal monochromatic wave.
Consider a wave packet with a Gaussian energy distribution, centered on wavenumber k 0 with a standard deviation σ. The implied injected fields are given by where φ k is a random phase, and the summation is carried out over all k modes up to the k − k 0 ≃ 7σ level around k 0 . The factor 1 2 introduced to the standard Gaussian exponent renders σ the standard deviation in the energy, rather than in the field.  Fig. 4). The spectrum of electromagnetic flux normalized to the incident flux, S(k) S 0 , is plotted (curves in left panel, with coincidental vertical lines in the right panel) at times κmaxt ≃ 10 (solid blue), 20 (dashed magenta), 35 (dash-dotted green), and 55 (dotted black). The spectrum has been smoothed using a moving average of ten k-modes, except in the shaded region, which contains the pumping, k = k 0 wave. The extended spectrum of the beam should have a significant effect on the evolution of the plasma if it is not masked by the thermal spread, i.e., when (σ k 0 ) 2 ≫ T m. Numerically, one can represent the beam only by a finite, discrete set of modes, due to the minimal, δk k 0 = λ 0 L x spacing in k space allowed by the periodic simulation box. To avoid numerical artefacts due to this non-continuum spectrum, we require that the discrete-ness is masked by thermal broadening, (δk k 0 ) 2 ≪ T m. Combined, we thus require all simulated broad beams to satisfy where N ≃ 2σ δk is the number of modes injected within 1σ of k 0 in the numerical representation of the Gaussian in Eq. (86). We perform a suite of nominal 1D simulations, but replacing a with a spectrum a ω carrying the same total energy as in the monochromatic case. This spectrum is chosen as the Gaussian beam (86), with a range of widths: σ k 0 ≃ 0. 04,0.07,0.1,0.13,0.16,0.18,and 0.21. Conditions (87) are then satisfied, except possibly for the first few, narrow Gaussians, where (σ k 0 ) 2 (T m) ≃ 1.2, 3.6, . . ., so the spectral extent of the beam may be partly washed out by the thermal broadening.
All these simulated broad beams yield growth rates that qualitatively agree with the theory and are modified with respect to the monochromatic case. The case σ = 0.1k 0 is demonstrated in Figure 9, which shows the growth rate of different reflected k-modes. Despite the noise obtained in the absence of any smoothing, the results are seen to be in good agreement with the asymptotic growth rate profile. This profile is computed numerically from Eq. (35), which convolves the monochromatic rate Eq. (26) in the slow-growth limit (23) with the spectrum of the beam.
To analyze the temporal evolution of the broad-beam simulations, we use the generalized definition of δS S 0 as the electromagnetic flux outside the injected, k −k 0 < 7σ beam, normalized to the initial injected flux. The timedependence of this quantity is shown in Figure 10 (left panel) for four representative values of σ k 0 . Extended exponential growth phases are resolved, so the growth  30)). The black error bars inside diamonds demonstrate the combination of both fitting and convergence uncertainties, seen to be quite small. rates can be accurately measured. Note the irregular changes in saturation phase behavior as a function of σ.
The measured exponential growth rates of this generalized δS S 0 are shown in the right panel of Figure 10, for all simulations, including both broad beams and the monochromatic, σ k 0 = 0 case. These simulated rates are seen to be in qualitative agreement with the maximal theoretical growth rates inferred from the asymptotic Eq. (35). As expected, the rate declines with in-creasing σ.

Filamentation instability
We now turn our attention to the filamentation instability. A strong electromagnetic wave is injected into the pair plasma as in the induced scattering case, but here it is necessary to resolve long-wavelength perturbations in the perpendicular, y-direction. Hence, we use a 2D box with a large, L y = 100λ 0 transverse extent. Computation resource limitations then dictate using lower values of N ppc and λ 0 ∆x (see §3.5), and the minimal L x = λ 0 possible still able to sustain the incident wave.
Normalized perturbations in the density profile N (y), after averaging along x, are shown in Fig. 11 (inset) for the nominal 2D simulation. An oscillatory mode of wavelength ∼ 100λ 0 6 ≃ 17λ 0 is seen to grow exponentially, evolving into dense, merging filaments separated by low-density plasma as saturation is approached. Indeed, a power-spectrum of N (y) (figure body) shows a fairly sharp peak at k ≃ 0.06k 0 early on, progressively broadening and shifting to lower k at late times. This wavenumber is broadly consistent with the k d ≃ (ω p a 0 2)(T m) −1 2 ≃ 0.05k 0 estimate of the fastest growing mode in Eq. (60), for our nominal physical parameters. Note that Eq. (60) should not be accurate here as it depends on the assumption ω d ≪ k d v T , or equivalently Eq. (61), which is only weakly satisfied for these parameters; solving Eq. (55) numerically gives an even lower, k ≃ 0.04k 0 . The figure also shows that even at late times, there is only noise at small, k ≳ k 0 3 scales, and that the plasma remains neutral at all times.
The temporal evolution of the RMS density fluctuations in a nominal 2D simulation is shown in Fig. 12. In order to lower the noise level and resolve a more prolonged period of exponential growth, here we first apply a high-k filter, removing all k > k 0 3 modes which are indicated as noise by Fig. 11. The exponential growth rate is seen to be in good agreement with the fastest growth rate found by solving Eq. (55) numerically (nearly paral- Right: The estimated growth rates (of δS S 0 ; red circles with fitting error bars) and asymptotic maximal rates (green dashed; obtained from Eq. (35)), as a function of σ. Fig. 11.-Density perturbations (thin curves) from a nominal 2D (filamentation) simulation, shown during early (t ≃ 5.7 × 10 3 ω −1 0 ≃ 4Γ −1 max ; dot-dashed black) and late (t ≃ 8.6 × 10 3 ω −1 0 ≃ 6Γ −1 max ; solid blue) exponential growth, and during saturation (t ≃ 1.1×10 4 ω −1 0 ≃ 8Γ −1 max ; dashed red). Filaments are seen forming in the density field (inset) and as the pronounced peak in the power spectrum (unfolded; figure body), consistent with the fastest growing mode given by Eq. (60) (magenta triangle and vertical line). Filament merger and spectral broadening are seen to develop at late times. Fluctuations in electric charge (Np − Ne) Ntot (thick cyan, for the later time) remain small and consistent with small-scale white noise.
lel measurement in solid red curve, theory in dot-dashed magenta, and automatic fitting result in dashed black).
Finally, we carry out a suite of simulations to investigate the filamentation instability as a function of physical parameters. Figure 13 shows the dependence of the estimated growth rate Γ of the instability upon T m (left panel) and a (right panel). These rates are compared to the fastest growth rate derived from the numerical so- Fig. 12.-Evolution of density perturbations along y in a nominal 2D (filamentation) simulation. The measured (δNrms N 0 ) is shown (solid red curve) after averaging over x and removing shortwavelength, k > k 0 3 modes to lower the noise level. The perturbations undergo a period of exponential growth, in agreement with the theoretical fastest growth rate (magenta dot-dashed slope; obtained by solving Eq. (55) numerically); our automatic fitting measurement is also shown (black dashed slope).
lution to the dispersion relation Eq. (55), as well as the simulations. The measured rates (red fitting error bars inside circles) are based on the RMS density perturbations, δNrms N 0 , after averaging along x and filtering out all k > k 0 3 modes. The results are compared to the theoretical maximal growth rates, calculated by solving Eq. (55) numerically (solid blue curve) or using the asymptotic estimates (dashed green; minimum of the two rates in Eqs. (59) and (64)). The black error bars inside diamonds, demonstrating the combination of both fitting and convergence uncertainties, are seen to be quite small. approximate asymptotic relations given by Eqs. (58) and (63). One sees that the growth rates obtained from the simulations agree very well with the theory, over a wide range of physical parameters.

SUMMARY
This is the first paper in a series devoted to the study of nonlinear interactions of powerful electromagnetic radiation with pair plasma. These effects could play a role in different astrophysical environments, including AGN, pulsars and interstellar masers, and drew considerable interest in recent years, after the discovery of fast radio bursts. In this paper, we study the case where the strength parameter of the wave, defined in Eq. (21), is small, a ≪ 1. Here, particles oscillate in the field of the wave with velocities well below the speed of light. This case can thus be studied analytically, facilitating the development of numerical tools by comparing the theoretical and numerical results.
At the initial stage, perturbations caused by the strong radiation beam grow exponentially, and the evolution may be considered as an instability of the beam. Therefore, the first question is what are the characteristics and growth rates of the relevant instabilities. The two main effects for a pair plasma are induced scattering and the filamentation instability. We derived dispersion relations for these instabilities in Eqs. (22) and (55), respectively, from which the corresponding growth rates are easily calculated numerically. The well-known asymptotic relations in different limiting cases are easily reproduced from these dispersion relations.
We performed extensive numerical simulations to study the temporal evolution of an initially homogeneous pair plasma exposed to powerful radiation, and measured the properties of the resulting instabilities. The simulations were carried out using both EPOCH and Tristan-MP, with good agreement between the two codes, and with small statistical and convergence uncertainties (see Appendix C and Figs. 14-15). By choosing long, narrow simulation boxes either parallel (L x ≫ L y , and nominally with L y = 0) or perpendicular (L y ≫ L x ) to the incident beam, we select for either induced scattering or filamentation growth, allowing us to study each effect separately. The results, assisted by an algorithm for identifying the exponential growth period ( §3.4) and a decomposition of fields in traveling waves (Appendix B), are shown in Figs. 4-11 for induced scattering, and in Figs. 12-13 for the filamentation instability.
The simulations show the initial exponential growth of perturbations in electromagnetic fields and in density, and the transfer of energy from the incident beam to lower frequencies and to particle oscillations, in good quantitative agreement with theory for a wide range of physical parameters (see in particular Figs. 8 and 13) and across multiple diagnostics. The simulations trace the evolution deep into the saturation phase, showing non-linear structure growth, the dissipation of the density perturbations, and substantial plasma heating. For the case of induced scattering, where the spectrum of the beam is important, we examine both monochromatic and broad beams. The growth becomes slower as the incident spectrum is broadened, in agreement with theory, implying that for astrophysically relevant beams, filamentation occurs well before the induced scattering comes into the play.
Having established consistent theory and simulations in the a ≪ 1 regime, with different codes and multiple diagnostics, future papers in this series will explore the strong a regime, the addition of background magnetic fields, the inclusion of ions, etc.

A. NONLINEAR CURRENT: MONOCHROMATIC WAVE
Here we briefly outline the derivation of the nonlinear current in a plane monochromatic wave (Montgomery & Tidman 1964;Sluijter & Montgomery 1965). The current is produced by electrons oscillating in the field of the wave, see Eq.
To first approximation, this yields the linear current (5). The next approximation is of the third order in a. Two effects contribute to the current in this order: the density variation in the field of the wave, and the nonlinear corrections to the particle oscillation velocity, j nl 2 = N 0 qv (3) .
Therefore, one can conveniently present the nonlinear current as a superposition of two terms. The first contribution arises from the Lorentz force (the second term in the RHS of Eq. (3)). Due to oscillations with velocity (A2), the magnetic field of the wave produces a longitudinal (along the direction of the pumping wave) force at the doubled frequency, qv × ( ( (∇ × A) ) ) = − 1 2 mk 0 a 2x sin 2η.
This force yields longitudinal oscillations in the second order in a: v (2) = 1 4 k 0 a 2 ω 0x cos 2η.
These oscillations produce, by virtue of the continuity equation, the density perturbation, δN (2) = (k 0 ω 0 )N 0 v (2) x . The beating between the density oscillations at the double frequency and the velocity oscillations (A2) contributes to the nonlinear current (A3). Taking into account that δN is a quadratic function of the particle charge, so that the electron and positron fluctuations are equal, one gets 4πj nl 1 = − ma 3 k 2 0 8eω 2 0 ω 2 p (cos 3η + cos η)ŷ.
The second contribution to the nonlinear current arises because relativistic corrections effectively raise the mass of an oscillating particle. In the plane wave, the component of the generalized momentum perpendicular to the propagation direction is conserved, p y +qA y = const.; one can check directly that this relation satisfies the equation of motion (3) exactly. In our case, this relation reduces to v y √ 1 − v 2 + q e a cos η = 0.
In the first approximation, this relation yields Eq. (A2). In the next approximation, one finds that v y has a correction third order in a: This yields the second contribution to the non-linear current, Eq. (A4): 4πj nl 2 = ma 3 8e ω 2 p (cos 3η + 3 cos η)ŷ.
Together, j nl 1 and j nl 2 give the nonlinear current. However, only resonant terms (those oscillating with the frequency ω 0 ) should be taken into account in the wave equation because only perturbations from these terms accumulate at the large scale whereas nonresonant terms produce small perturbations that do not grow. Retaining only resonant terms yields which gives Eq. (43) in the high frequency limit, k 0 = ω 0 .

B. EXPANDING THE FIELDS IN TRAVELLING WAVES
Simulation snapshots provide the discretized fields E(r) and B(r) at each time step. One could expand the fields in Fourier series, E(r) = k e k e ik⋅r ; where e k = e * −k and b k = b * −k because E and B are real. However, such an expansion does not facilitate the separation of waves propagating in opposite directions. In order to analyze the results of simulations of the scattering process, we expand the fields in travelling waves, E(r, t) = k E k e i(k⋅r−ωt) + c.c.; B(r, t) = k B k e i(k⋅r−ωt) + c.c..

(B2)
In this case, the direction of k gives the direction of the wave. The Poynting flux is presented as so that the flux of each wave is By comparing Eqs. (B1) and (B2), we find Here we take, without loss of generality, t = 0. These two equations are supplemented by two Maxwell equations; for the travelling waves (B2) they are written as Now we have four equations, so we may express the amplitudes of the travelling waves via the Fourier components of the snapshots: wherek is the unit vector in the direction of the wave, and we take into account that ω = k for high-frequency waves. Substituting these relations in Eq. (B4), we get the Poynting flux in each wave as The numerical parameters associated with our PIC simulations are N ppc , λ 0 ∆x, and either L x λ 0 (for 1D simulations) or L y λ 0 (for 2D). Convergence tests based on nominal 1D simulations demonstrate the convergence of the growth rates inferred from our automatic fitting procedure (see §3.4), with respect to all numerical parameters, using a large suite of simulations. For example, one test examines the convergence of the results with a control parameter C, which scales both N ppc and λ 0 ∆x simultaneously, thus multiplying the number of particles per wavelength by C 2 .
We test the agreement between EPOCH and Tristan-MP results using a suite of 1D simulations, based on the nominal parameters but modified by different values of C, where C = 1 is the nominal case. Figure 14 shows the evolution of the fractional electromagnetic perturbations for four different choices of C, in both EPOCH (lines) and Tristan-MP (markers). The growth rates inferred from the different simulations are consistent with each other and with the theoretical value (magenta dot-dashed slope). Figure 15 demonstrates the convergence of the fitted growth rates obtained as a function of the scaling factor C in our otherwise nominal 1D (left panel) and 2D (right panel) simulations. Richardson extrapolation to C → 0 suggests convergence uncertainties of order ≲ 15%.