Numerical Implementation of Streaming Down the Gradient: Application to Fluid Modeling of Cosmic Rays and Saturated Conduction

The equation governing the streaming of a quantity down its gradient superficially looks similar to the simple constant velocity advection equation. In fact, it is the same as an advection equation if there are no local extrema in the computational domain or at the boundary. However, in general when there are local extrema in the computational domain it is a non-trivial nonlinear equation. The standard upwind time evolution with a CFL-limited time step results in spurious oscillations at the grid scale. These oscillations, which originate at the extrema, propagate throughout the computational domain and are undamped even at late times. These oscillations arise because of unphysically large fluxes leaving (entering) the maxima (minima) with the standard CFL-limited explicit methods. Regularization of the equation shows that it is diffusive at the extrema; because of this, an explicit method for the regularized equation with $\Delta t \propto \Delta x^2$ behaves fine. We show that the implicit methods show stable and converging results with $\Delta t \propto \Delta x$; however, surprisingly, even implicit methods are not stable with large enough timesteps. In addition to these subtleties in the numerical implementation, the solutions to the streaming equation are quite novel: non-differentiable solutions emerge from initially smooth profiles; the solutions show transport over large length scales, e.g., in form of tails. The fluid model for cosmic rays interacting with a thermal plasma (valid at space scales much larger than the cosmic ray Larmor radius) is similar to the equation for streaming of a quantity down its gradient, so our method will find applications in fluid modeling of cosmic rays.

💡 Research Summary

The paper addresses the numerical challenges associated with the “streaming down the gradient” equation, which describes the transport of a scalar quantity (or cosmic‑ray pressure) along the direction of its own gradient. While the equation reduces to a simple constant‑velocity advection equation in regions without local extrema, the presence of maxima or minima makes the problem intrinsically nonlinear because the streaming speed changes sign at those points. Traditional upwind explicit schemes, constrained by the Courant–Friedrichs–Lewy (CFL) condition, generate spurious, grid‑scale oscillations when extrema are present. These oscillations originate from unphysically large fluxes leaving a maximum (or entering a minimum) and propagate throughout the domain without damping, contaminating the solution even at very late times.

To understand the source of the instability, the authors regularize the governing equation by adding a small diffusion term, D∇²q. This regularization reveals that the equation behaves diffusively at extrema, which naturally damps the problematic fluxes. With the regularized form, an explicit scheme becomes stable provided the time step obeys Δt ∝ Δx²/D, a much stricter condition than the usual CFL limit but one that eliminates the oscillations entirely.

The study then explores implicit integration strategies. By treating the nonlinear streaming speed implicitly (via fixed‑point iteration and Newton–Raphson linearization), the authors find that the method can remain stable with a time step scaling as Δt ∝ Δx, a significant relaxation compared with the explicit case. However, contrary to the common belief that implicit schemes are unconditionally stable, the paper demonstrates that if Δt is taken too large the method still fails to converge, leading to numerical blow‑up. This failure is traced back to the same non‑linear sign‑switching of the streaming speed that destabilizes the explicit method.

Beyond the numerical analysis, the paper reports several novel physical features of the streaming equation’s solutions. Starting from smooth initial profiles (e.g., Gaussian or sinusoidal), the evolution produces non‑differentiable “kinks” at points where the gradient changes sign. Moreover, the solution develops extended tails that transport material over distances far larger than the initial width of the profile. These characteristics are absent in pure diffusion or constant‑velocity advection problems and reflect the unique physics of cosmic‑ray streaming, where particles move at the Alfvén speed down their pressure gradient.

The authors argue that the fluid model for cosmic rays interacting with a thermal plasma—valid on scales much larger than the cosmic‑ray Larmor radius—is mathematically identical to the streaming‑down‑gradient equation studied here. Consequently, the regularized explicit and carefully time‑stepped implicit schemes presented are directly applicable to large‑scale astrophysical simulations that aim to capture cosmic‑ray feedback, galactic wind driving, or interstellar medium heating. The paper concludes by emphasizing that proper handling of extrema, appropriate diffusion regularization, and judicious choice of time step are essential for obtaining physically accurate and numerically stable solutions in any fluid model that incorporates gradient‑driven streaming.