dHybridR: Hybrid PIC for Relativistic Ions

• dHybridR is a hybrid particle-in-cell code that combines fluid electrons with fully relativistic kinetic ions to model collisionless plasmas.
• It employs the Darwin approximation to eliminate light waves, enabling efficient simulations of cosmic-ray streaming, collisionless shocks, and turbulence.
• The code has been validated through benchmarks showing consistent non-thermal ion acceleration and effective reproduction of key plasma instabilities.

Searching arXiv for papers on dHybridR and related uses. dHybridR is a hybrid particle-in-cell code for collisionless plasmas that combines fluid electrons with kinetic ions whose dynamics can remain fully relativistic, while retaining the radiation-free Darwin approximation for Maxwell’s equations. It was designed for problems in which a small population of energetic non-thermal ions interacts with an otherwise non-relativistic plasma, including cosmic-ray streaming instabilities, collisionless shocks, magnetic reconnection, and turbulence. In that sense, dHybridR extends earlier non-relativistic hybrid PIC methods into a regime where thermal and energetic ions can be evolved self-consistently from sub-relativistic to relativistic energies without the cost of full relativistic PIC (Haggerty et al., 2019).

1. Origin and scientific scope

dHybridR was introduced as a hybrid PIC code with fluid electrons and both thermal and energetic ions that retain relativistic dynamics. Its target regime is one in which the bulk plasma is non-relativistic, but a low-density population of energetic ions can be relativistic and can still affect the overall dynamics of the system. The method paper explicitly frames the code around astrophysical and space-physics problems such as cosmic-ray-driven instabilities, collisionless shocks, magnetic reconnection, and turbulence, and demonstrates applications to linear resonant and non-resonant cosmic-ray streaming instabilities and to strongly non-linear parallel shocks (Haggerty et al., 2019).

The scientific motivation is the gap between standard hybrid codes and full relativistic PIC. Traditional hybrid schemes treat ions kinetically and electrons as a fluid, which is efficient at ion scales but cannot correctly follow ion dynamics once the ion population becomes relativistic. Full relativistic PIC can capture that regime, but at far higher computational cost because electron scales and electromagnetic radiation must also be resolved. dHybridR occupies the intermediate regime: the background plasma remains non-relativistic, electrons remain a charge-neutralizing fluid, yet energetic ions are advanced with relativistic dynamics.

The original demonstrations were chosen to make that scope explicit. The code produced the first self-consistent hybrid runs showing the acceleration of relativistic ions at non-relativistic shocks, and it also provided examples of 2D runs relevant for fast shocks in radio supernovae and 3D runs of low-Mach-number heliospheric shocks that can be compared with in-situ spacecraft observations (Haggerty et al., 2019).

2. Hybrid-relativistic plasma model

The physical model is a hybrid decomposition of the plasma. Ions are represented by PIC macro-particles sampling the ion distribution function, while electrons are treated as a massless, charge-neutralizing fluid. Quasi-neutrality is imposed in the form $n_e \approx n_i \equiv n$, and the current is written as

$\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$

Neglecting electron inertia reduces the electron momentum equation to a generalized Ohm’s law. In terms of ion bulk velocity and current, the electric field is

$$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$

The electron pressure is prescribed through an isotropic polytropic closure $P_e \propto n^{\gamma_{\rm eff}}$. In shock applications, $\gamma_{\rm eff}$ is chosen to enforce approximate downstream equipartition between electrons and ions, whereas in cosmic-ray streaming simulations an adiabatic $\gamma_{\rm eff}=5/3$ is used (Haggerty et al., 2019).

The field evolution adopts the Darwin approximation, so the displacement current is neglected in Ampère’s law: $$\nabla\times\mathbf{B} = \frac{4\pi}{c}\mathbf{J}, \qquad \frac{\partial \mathbf{B}}{\partial t} = -c\nabla\times\mathbf{E}.$$ This removes light waves from the system and preserves the ion-scale focus of the hybrid method. The approximation remains valid provided

$$\left(\frac{V}{c}\right)^2 \ll 1, \qquad \frac{V v_A}{c^2}\ll 1,$$

with $V$ a characteristic bulk speed and $v_A$ the Alfvén speed. The code is therefore intended for systems with non-relativistic bulk flow and $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$0, even if a subset of ions has $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$1.

Its defining extension beyond earlier hybrid codes is the ion equation of motion. Ions obey the fully relativistic Lorentz force,

$\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$2

so the same ion population can be followed from thermal energies to relativistic energies within one simulation (Haggerty et al., 2019).

3. Numerical architecture and normalization

dHybridR is built on the earlier dHybrid code and retains the standard hybrid-PIC separation between particles and fields, but replaces the Newtonian ion advance with a relativistic Boris pusher. The magnetic field is updated with a two-step Lax-Wendroff scheme applied to Faraday’s law. The workflow is therefore conventional at the hybrid-PIC level—deposit moments to the grid, compute $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$3, evaluate $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$4 from Ohm’s law, advance $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$5, and push particles—but the particle pusher uses relativistic momentum $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$6 rather than Newtonian velocity (Haggerty et al., 2019).

The normalization is hybrid-standard and tied to a reference magnetic field $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$7 and density $\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$8. Length is normalized to the ion inertial length

$\mathbf{J} = en(\mathbf{V}_i - \mathbf{V}_e).$9

time to the inverse cyclotron frequency

$$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$0

and velocity to

$$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$1

The electric field normalization is $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$2. In many runs the upstream plasma is initialized with $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$3, $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$4, and $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$5, so that $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$6.

The code supports 1D, 2D, and 3D geometries. In shock applications the standard setup is a reflecting-wall problem: upstream plasma flows toward a conducting wall, reflects, and forms a shock that propagates upstream. Streaming-instability tests use periodic boundaries. The timestep is chosen so that the fastest ions move no more than about one grid cell per step, which becomes a direct constraint on $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$7 once relativistic ions with $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$8 are present. Example production runs used resolutions of 2 cells per $$\mathbf{E} = -\frac{\mathbf{V}_i}{c}\times\mathbf{B} + \frac{\mathbf{J}}{en c}\times\mathbf{B} - \frac{1}{en}\nabla P_e.$$9, with domains ranging from $P_e \propto n^{\gamma_{\rm eff}}$0 in streaming tests to $P_e \propto n^{\gamma_{\rm eff}}$1 along the shock normal in long parallel-shock calculations (Haggerty et al., 2019).

4. Benchmark physics and core demonstrations

The method paper validated dHybridR on both linear and strongly non-linear problems. In the cosmic-ray streaming problem, the code reproduced both the non-resonant Bell instability and the resonant streaming instability. For the non-resonant regime, the fastest-growing mode and growth rate were found at

$P_e \propto n^{\gamma_{\rm eff}}$2

while in the resonant regime the fastest-growing modes satisfied

$P_e \propto n^{\gamma_{\rm eff}}$3

Spectral analysis of $P_e \propto n^{\gamma_{\rm eff}}$4 and $P_e \propto n^{\gamma_{\rm eff}}$5 confirmed that the simulated $P_e \propto n^{\gamma_{\rm eff}}$6 and $P_e \propto n^{\gamma_{\rm eff}}$7 agreed with linear theory in both regimes.

The shock calculations established the code’s central result. In parallel shocks, downstream ions develop a non-thermal tail with

$P_e \propto n^{\gamma_{\rm eff}}$8

which is the canonical strong-shock diffusive-shock-acceleration spectrum in momentum space. When mapped into energy space using the relativistic relation

$P_e \propto n^{\gamma_{\rm eff}}$9

the same population becomes a broken power law: $\gamma_{\rm eff}$0 in the non-relativistic regime and $\gamma_{\rm eff}$1 once $\gamma_{\rm eff}$2. The simulations therefore showed directly that a non-relativistic shock can generate a nearly fixed power law in momentum while producing a steepening in energy around the ion rest mass.

The same runs quantified acceleration. The maximum energy $\gamma_{\rm eff}$3 grows approximately linearly in time, but its growth rate decreases by about a factor of two once particles become relativistic, consistent with a change in the upstream cosmic-ray current and in the effective diffusion coefficient. The acceleration efficiency,

$\gamma_{\rm eff}$4

reached about $\gamma_{\rm eff}$5 and remained approximately constant throughout the simulations. By contrast, quasi-perpendicular shock runs heated ions to supra-thermal energies but did not produce a significant diffusive-shock-acceleration tail, confirming the difficulty of thermal injection in highly oblique geometry. In 3D low-Mach-number heliospheric-shock runs, dHybridR also produced downstream shock rippling and a self-generated $\gamma_{\rm eff}$6 component, permitting direct comparison with in-situ probe-like diagnostics (Haggerty et al., 2019).

5. Later use as a research platform

Subsequent studies used dHybridR as a simulation engine for substantially more elaborate analyses. In fully three-dimensional $\gamma_{\rm eff}$7-$\gamma_{\rm eff}$8 simulations of non-relativistic, oblique collisionless shocks, the code resolved a rippled shock surface, complex ion distributions in the foot and ramp, and a dominant local mode with

$\gamma_{\rm eff}$9

which was identified on the Alfvénic branch as an Alfvén ion cyclotron instability. The same work combined dHybridR output with an instability-isolation procedure and field-particle correlations to separate steady-state shock energization from instability-driven energization, finding the latter to be more than two orders of magnitude smaller than the former while still leaving distinct velocity-space signatures (Brown et al., 2022).

The code was also used in $\gamma_{\rm eff}=5/3$0-dimensional hybrid-kinetic simulations of decaying, supersonic, non-relativistic turbulence. In that setting, dHybridR captured strongly compressible turbulence, density clumping, and non-thermal ion acceleration without a large-scale shock. After accounting for compressibility with the density-weighted velocity field, the turbulent spectrum showed inertial-range slopes of $\gamma_{\rm eff}=5/3$1 at low Mach number and $\gamma_{\rm eff}=5/3$2 at high Mach number. The ion energy spectrum developed a power-law tail $\gamma_{\rm eff}=5/3$3 with $\gamma_{\rm eff}=5/3$4 in non-relativistic energy, and the acceleration efficiency reached values comparable to the $\gamma_{\rm eff}=5/3$5 efficiencies familiar from hybrid shock simulations (Gootkin et al., 22 Sep 2025).

These later studies indicate that dHybridR has been used not only to simulate ion-scale plasma dynamics but also to support advanced post-processing frameworks, including wavelet-based mode identification, local dispersion analysis, and phase-space energization diagnostics. That role depends directly on its ability to provide self-consistent field and particle data over large domains and long integration times.

6. Resolution requirements, strengths, and limitations

Later work on the ion Weibel instability clarified an important numerical constraint for hybrid simulations performed with dHybridR-like massless-electron physics. In weakly magnetized regimes, the dominant Weibel mode is recovered provided that the relevant ion-scale wavelength is resolved; the minimum spatial resolution was summarized as

$\gamma_{\rm eff}=5/3$6

where $\gamma_{\rm eff}=5/3$7 is the number of grid cells per $\gamma_{\rm eff}=5/3$8. The same analysis showed that excessive resolution is also problematic, because it introduces unphysical small-scale whistler modes inherent to the massless-electron approximation. The recommended resolution band was therefore

$\gamma_{\rm eff}=5/3$9

and the comparison with full PIC indicated that massless-electron hybrid models remain reliable for Weibel-mediated shocks only within that window and only up to $$\nabla\times\mathbf{B} = \frac{4\pi}{c}\mathbf{J}, \qquad \frac{\partial \mathbf{B}}{\partial t} = -c\nabla\times\mathbf{E}.$$0 (Orusa et al., 6 Apr 2026).

The principal strengths of dHybridR follow from the hybrid formulation itself. It retains ion kinetic fidelity—reflection, streaming, ring distributions, non-Maxwellian tails, and relativistic particle transport—while avoiding electron kinetic scales and light waves. This makes possible very large domains, long runs over thousands of cyclotron times, and direct simulation of ion acceleration from thermal to relativistic energies. The original method paper’s shock calculations and later turbulence and shock-ripple studies all rely on that separation of scales (Haggerty et al., 2019).

Its limitations are equally structural. Electrons are a fluid, so electron-scale instabilities, non-adiabatic electron heating, and explicitly kinetic electron acceleration are outside the model. The Darwin approximation excludes electromagnetic radiation and requires non-relativistic bulk flow and $$\nabla\times\mathbf{B} = \frac{4\pi}{c}\mathbf{J}, \qquad \frac{\partial \mathbf{B}}{\partial t} = -c\nabla\times\mathbf{E}.$$1, even though a minority population of ions may be relativistic. The 2026 Weibel analysis adds a further practical limitation: because the massless-electron approximation creates an unphysical whistler branch at sufficiently high $$\nabla\times\mathbf{B} = \frac{4\pi}{c}\mathbf{J}, \qquad \frac{\partial \mathbf{B}}{\partial t} = -c\nabla\times\mathbf{E}.$$2, refinement alone does not monotonically improve fidelity in all regimes (Orusa et al., 6 Apr 2026). These constraints define dHybridR as a specialized but powerful tool for ion-scale, collisionless-plasma problems in which relativistic ions matter dynamically but the background plasma remains non-relativistic.