Introduction-Modeling cosmic-ray (CR) transport in different environments requires a fundamental understanding of the interplay of various physical processes. Based on measurements of secondary-to-primary CR flux ratios [1][2][3], CRs seem to diffuse due to scattering off turbulent magnetic fields. Whether these fields are self-generated or are injected by supernova explosions on large scales and cascade from large to small scales is unclear. The discovery of a break in secondary-to-primary ratios [1][2][3] might hint at CRs generating these fields via the excitation of streaming instabilities [4][5][6], while preexisting turbulence may take over at energies ≳TeV.

The problem of CR transport in the presence of streaming instabilities excited by the same CRs is intrinsically nonlinear as the growth rate depends on the CR distribution, which in turn is determined by the way particles propagate. Especially within the context of Galactic transport on kpc scales, the resonant streaming instability seems a promising candidate to describe CR confinement, at least up to ≲TeV energies [4,5] and perhaps even at all energies [7].

A crucial role in the assessment of the role of selfgeneration for CR transport is played by damping processes, which are thought to be responsible for saturating the growth of the instability and eventually determining the energy dependence of the diffusion process.

If growth occurs fast enough compared with damping, then the CR population has time to isotropize in the reference frame of the waves. Hence, their bulk speed v D , initially super-Alfvenic to excite the instability, is reduced to the Alfvén speed v A . While this may be important for low energy CRs, around GeV energies, the decreasing energy dependence of secondary-to-primary ratios suggests that higher energy particles move with a super-Alfvénic drift speed, v D > v A , that increases with energy. This is in fact expected if the growth of the resonant streaming instability is saturated due to damping processes [8]. The common expectation is that waves of wave number k grow due to the super-Alfvénic drift but that the level of perturbations is determined by the equality of the growth rate and the damping at the same wave number. Depending on the phase of the interstellar medium, the most relevant mechanisms are typically ion-neutral damping [8] and nonlinear Landau damping (NLLD) [9], although many other processes might be at work [10][11][12][13]. In the hot ionized medium, the absence of neutrals renders NLLD the dominant process to limit the wave growth both on galactic scales [9] and even near sources [14][15][16].

The process can be pictured as follows: two circularly polarized Alfvén waves of slightly different wave number combine into a beat wave that can resonantly interact with particles of the background plasma, transferring energy from the waves to the plasma in the form of heat [9]. For two waves of opposite polarizations, both waves are damped, and the background plasma is heated. For the same polarization, the shorter wavelength wave transfers energy to both the other wave and the background plasma.

In the case of a spectrum of waves, every wave can combine with all the others. Thus, a given wave should be damped or grow due to interactions with all other waves, an effect often overlooked in the literature [15][16][17][18][19]. If there is an equal amount of both polarizations (linearly polarized spectrum), the growth due to the interaction with a shorterwavelength wave of the same polarization will be canceled by the interaction with the same wave of opposite polarization and the two waves of the same wavelength will damp each other. Naively, this configuration should apply to the resonant streaming instability of an almost isotropic CR distribution generating forward-moving left-and righthanded circularly polarized waves in nearly equal amounts [20]. On the other hand, if there is an excess of a particular polarization, energy could be transferred to larger wavelength waves, which might result in an inverse cascade.

As discussed above, having in mind CR transport in our Galaxy, the streaming instability is assumed to be balanced by NLLD. In particular, the phenomenon of nonlinear transport was studied on Galactic scales [4][5][6], around SNRs [15,16,21], in the multiphase ISM during Galaxy formation [22], in Galaxy clusters [19], in heating the Galactic gas [18] and in launching Galactic winds [23][24][25]. On the other hand, numerical attempts at simulating the saturation of the instability focused mainly on the β ≪ 1 regime in which NLLD is negligible or on ion-neutral damping, employing particle-in-cell (PIC) [26] or magnetohydrodynamics ðMHDÞ þ PIC [27,28] methods.

In this Letter, we focus on the saturation of the resonant streaming instability in the astrophysically relevant regime of β ≳ 1 utilizing hybrid PIC simulations. We illustrate the first conclusive evidence that in such plasma NLLD is active and that it leads to a reduction of the drift speed of CRs driving streaming instability. Furthermore, contrary to previous studies we keep driving the instability instead of having a finite energy available in CRs with periodic boundaries. This situation is reminiscent of a halo plasma element adjacent to the Galactic disk, that is constantly bombarded with freshly accelerated particles entering the plasma element from the disc region. Moreover, we show that, contrary to common approximations, modes on a given scale λ are damped by all the power contained on larger scales. Finally, we find that there is a small anisotropy of left-to right-handed modes created by the streaming instability which eventually drives an inverse cascade, as mentioned above.

Next we describe the code and our simulation setup, then we present our results and conclusions.

PIC simulations-We study the self-confinement of CRs resulting from excitation of a resonant streaming instability and NLLD limiting the growth, using dHybridR, a relativistic hybrid code with kinetic ions and (massless, chargeneutralizing) fluid electrons [29,30]. A detailed discussion of the hybrid limit and of our code can be found in Sec. 2 of [30]. We consider an adiabatic closure for the electron pressure, i.e., P e ∝ ρ 5=3 [31]. dHybridR is suited to properly simulate CR streaming instabilities [30,[32][33][34] and should capture all the relevant mechanisms for saturation of the resonant streaming instability, including NLLD (see also [35]). By choosing a hybrid code, we are able to simulate longer timescales compared to PIC codes, essential for exploring low growth rates of the resonant streaming instability.

In our simulations, velocities, timescales, and length scales are normalized to the Alfvén speed v A , the ion cyclotron time Ω 0 , and the ion skin depth d i ¼ v A =Ω 0 , respectively. Magnetic fields are normalized to the background magnetic field B 0 . Particle number densities are normalized to the initial density of the background plasma, n 0 .

Our simulation box is quasi-1D along x and retains all three components of the momenta and electromagnetic fields. The total size is 60 000 d i × 5d i , divided into 120 000 × 10 cells. The speed of light is fixed to c ¼ 100v A and the background magnetic field is directed along the x axis with strength B 0 . The background plasma is periodic in all directions and is sampled from a Maxwellian with an ion plasma β ¼ 2v 2 th;i =v 2 A ¼ 2 and 128 particles per cell, in order to reduce the magnetic field noise floor. The electromagnetic fields are allowed to leave the simulation box, having open boundaries along the x direction in order to prevent the recycling of the same waves. CRs have periodic conditions in y, open boundaries in x, and we use 16 particles per cell.

The injected CR distribution is isotropic in a CR rest frame with p iso ¼ 100mv A (Lorentz factor ffiffi ffi 2 p ) and boosted into the simulation frame with p bst ¼ 5mv A leading to a drift speed in the background plasma of 4v A , i.e., mildly super-Alfvénic. The time step is 0.005Ω -1 ci . Though some of the chosen parameters are not strictly equal to the values in the Galaxy (e.g., c=v A ), our setup respects the hierarchy of physical scales, which allows us to properly study the problem.

Our simulation box is 600 CR gyroradii large, r L ¼ 100d i , to accurately resolve the maximally growing wave number k res ¼ r -1 L and to capture the scattering mean free path for magnetic fields as low as δB 2 =B 2 0 ≈ 2 × 10 -3 . We set n CR ¼ 1.8 × 10 -4 n 0 , so that the growth rate expected for our distribution is

An important detail in our setup is the implementation of the injector for the CR distribution. Because of CRs having v iso ≫ v D , there are CRs with positive and negative v x in our box. Hence, we have to inject particles with v x > 0 at the left and v x < 0 at the right boundary. When particles start to scatter, the amount of particles moving to the left and right changes relatively, leading to a decrease in v D . Hence, one could either keep injecting particles with the initial distribution or readjust the injector to inject the new distribution inside the box. For simplicity, the runs presented here adopt the former prescription, though we tested PHYSICAL REVIEW LETTERS 134, 045201 (2025) also the latter and concluded that our boundary conditions do not qualitatively affect our main conclusions.

Results-Evolution of the magnetic field: The overall evolution of the CRs and the magnetic field is illustrated in Fig. 1. We can identify four phases, that indicate different interplay of CRs and fields. In stage I, CRs stream freely through the box with the initialized drift speed. The magnetic field's growth rate agrees well, within 10%, with the prediction derived from [20]. Around t ¼ 3γ -1 max , the drift speed starts to drop as a result of CRs scattering off the self-generated magnetic perturbations (Stage II). The transition occurs when the mean free path is comparable to our box size, i.e., δB 2 =B 2 0 ≈ 10 -3 , see Eq. (1). At this point, scattering influences the evolution of B 2 ⊥ by reducing the growth rate from the unperturbed one of our initial, pristine CR distribution to a new one given by the self-scattered distribution.

After ∼17γ -1 max , stage III starts in which the field growth slows down significantly. In the absence of damping, the fields would keep growing and scattering would increase, until CRs become isotropic in the reference frame of the waves and get advected with v D ¼ v A . Hence, the super-Alfvénic drift throughout this stage with v D ≳ 3v A , indicates that other processes are opposing the growth of the waves. This point is even more evident when looking at the top panel of Fig. 1.

In stage III, the magnetic field contained on all scales that are able to resonantly scatter CRs plateaus and then slightly damps away, despite v D > v A . After t ≈ 60γ -1 max , the total field and the drift speed reach a plateau indicating the final phase of saturation (stage IV).

Nonlinear Landau damping: To understand the increasing drift speed observed in stage III, it is instructive to first determine the reason why the power on resonant scales saturates. Figure 2 shows the background ion distribution function at three different times. At momenta around p x ¼ mv A , the whole distribution flattens due to particles in Landau resonance with the beat wave [36] which is a clear signature of NLLD. The feature appears only in the p x component and not the perpendicular ones, which supports our interpretation. As illustrated, the feature is already present around t ≈ 7γ -1 max , meaning that damping started to transfer energy from waves to the plasma before that, but remains too weak to compensate the growth due to CR streaming during stage II. However, as the field continues to grow, damping can eventually saturate the field at different scales, setting the beginning of stage III.

Naively, we would expect fields to grow only on scales that can be resonant with the driving particles. This expectation holds true initially, as can be seen from the power spectra of right-and left-handed waves in Fig. 3 and the large-scale fields being flat in Fig. 1. However, once NLLD starts to operate, it is evident that the field starts growing on scales that are not in resonance with any of the particles. Furthermore, contrary to common approximations [15,18,37,38], we confirm that small scales are damped by all the power contained on larger scales. This renders the damping rate maximal at small scales and gradually decreasing with k, compatible with earlier predictions [39]; more precisely, the power at small scales saturates first and starts partially damping away during stage III (see Fig. 3). One subtle point is that the magnetic field on large scales keeps growing throughout stage III. Hence, it is far from trivial that the damping rate exactly compensates the growth, and the resonant field is not entirely damped away. Since the growth rate in our simulation changes due to scattering, it is impossible to determine the exact functional form of the damping rate and whether nonlinear effects, such as particle trapping [39], modify the predictions of linear analyses. Despite these issues, that deserve further FIG. 2. Background ion distribution function as a function of p x at three different times, showing a flattening at p x ¼ mv A as a signature of NLLD.

investigation, there is enough evidence that NLLD at work in our simulation and that it plays an important role to limit the growth of the waves. Overall, it is clear that magnetic power is transferred from small scales to large scales; a similar effect was observed in PIC simulations by [40].

Interestingly, such an inverse cascade is a feature of NLLD that was discussed in the early literature on this topic [9], although it only occurs when two circularly polarized waves of the same polarization interact; to our knowledge, it has never been shown in hybrid simulations of the resonant streaming instability. Note that, in principle, the resonant streaming instability should excite linearly polarized waves, so that this effect should be absent. However, there is a small anisotropy of left-to right-handed modes of the order of v D =c, since there are more particles with positive pitch angle than with negative one. In Fig. 3, one can see that initially the two helicities grow almost identical, i.e., the field is linearly polarized. However, the small anisotropy is sufficient to halt the growth of the right-handed modes and render the final spectrum partially circularly polarized. This kickstarts an inverse cascade that becomes efficient enough to remove power from small scales and relocate it to scales larger than the CR gyroradius. As a result, opposite to what is always assumed, the energy given by CRs to the fields is not entirely dumped into heating the background gas, but rather it is partially transferred to large-scale modes for which damping is inefficient. In terms of CR scattering, this opens the possibility for low energy CRs to influence the transport properties of higher-energy particles, with potentially profound implications for our description of CR transport. We defer to a future paper the study of this effect.

Drift velocity: The removal of power from small scales is reflected in the drift speed because it effectively worsens the so-called 90°-scattering problem, meaning that particles cannot cross the pitch angle barrier at μ ¼ 0 and do not manage to isotropize. Overall, the increase in v D can then be understood as the pitch angle barrier becoming increasingly difficult to cross due to damping removing power from small scales. We can estimate the mean free path of a particle with a given μ ¼ ð1=kr L Þ by calculating the product of its velocity along x and the quasilinear theory scattering time resulting in

The particles efficiently scatter if λ is smaller than the box size L which means we need a magnetic field at the scale resonant with a given μ of ∼10 -3 . From Fig. 3, it is clear that during stage II the field at small scales grows sufficiently large to ensure pitch angle scattering over a large range of values of μ. However, at late times it gets damped away below this threshold, which is why the drift speed starts to increase, as the pitch angle barrier becomes effectively too large to cross.

After the removal of power on small scales, the increase in drift speed slows down and approaches a new equilibrium value v D ∼ 3.7v A , given by diffusion in the almost constant magnetic fields at the resonant scales. However, the power in large-scale fields is unexpectedly still increasing during that stage and even surpasses the power on resonant scales around t ∼ 50γ -1 max . Since waves are allowed to leave the box, we expect that an Alfvén wave leaves the box at most after a time L=v A ∼ 26γ -1 max . On the other hand, the energy transfer of CRs to magnetic fields and subsequent redistribution as heat and large-scale fields is a highly nonlinear phenomenon depending on the fields at resonant scales and the CR drift speed among other parameters. As such, saturation in our box can only be achieved with a delay with respect to all these other processes, a delay that can be estimated as roughly one crossing time. Thus, we expect the green line to flatten with a delay of 26γ -1 max with respect to the drift speed and the peak of the orange line which would be around t ∼ 56γ -1 max . Indeed, the growth halts around t ∼ 60γ -1 max initiating stage IV and a steady state is reached throughout our box. At this point, the inverse cascade ends due to the waves leaving the box. In a closed periodic box, where waves cannot leave the box and are recycled on the other side, saturation is impossible to achieve, as we will discuss in detail in future work.

The configuration in stage IV is qualitatively what we would expect to happen in the Galaxy: A CR distribution that is shaped by diffusion due to self-generated turbulent magnetic fields.

A crucial new ingredient is the appearance of the largescale fields, a by-product of self-generation. At present, it is not clear what could be the fate of such large scale modes: in principle, these waves could be absorbed resonantly by Power spectra of the magnetic field δB y ∓ iδB z at different k for left-handed (top panel) and right-handed (bottom panel) modes, respectively. Small-scale power is transferred to large scales, indicating an inverse cascade.

higher energy particles, that in this way may limit their growth. The spread in momentum in our simulation is too small to address this issue. It is however a point of major since in models of Galactic CR transport, largescale turbulence is typically assumed to be injected only by supernova explosions [4,5] or produced by high-energy CRs in resonance with these scales [6].

Instead, our simulations indicate that large-scale fields naturally arise from NLLD even in the absence of these effects and can grow magnetic fields larger than the ones on resonant scales. In principle, a new intriguing picture is outlined: low-energy CRs can grow fields responsible for scattering high-energy CRs. The phenomenological consequences of these findings will be explored in a forthcoming work.

Conclusions-We discuss the first confirmation of the role of NLLD as the main process regulating the growth of the resonant streaming instability in a high β plasma. This phenomenon has been discussed for several decades now as possibly responsible for the diffusion of CRs on Galactic scales, and the self-confinement around sources.

In the absence of damping processes, it is expected that the growth of the resonant modes continues until full isotropization of the particles is reached in the wave frame, a phenomenon resulting in advection of CRs with the waves. On the other hand, in the Galaxy the diffusion of CRs with E ≳ 10 GeV corresponds to super-Alfvénic drifts, which can be accommodated if the growth of the resonant instability is limited by damping.

In the Galactic halo, the dominant damping process is most likely NLLD, and previous works have shown that a self-generated halo can form under these assumptions [4][5][6].

In this Letter, we scrutinize this picture with hybrid simulations of a periodic box of background gas in which a continuous hot, slowly drifting CR distribution is injected. To mimic a fluid element close to the Galactic disk, we let the waves escape our box and study the final state of the system. We find that saturation of the magnetic fields occurs while the drift speed is super-Alfvénic and we observe the typical gas heating that is associated with this damping process.

However, contrary to the standard assumption, our simulations demonstrate that not all of the energy put into resonant magnetic fields is damped in the form of heat to the background gas: NLLD initiates an inverse cascade that grows fields on nonresonant, large scales, a novel effect that might have profound phenomenological implications. Although the fate of these large scale modes is unclear based on our simulations alone (the cascade could proceed in a different way in higher dimensions), it is an intriguing suggestion that low energy CRs can generate fields able to resonantly scatter higher-energy CRs.

It is also worth stressing that in our simulations lower energy particles appear to be harder to scatter, as smallscale modes are damped stronger in NLLD and a pitch angle barrier appears. The arising picture potentially constitutes a paradigm shift in how we understand selfgenerated CR transport.

The details of these phenomena and how they are realized in our Galaxy are beyond the scope of this Letter and will be explored in future work.