title: Kinetic Simulations of Cosmic-Ray-Modified Shocks I: Hydrodynamics
authors: Colby C. Haggerty Damiano Caprioli
text: 
  
  
  
  Colby C. Haggerty
  chaggerty@uchicago.edu
  
  0000-0002-2160-7288]Colby C. Haggerty
  0000-0003-0939-8775]Damiano Caprioli
  Department of Astronomy and Astrophysics, University of Chicago, 5640 S Ellis Ave, Chicago, IL 60637, USA
  
  
  
  
  Collisionless plasma shocks are efficient sources of non-thermal particle acceleration in space and astrophysical systems. 
  We use hybrid (kinetic ions – fluid electrons) simulations to examine the non-linear feedback of the self-generated energetic particles (cosmic rays, CRs) on the shock hydrodynamics.
  When CR acceleration is efficient, we find evidence of both an upstream precursor, where the inflowing plasma is compressed and heated, and a downstream postcursor, where the energy flux in CRs and amplified magnetic fields play a dynamical role.
  For the first time, we assess how non-linear magnetic fluctuations in the postcursor preferentially travel away from the shock at roughly the local Alfvén speed with respect to the downstream plasma.
  The drift of both magnetic and CR energy with respect to the thermal plasma substantially increases the shock compression ratio with respect to the standard prediction, in particular exceeding 4 for strong shocks.
  Such modifications also have implications for the spectrum of the particles accelerated via diffusive shock acceleration, a significant result detailed in a companion paper, <cit.>, <https://arxiv.org/abs/2009.00007>.
  
  
  
  
  
  
  § INTRODUCTION
  
  Non-relativistic shocks are abundant in space and astrophysical systems, such as the Earth's bow shock, interplanetary shocks associated with coronal mass ejections, supernovae (SNe) and supernova remnants (SNRs), and galaxy clusters;
  they are typically associated with non-thermal particles and emission.
  Such shocks occur on length scales where collisions are too infrequent to efficiently dissipate the energy flux;
  in fact, they are referred to as collisionless shocks because the energy dissipation is regulated by the collective interactions between the charged particles and the electromagnetic fields.
  Despite this, much of our understanding and predictions for the macroscopic structure of plasma shocks comes from the (collisional) fluid Rankine–Hugoniot jump conditions, which cannot accurately model the dynamics of non-thermal particles and magnetic fields.
  Non-relativistic shocks can transfer as much as 10% – 20% of their ram energy into energetic particles (hereafter, cosmic rays, CRs), through repeated first-order Fermi reflections, a process referred to as diffusive shock acceleration <cit.>.
  Additionally, CRs drive plasma instabilities upstream of the shock, which  amplify magnetic fluctuations and produce turbulent, large-amplitude magnetic fields, which are further compressed in the downstream region <cit.>.
  
  Previous theoretical works have included CRs in modeling the hydrodynamics of collisionless plasma shocks using two-fluid theory <cit.>, analytical kinetic theory <cit.>, and Monte Carlo numerical approaches <cit.>. 
  Further numerical approaches have also been employed to understand the time dependence of this effect on different astrophysical shocks <cit.>. 
  For detailed reviews of CR-modified shocks, readers can refer to <cit.>.
  
  When even a relatively small fraction of shock energy (∼ 10%) is channeled into CRs, an upstream precursor is formed and the speed and compressibility of the shock are affected at the zero-th order.
  The standard prediction is that the CR contribution should lead to a larger total shock compression ratio, possibly as large as 10-100 for strong shocks <cit.>.
  The effect of self-generated, amplified magnetic fields has been predicted to limit such a compression to values ≲ 10 <cit.>, but always greater than the fiducial value of 4 for shocks with Mach number ≫ 1.
  While there is general agreement that self-generated CRs and magnetic fluctuations should modify the hydrodynamics of a shock, the process through which this occurs has not been assessed from first principles, yet.
  
  In this paper we present the first numerical evidence of CR-modified shocks in ab-initio plasma simulations in the hybrid limit (kinetic ions – fluid electrons).
  We detail how CRs and CR-driven magnetic fluctuations affect the hydrodynamics of the shock, while in a companion paper, <cit.>, we discuss how the modified hydrodynamics affect the spectrum of the particles accelerated via DSA.
  
  The layout of this paper is as follows: in <ref>, we detail the hybrid code and the shock simulation setup.
  In <ref> and <ref>, we identify and discuss the formation of a precursor region in the upstream and a postcursor region in the downstream, respectively. 
  In <ref>, we present CR modified shock jump conditions and show they are in good agreement with simulations.
  Finally, in <ref>, we discuss the implication of these corrections to the fluid dynamics of collisionless plasma shocks.
  
  
  
  § HYBRID SIMULATIONS
  
  To study the non-linear effects of self-generated CRs on the hydrodynamics of the shock, we perform self-consistent simulations using , a relativistic hybrid code with kinetic ions and (massless, charge-neutralizing) fluid electrons <cit.>.
   is the generalization of the Newtonian code dHybrid <cit.>, which has been already widely used for simulating collisionless shocks <cit.>.
  Hybrid codes are better suited to self-consistently simulate the long-term shock evolution than full particle-in-cell codes because they do not need to resolve the small time/length scales of the electrons, which are usually dynamically negligible  <cit.>. 
  
  All physical quantities are normalized to their far upstream values, namely: mass density to ρ_0 ≡ m_i n_0 (with m_i the ion, namely proton, mass), magnetic fields to B_0, lengths to the ion inertial length = c/ω_pi (with c the speed of light and ω_pi the ion plasma frequency), time to the inverse ion cyclotron frequency , and velocity to the Alfvén speed 0 = B_0/√(4πρ_0).
  The ion temperature is chosen such that the thermal gyroradius is 1, corresponding to an ion thermal to magnetic pressure ratio of β_i = 2.
  The system is 2D in real space (in the x-y plane), and all the three components of momenta and electromagnetic fields are retained.
  The hybrid model requires an explicit choice for the electron equation of state, and in this work, electrons are assumed to be adiabatic, i.e., the electron pressure is P_e∝ρ^5/3.
  We remark on the motivation and consequences of choosing this equation of state in <ref>.
  
  
  The simulations are initialized with a uniform magnetic field B_0=B_0  x and with a thermal ion population with a bulk flow v_x = - M  x). 
  The simulations are periodic in the y direction, the right boundary is open and continuously injecting thermal particles, and the left boundary is a reflecting wall;
  after tens of cyclotron times, the ion population closest to the wall reflects, becomes unstable, and forms a shock.
  This shock then travels in the + x direction, parallel to the background magnetic field (parallel shock), with the downstream plasma at rest in the simulation reference frame.
  Note that M defines both sonic and Alfvénic Mach numbers M ≡ M_A = √(γ)M_s in the downstream (simulation) frame, while the Mach number that enters the stationary jump conditions is the one in the shock frame, which is 20-30% larger, depending on the shock compression ratio (see <ref>).
  Most of the analysis in this work is done on a benchmark simulation with M = 20 and a domain of size [L_x,  L_y] = [10^5, 200], wide enough to account for 2D effects and long enough so that the simulation could be run for more than 1000 without energetic particles escaping the box.
  The simulation has two grid cells per , and each grid cell is initialized with four particles per grid.
  The speed of light is set to be much larger than the Alfvén and thermal speeds (c/0 = c/v_thi = 100), as discussed in <cit.>; the time step is set as cΔ t = d_i/10. 
  
  Simulations with Mach numbers in the range M=10 - 80 were performed;
  their parameters are detailed in Appendix A. 
  All of the following analysis and figures in this work are preformed with the M=20 benchmark simulation, unless stated otherwise.  
  The choice stems from the fact that a M=20 simulation is representative of a strong shock (M ≫ 1) and can still be run for a long time without CRs escaping the simulation domain; 
  the computational cost of a  simulation scales ∝ M^2 when keeping the box size fixed in units of the ion gyroradius (for a Newtonian hybrid code it would scale ∝ M^3, since also dt is inversely proportional to M if not normalized to c.) 
  
  Moreover, for M = 20 the fastest growing modes, driven by streaming CRs, are in the resonant regime <cit.>;
  this means that the amplified fields should be Alfvénic and remain quasi-linear upstream, which simplifies the theoretical interpretation.
  We have also performed runs with larger M (<ref>), for which magnetic field amplification occurs in the Bell regime, finding consistent results for all strong shocks.  
  
  A shorter M=20 simulation with a wider box L_y = 10^3  was run for assessing convergence with the transverse size; quantities from this simulation are shown in Figure <ref>.
  With these parameters, strong parallel shocks can channel as much as 10 - 15% of their kinetic energy into CRs and effectively amplify magnetic fields <cit.>.
  
  
  
  § THE UPSTREAM PRECURSOR
  
  
  
   §.§ Hydrodynamics
  
  
  When a large enough fraction of the shock kinetic energy and pressure are deposited into energetic particles, the standard shock hydrodynamics are altered from the ones described by the gaseous Rankine-Hugoniot jump conditions.
  CRs break the causality wall represented by the shock, and thus, can transfer momentum and energy back upstream, effectively slowing the shock and pre-compressing the incoming plasma, forming a CR-induced precursor <cit.>. 
  For the quantities measured in this work, we use the subscripts 0, 1 and 2  for the far upstream, the precursor immediately upstream of the shock, and downstream respectively.
  
  The formation of the CR precursor is evident in our simulations:
  the upper panel of Figure <ref> shows the x-v_x phase space distribution function integrated/averaged over the y, v_y and v_z directions at 1400. Cuts of the ion v_x distribution are shown in the lower panel. Each cut is taken at different distances from the shock (color coded) and averaged over t = 1350 - 1550.
  The bulk flow speed (i.e., the mean of the v_x distribution) upstream of the subshock is u_1∼ 0.9 u_0, consistent with pressure in CRs being ∼ 10% of the bulk ram pressure.
  In the precursor, thermal ions warm up (i.e., the v_x distributions become wider) well beyond what is predicted by adiabatic heating due to the observed compression ρ_1/ρ_0∼ 1.09[Mass flux conservation implies that ρ u is constant.].
  A temperature increase consistent with adiabatic heating would yield T_1/T_0 = (ρ_1/ρ_0)^2/3≈ 1.06, while Figure <ref> indicates a much larger temperature increase, T_1/T_0∼ 2.7.
  Such a non-adiabatic heating is due to the damping of the waves produced by CR-driven instabilities <cit.> and corresponds to maintaining a constant thermal/magnetic pressure ratio in the shock precursor, despite the effective magnetic field amplification. 
  
  
  
   §.§ Magnetic Field Amplification
  
  
  
  The CR population is intrinsically anisotropic in the upstream fluid frame, which leads to the development of several plasma instabilities that  have been investigated theoretically and numerically in numerous works <cit.>.
  The presented kinetic simulations capture self-consistently both the formation of the precursor and the CR-driven magnetic fluctuations.
  
  The power spectrum of self-generated magnetic fluctuations in the upstream is shown in Figure <ref> for our benchmark run. 
  The spectrum is calculated using the transverse components of the magnetic field with |B_⊥|^2 ≡ |B_y|^2+|B_z|^2, where B_i is the Fourier transform of B_i(x,t), as a function of wave number k and angular frequency ω[In practice, we calculate the discreet Fourier transform of B_i(x,y,t) in the x and t directions, average over the y direction, and normalize to the number of grid points in x and t].
  Also, it is evaluated over a relatively large window [60000, 200], embedded in the upstream flow, with the left side of the window beginning at x = 23600 at t = 600 and ending at x = 8000, just in front of the shock, at t = 1780.
  The phase velocity of an Alfvén wave, based on the initial upstream magnetic field and density is shown by the white dashed lines, with the upper and lower lines corresponding to waves traveling  towards (leftward) and away from (rightward) the shock respectively.
  Most of the magnetic power is in modes with phase speeds comparable to the initial upstream Alfvén speed and with ω < 0, i.e., modes that are propagating away from the shock in the fluid rest frame and down the CR pressure gradient.
  
  Furthermore, most of the power is in modes with wave numbers between kd_i ∼ 0.002 - 0.02, which are resonant with the CR population upstream of the shock.
  The resonant CR population can be seen in the bottom panel of Figure <ref>, which shows a 2D map of the ion energy spectrum, Ef(E), integrated over the transformation window. The spectrum is plotted as a function of time and inverse momentum (m_i0/p), which, for a constant magnetic field strength, can be interpreted as the corresponding resonant wave number. 
  Note that the CR spectrum is cut off for values smaller than the injection momentum (, vertical dashed line), so that the signal is not overpowered by the inflowing beam of ions with p = 20m_i0 <cit.>.
  By comparing the x-axis of the energy spectrum and magnetic power spectrum, it can be concluded that, for the benchmark simulation, the amplified modes are driven by the resonant CR streaming instability, <cit.>;
  this result is in agreement with theoretical expectations <cit.> and previous hybrid simulations <cit.>, which suggest that the non-resonant (Bell) instability should become dominant only for M≳ 30.
  
  
  
   §.§ The Quasi-periodic Nature of the Precursor
  
  
  
  While the standard predictions for the precursor are recovered in time-averaged simulation data, a high-cadence analysis shows that the precursor varies significantly on intermediate time scales, around 50 - 200.
  These time scales roughly match the period of the Doppler-shifted Alfvénic fluctuations in the precursor[Most of the power in resonant modes are at wavelengths of λ∼ 2π/k ∼ 2000 (Figure <ref>).
   With a bulk flow of ∼ 200, the period of these oscillations in the shock frame are ∼λ/(200) ∼ 100.].
   
  A global time variability originates from the CR-driven instabilities and the spatially inhomogeneous amplification of the upstream magnetic field. 
  In fact, the effective shock inclination, (t), is modulated on time-scales comparable to the wave period. Furthermore, the shock inclination controls particle injection into DSA, which is more prominent for quasi-parallel shock configurations and suppressed for oblique (≳ 50) ones.
  
  This modulation can be seen in Figure <ref>, where the dark green line (left axis) shows (t)≡cos^-1( <B_x/|B| >), with the magnetic field averaged over a region between 350 and 450  upstream of the shock.
  Sufficiently far upstream the field is mostly along its initial direction, but (t) oscillates around its average of ≈ 48^∘ (dashed line) with a period comparable to that of the self-generated waves.
  
  Figure <ref> also shows the transverse CR pressure, defined as P_yy + P_zz, where P_ii = ∫_^∞ p_iv_i f d^3p in units of the upstream thermal pressure P_0 (light blue, right axis), calculated over the same region as .
  Such a CR pressure varies with a period comparable to the shock inclination too, but the two quantities appear to be generally anti-correlated.
  We conclude that CR injection, the prominence of the precursor, and the shock inclination are all modulated by the period of the waves generated via CR streaming instability in the precursor.
  The imprint of upstream magnetic fluctuations on the shock dynamics is expected to survive also in more realistic (much wider and 3D) systems, though as a local phenomenon that, arguably, causes shock rippling.
  This quasi-periodic nature of precursors may be measurable in-situ at heliospheric shocks and must be reckoned with when determining the fundamental inclination of a shock based on an instantaneous measurement  .
  
  
  
  
  § THE DOWNSTREAM POSTCURSOR
  
  
  In contrast to the formation of the shock precursor, what was not predicted by any theory (that we are aware of) and what we found for the first time in self-consistent plasma simulations is the formation of a non-linear structure downstream of the shock, which we call a postcursor.
  
  The postcursor manifests itself as an extended downstream region where the self-generated magnetic fluctuations have an important dynamical role and the dynamics of thermal plasma and CRs are decoupled, in the sense that CRs and magnetic perturbations have a sizable drift speed with respect to the bulk plasma.   
  
  
  
  A standard assumption in DSA theory is that CRs are quickly isotropized in the frame of the magnetic fluctuations, usually assumed to be Alfvén waves <cit.>, and that the Alfvén speed is much smaller than the flow speed in the shock frame, so that eventually CRs are isotropic in the flow frame at the order of v_A/u.
  Deviations from such an isotropy may arise if the shock is oblique and trans-relativistic <cit.> or in general in the presence of anisotropic transport <cit.>.
  
  In the presence of magnetic field amplification, the condition v_A/u≪ 1 may be violated: it has been suggested <cit.> that retaining the upstream Alfvénic drift of CRs with respect to the background fluid may also affect the resulting CR spectrum  appreciably  <cit.>.
  Retaining this effect upstream is natural, since self-generated waves travel against the CR gradient, i.e., towards upstream infinity resulting in a net CR drift ≃ u-v_A, as shown in the top panel of Figure <ref>. 
  However, it has always been assumed that downstream magnetic fluctuations should not have a preferential direction, hence canceling any net drift <cit.>.
  
  
  
   §.§ CR Drift Downstream
  
  In this work we find, for the first time, evidence for a net CR drift away from the shock, as shown in Figure <ref>.
  In all the panels, the color code corresponds to time evolution.
  The top two panels show the bulk velocity normal to the shock, for all the ions (thermal + CRs, v_x, top panel) and for the CRs alone ( bottom panel). The speeds are averaged along y and shown as a function of x. 
  The CR velocity is defined as the mean velocity of particles with momentum larger than , i.e.,
  
      .  v_ CR≡∫_^∞ vp^2 f dp 
              / ∫_^∞ p^2f dp .,
  
  where ≈√(10)mu_0 according to the injection theory derived from hybrid simulations in <cit.>.
  
  While the downstream plasma is –as expected– at rest behind the shock in the simulation frame (top panel of Figure <ref>), the velocity of the CR population has a net negative value of about 2-3 0 in the postcursor, about 500 wide just downstream, and reduces in magnitude further downstream until vanishing close to the left wall (second panel).
  Ultimately, the drift of CRs is controlled by their interaction with the magnetic fluctuations, and so, this flux is expected to be tied to the dynamics of the downstream magnetic field.
  
  The third panel of Figure <ref> shows the local Alfvén speed | v_A(x)| =  <| B(x,y)|/√(4πρ(x,y)) >_y, averaged in the y direction (see also the bottom panel of Figure <ref>).
  It is important to stress that this quantity is only informative about the typical speed of the magnetic perturbations, which are not simple Alfvén waves but rather large-amplitude structures created by CR-driven instabilities in the precursor and further compressed at the shock. 
  In the paper we will commonly use “wave” for brevity, bearing in mind that magnetic structures are often nonlinear.
  Finally, the bottom panel of Figure <ref> compares directly the average CR drift velocity and the postcursor Alfvén speed: they are consistently similar in time and approach | v_A| ≈ 3 0 at later times.
  
  In summary, in the postcursor, the  downstream CRs drift away from the shock with respect to the thermal plasma at a speed comparable to the local Alfvén speed;
  the extent of the CR postcursor is determined by the spatial extent of the amplified magnetic field.
  
  
  
   §.§ Motion of Post-shock Magnetic Structures
  
  
  
  Even if there is agreement in the magnitude of CR drift and Alfvén speeds, we are left with the task to demonstrate that the frame in which CRs are isotropic is actually the wave frame.
  Note that in simulations we observe CRs drift away from the shock relative to the thermal plasma, and it is not obvious that  magnetic fluctuations should have a preferred direction of propagation downstream.
  To quantify this, we present two distinct diagnostics: a morphological analysis of the magnetic structures behind the shock and a calculation of the wave dispersion relations.
  
  The evolution of the magnetic field (B_z and | B|) is shown in Figure <ref>: 1D cuts are taken along the x direction, averaged between y = {47.5, 52.5}, and plotted from x = 2500 to 7000 over the last 1200 of the simulation.
  In each panel the shock front makes a nearly straight line with a positive slope corresponding to the inverse of the shock speed;
  the downstream region lies above this line.
  The inclination of the post-shock, finger-like, magnetic structures provides an estimate of their velocity in the simulation frame: just behind the shock, waves have a negative velocity consistent with dashed lines of -2.50, while further downstream they are almost at rest (i.e., vertical).
  This suggests that, in the postcursor, magnetic fluctuations are traveling away from the shock faster than the thermal plasma, consistent with the net CR drift.
  
  
  
  The Fourier analysis of the post-shock magnetic structures shown in Figure <ref> also supports this claim. 
  The transform is calculated the same as in Figure <ref>, and over the last 1200 of the simulation, in a window [500 , 200], positioned just downstream and moving with the shock;
  the angular frequency versus wave number diagram is then boosted back into the downstream frame.
  Most of the power is in modes with ω>0, i.e., modes that move away from the shock with respect to the background plasma, at  speeds ranging between -40 and -10, corresponding to the different dashed lines in Figure <ref> and consistent with the 2D time stack plots in Figure <ref>.
  
  While the dispersion relationship is not as clear as Figure <ref>, due to the non-linear nature of the modes, we can still conclude that the postcursor magnetic structures preferentially move away from the shock, with phase speeds consistent with the Alfvénic and CR drift speeds in Figure <ref>.
  We point out that accounting for standard reflection/transmission of the  quasi-linear Alfvén waves observed upstream <cit.> would return modes with ω< 0 in the downstream, instead. 
  Most likely, the failure of such a linear theory is due to the non-linear, collective, reaction of the shock transition itself to the hammering of the upstream fluctuations. 
  The shock would therefore behave as an antenna that sends magnetic perturbations towards downstream, at a speed comparable to the “natural” speed for a magnetized plasma, namely the Alfvén speed[We thank S. Schwartz for pointing out such an analogy.].
  
  From all of these diagnostics a picture arises where magnetic structures generated upstream, via CR-driven instabilities, are advected and compressed through the shock and retain their inertia over a sizable region downstream, forming the postcursor.
  In the postcursor both the magnetic fluctuations and the CRs, which are isotropized in the wave frame, drift away from the shock faster than the thermal plasma, with a velocity on the order of the local Alfvén speed, 2.
  Since the post-shock magnetic pressure may become of order ≲ 10 of the shock bulk ram pressure, such a drift is not negligible with respect to the downstream fluid velocity;
  for our benchmark run, 2∼ 0.5 u_2.
  
  This extra energy flux in both energetic particles and magnetic fields behind the shock affects the shock hydrodynamics at the zero-th order.
  We present a discussion of the modifications to standard  Rankine-Hugoniot jump conditions necessary to account for this non-linear feedback in the next section.
  
  
  
  § MODIFIED SHOCK HYDRODYNAMICS
  
  
  
   §.§ Theory
  
  In the previous sections, we have illustrated the appearance of two non-linear features in the shock profile: the precursor and the postcursor. 
  Here we address how they modify the shock compression and how energy is partitioned among different components;
  for the effects that such modifications have on the spectrum of the accelerated particles, we refer to the companion paper, <cit.>.
  
  Several works have developed a theory of CR-modified shocks including the dynamical role of energetic particles  <cit.>, as well as of the self-generated magnetic field <cit.>, but all of these studies only accounted for the CR precursor.
  The main effect of the precursor is to produce a weaker subshock with compression ratio ≡ u_1/u_2<4 and a total compression ratio ≡ u_0/u_2>4.
  The actual value of the total compression ratio depends on: 
  1) the softer equation of state of CRs, which have an effective adiabatic index of γ_c=4/3; and 
  2) the escape of CRs from far upstream, which makes the shock behave as partially radiative <cit.>.  
  The former effect would saturate the total compression to 7, but with CR escape,  may become very large.
  
  To understand the role of the postcursor, we consider the conservation equations for mass, momentum, and energy in a 1D, non-relativistic, stationary, shock:
  
      [ ρ u ] = 0 
       [ ρ u^2 + P_g + P_c + P_B] = 0 
       [ 1/2ρ u^3 + F_g + F_c +F_B ] = 0
  
  where γ_i, P_i, and F_i are the adiabatic index, pressure, and energy flux of thermal gas, CRs, and magnetic fields (i=g,c,B, respectively).
  We define the bulk flow velocity as u≡ -u  x, such that both u > 0 and F > 0 even if x points from downstream to upstream, as in the simulations; 
  the square brackets denote the difference between two arbitrary x locations.
  With appropriate prescriptions for F_i, this set of equations can be used to solve for the shock jump conditions.
  
  The thermal gas energy flux has the canonical form:
  
      F_g (x)= γ_g/γ_g -1u(x)P_g(x)
  
  and the CR energy flux is linked to the CR pressure in a similar way, i.e.:
  
      F_c(x) = γ_c/γ_c - 1 u_c(x) P_c(x),
  
  where u_c is the CR bulk speed, including drifts.
  Typically, most of the CR energy is in relativistic particles and γ_c ≈ 4/3 but, for ab-initio simulations of non-relativistic shocks, non-relativistic CRs carry a sizable fraction of the energy <cit.>;
  in our benchmark run, we measure γ_c as the ratio of the enthalpy density to internal energy density and obtain γ_c ≈ 1.5 (see Appendix <ref> for more details).
  
  The “magnetic" energy flux, in general, has contributions from both the magnetic pressure and the kinetic energy associated with the plasma fluctuations, which for Alfvénic fluctuations reads <cit.>:
  
      F_B(x) = [2u_B(x) + u(x)]P_B(x),
  
  where u_B= u ± v_A is the local velocity of the magnetic fluctuations.
  Equation <ref> encompasses the effective equation of state for the magnetic fluctuations, which, in general, cannot be expressed as polytropic; 
  such an equation of state depends on the nature of the fluctuations, and Equation <ref> strictly holds only for Alfvénic perturbations.
  In the absence of a general theory, we adopt the same prescription, even when non-resonant modes dominate, relying on the fact that in its nonlinear stage, Bell-driven turbulence looks quasi-Alfvénic <cit.>.
  
  
  Using the results from Figure <ref> and <ref>, we have that in the precursor u_B≃ u-v_A and in the postcursor u_B≃ u+v_A;
  far upstream u_B,0≃ u_0 because waves have no preferential direction. 
  
  
  
  Considering Equations <ref> and <ref> between 0 (far upstream) and 2 (postcursor), we normalize the momentum(energy) flux equation by dividing by the ram pressure (energy) ρ_0u_0^2, (ρ_0u_0^3/2), introduce the normalized pressure ξ_i≡ P_i,2/(ρ_0u_0^2)  and  η_i≡ 2γ_i/(γ_i-1), obtaining:
  
      ξ_g ≃  1 + 1/γ_g M_s^2 + 1/2 M_A^2- 1/ - ξ_c - ξ_B
  
  and
  
      1    + η_g/γ_g M_s^2 + 3/M_A^2≃
         1/^2 + η_gξ_g/ 
          + u_c,2/u_2η_c ξ_c/
          + (2u_B,2/u_2+1) 2ξ_B/,
  
  where  M_s ≡ρ_0u_0^2/γ_g P_0 and M_A ≡ 4πρ_0u_0^2/B_0^2 are the far upstream sonic and Alfvénic Mach numbers, respectively.
  Both the far upstream CR pressure and escape flux are negligible with respect to all the other terms for moderate ξ_c <cit.>, so we pose both P_c,0≃ 0 and F_c,0≃ 0.
  
  As we saw in <ref>, u_B,2≃ u_c,2≃ u_2+2 and hence:
  
      u_c,2≃ u_B,2 = u_2(1 + √(2ξ_B)).
  
  
  Equations <ref>, <ref>, and <ref> can be rearranged into one equation for  as a function of M_s, M_A and the post-shock pressures ξ_c and ξ_B. 
  This equation —Equation <ref>, detailed in Appendix <ref>— is quartic with respect to √() and its physical solution is shown in the left panel of Figure <ref>, where  is plotted as a function of ξ_B  and ξ_c, for the values of M_s, M_A and η_c from our fiducial simulation. 
  The striking conclusion from Figure <ref> is that –even if just a few percent of the pressure is channelled into CRs and magnetic fields–  becomes larger than the test-particle prediction of 4 (light-green line). 
  Note that, in contrast to the standard non-linear DSA theory without the postcursor, ≳ 4 can be realized even if CRs are non-relativistic (γ_c=5/3) and do not escape upstream (see upper left subplot of Figure <ref> in Appendix <ref>).
  
  In fact, on top of the usual compressibility enhancement due to the presence of both CRs and magnetic field, which have a softer equation of state than the thermal gas, a key role is played by the drift of the non-thermal components in the postcursor, which effectively acts as an energy sink.
  Note that  is more sensitive to ξ_B than ξ_c, as the amplified magnetic field contributes to both the non-thermal energy downstream as well as the rate at which both CR and magnetic energy travel away from the shock.
  
  These non-linear effects reduce the kinetic energy that is dissipated into heat at the shock. 
  This can be seen in the right panel of Figure <ref>, which shows the fractional reduction in the downstream temperature, T_2/T_2(ξ_B = ξ_c = 0), as a function of ξ_B and ξ_c for the values of M_s, M_A and η_c from our fiducial simulation. Again, even if relatively little energy is diverted to the CRs and magnetic fields, the predicted downstream temperature can decrease by as much as 50%. 
  
  
  
  
   §.§ Simulation Comparison
  
  Let us now compare the prediction for  from the modified jump conditions with our benchmark simulation.
  First of all, we need to find the (instantaneous) shock speed and boost all the velocities in that frame, which is the only stationary frame in which Equations <ref> through <ref> are valid;
  this defines the actual value of u_0 =  +||=  + u_2 used for normalization purposes, i.e., 
  
      u_0= M/ -10.
  
  As in <ref>, we calculate ξ_c and ξ_B by averaging over a region 500 long behind the shock; 
  the CR pressure is obtained by taking the appropriate moment of the ion distribution function above  in the downstream frame[This definition assumes that the CR distribution is isotropic in the downstream frame. This differs from the actual CR pressure on the order of m_i2/p_c where p_c is the momentum where most of the energy in the spectra resides, i.e., p_c ∼ for steep spectra and in our simulation m_i2/p_c ∼ 4%]: 
  
      ξ_c = 4π/3ρ_0u_0^2∫_^∞  v(p) f(p) p^3 dp.
  
  
  
  
  ξ_c and ξ_B are plotted as a function of time in the first panel of Figure <ref> (crosses and stars, respectively);
  the color code corresponds to the time in the simulation.
  Together, the normalized CR and magnetic pressure make up about 15-20% of the pressure budget in the postcursor region.
  ξ_c increases quickly to a value just above 0.1 and remains nearly constant throughout the simulation, whereas the magnetic pressure rises more slowly up to a value of 0.05-0.075 towards the end of the simulation.
  The order of magnitude of the normalized pressures in the simulation makes it clear that, based on the predictions from the modified jump condition in Equation <ref>, the total compression ratio should be larger than 4.
  The second panel of Figure <ref> shows the y-averaged density profile as a function of time:  
  at early times (blue lines) the compression ratio ≈ 4, but it increases with time up to ≳ 5 - 6 toward the end of the simulation (red lines).
  
  The agreement between the prediction based on the modified jump conditions and the simulation can be quantified by determining a time-dependent (t), averaged over the postcursor.
  The solution of Equation <ref> for the actual values of ξ_c(t) and ξ_B(t) (top panel) is compared with the measured value of  (diamonds) in the bottom panel of Figure <ref>; the grey-shaded area corresponds to a fiducial 10% error on the prediction that encompasses the uncertainty in measuring pressures and velocities in a profile with small-scale spatial variation, as well as the assumption of Alfvénic-like magnetic turbulence (Equation <ref>). 
  We find a general agreement between theory and simulations in the value, trend, and periodic variation of , with minor deviations that can potentially be attributed to time evolution (not captured by Equations <ref>-<ref>) and transient shock features.
  To stress the importance of the postcursor in the shock dynamics, the bottom panel of Figure <ref> includes as a dashed line the prediction with no CR/magnetic drift (u_c = u_b = u_2).
  As mentioned above, CR and magnetic pressure terms alone are not sufficient to account for the strong shock modification that we observe. 
  
  The predicted values of  can also account for the quasi-periodic behaviour of the shock driven by the variations in ξ_b (top panel of Figure <ref>).
  Such oscillations hinge on the same precursor physics discussed in Section <ref> and will be explored in greater detail in future works.
  
  
  
   §.§ A Critical Review of Previous Results
  
  Such modifications on the shock hydrodynamics  were not seen in previous hybrid campaigns that covered similar parameter space <cit.> because of the choice of the effective electron polytropic index, γ_e.
  When an effective γ_e is chosen in order to mimic temperature equilibration and pressure balance between post-shock electrons and ions, a large value γ_e≈ 3-4 must be chosen for strong shocks <cit.>. 
  However, enforcing electron-ion equilibration via an effective γ_e requires assuming a priori the value of the realized compression ratio; 
  therefore, if one assumes that ≈ 4, the artificially-stiff electron equation of state ends up limiting the total compression that can be achieved. 
  Using such a prescription, <cit.> found ∼ 4.4, underestimating the shock modification that we report here with the adiabatic equation of state.
  Note that as long as the post-shock pressure does not become dominated by electrons (which is unphysical but may happen for γ_e≫5/3), the shock modification does not depend on whether electrons are adiabatic or in equipartition with the ions.
  
  Very recently <cit.>  pointed out how standard MHD jump conditions are not commonly realized in many papers based on PIC simulations of shocks; 
  our formalism, in which the kinetic backreaction of accelerated particles and self-generated magnetic fields is accounted for, naturally presents itself with a physically-motivated framework in which to study the agreement between simulations and theory.
  
  
  
   §.§ Mach Number Dependence
  
  
  Beyond our benchmark M = 20 simulation, we have run parallel shocks with M=10,  40 and 80, and the results are shown in Figure <ref>.
  As already reported by <cit.>, the normalized CR pressure —a proxy for acceleration efficiency— is generally about 10%; 
  the normalized magnetic pressure, instead, is typically 2-5%, values commonly inferred in multi-wavelength analysis of SNRs <cit.>.
  
  Consistent with our predictions, the shocks in each of these simulations have compression ratios that exceed the gaseous value. 
  This enhancement is shown in the bottom panel of Figure <ref>, in which the y-averaged density profiles are shown. 
  The color-coded shaded regions correspond to the prediction of  from Equation <ref> and are in good agreement with the simulated values for each Mach number, which strengthens the applicability of these results to many heliospheric/astrophysical systems.
  
  At larger Mach numbers M≳ 30 the magnetic field amplification in the precursor is controlled by the growth of the Bell instability  <cit.>, whose fastest growing modes are purely growing (i.e., with almost zero phase speed) and with wavelengths much shorter than the CR gyroradii.
  Nevertheless, global hybrid simulations of shocks with M=60,80 and 100 have shown that small-wavelength modes saturate rather quickly far upstream, and that in the precursor most of the power in self-generated fields is still in modes quasi-resonant with the accelerated particles <cit.>.
  In the non-linear stage of the Bell instability waves start propagating with a phase speed close to the Alfvén speed in the amplified field <cit.>, so we expect that the general phenomenology outlined in this work should apply also at very strong shocks, such as those in SNRs.
  
  
  
  
  § CONCLUSIONS
  
  In this work we use self-consistent hybrid simulations to study the modifications that self-generated CRs and associated magnetic turbulence induce on the dynamics of a collisionless plasma shock.
  The efficient acceleration of CRs leads to a pre-compression and deceleration of the plasma upstream of the shock: 
  in such a precursor region, there is non-adiabatic heating of the inflowing plasma, likely a byproduct of the CR-driven magnetic fluctuations <cit.>. 
  Additionally, we identify quasi-periodic fluctuations in the magnetic field strength and CR pressure in the precursor, which are attributed to the shock geometry transitioning back and forth between quasi-parallel and oblique/quasi-perpendicular configurations.
  
  Here, we report for the first time the formation of a characteristic region downstream of the shock, which we call the postcursor, where the enhanced magnetic fluctuations generated upstream and then compressed by the shock play a crucial dynamical role.
  Such magnetic fluctuations are found to propagate away from the shock in the downstream rest frame, with a velocity comparable to the local Alfvén speed in the amplified magnetic field (Figures <ref> and <ref>).
  In the postcursor, CRs become isotropic in a frame moving with the magnetic fluctuations, rather than in the downstream fluid frame, resulting in a peculiar drift with respect to the background plasma (Figure <ref>).
  
  Such non-linear features, ultimately driven by CR physics, lead to an enhanced shock compression ratio.
  More precisely, the total compression ratio becomes larger than the standard value of =4, due, not only to the compressibility of relativistic CRs and magnetic fields, but mainly because of the larger rate at which the non-thermal populations are advected away from the shock.
  The solution of the modified jump conditions (Equation <ref>) is presented and compared with simulations proving general agreement between the predicted and measured compression ratio as a function of time (Figure <ref>).
  Even a moderate CR acceleration efficiency, ξ_c∼ 10%, is able to increase the shock compression ratio for a large Mach number shock by nearly 50% from the standard fluid prediction to ∼ 5.5.
  
  
  We have tested that this behaviour is not limited to our benchmark case by running simulations with M=10,  40, and 80, which present a very similar phenomenology. 
  While it is computationally very challenging to run kinetic simulations of strong shocks for much longer than we did, we do not see strong evolution of the shock modification in the last few hundreds  and we achieve quite large values for ξ_c and ξ_B, which may suggest that in realistic shocks the compression ratio should not be much greater than ≳ 6. 
  Provided that a reliable prescription for CR injection were implemented, the long-term shock evolution could be followed with hybrid+MHD codes, in which thermal particles are described as a fluid <cit.>.
  
  The enhanced compression ratios found in these hybrid simulations can be regarded as the first ab-initio evidence of the existence of CR-modified shocks, which had been suggested almost 40 years ago <cit.> but never verified in kinetic simulations.
  Observational hints of shock compression ratios larger than 4 have been reported for young supernova remnants, such as Tycho  <cit.> and SN1006 <cit.>;
  in particular, in SN1006 the distance between the forward shock and the contact discontinuity is inferred to be modulated with the azimuth, being smaller (corresponding a larger compression ratio) where the shock is quasi-parallel <cit.>, i.e., the region where CR acceleration is expected to be more prominent <cit.>.
  
  The values of ξ_c∼ 5-15% and ξ_b∼ 2-10% required to produce ≳ 4 are  consistent with the values inferred from multi-wavelength observations of young SNRs <cit.>, so we expect CR-induced shock modification to be a ubiquitous phenomenon. 
  
  An important result that follows from our simulations is that the main driver of the shock modification, i.e., the magnetic drift in the postcursor becoming comparable to the local Alfvén speed, can be inferred observationally.
  In fact, the postshock magnetic field can be constrained with non-thermal X-rays <cit.> and the gas density and temperature by X-ray <cit.> and Balmer line emission <cit.>, while the shock speed can be estimated from proper motion of X-ray and/or Balmer features.
  In principle, high-resolution X-ray observations can even test the presented theory as a function of the local shock inclination, e.g., in bilateral SNRs such as SN1006, probing whether CR-modified shocks manifest themselves in quasi-parallel regions.
  A corollary of our results is that quasi-perpendicular shocks, which are generally poor ion injectors <cit.>, should not exhibit deviations from standard Rankine–Hugoniot conditions. 
  
  A natural question that arises is what is the spectrum of the particles accelerated in a CR-modified shock? 
  For a dedicated discussion of such a critical question we refer to a companion paper, <cit.>.
  
  
  
  
   <cit.>, 
  
  
  
  We would like to thank  P. Blasi, E. Amato, A. Spitkovsky, D. Eichler, L. O'C. Drury, S. Schwartz, and L. Wilson III for stimulating and constructive discussions. 
  This research was partially supported by NASA (grant NNX17AG30G, 80NSSC18K1218, and 80NSSC18K1726), NSF (grants AST-1714658, AST-1909778, PHY-1748958, PHY-2010240), and by the International Space Science Institute’s (ISSI) International Teams program.
  Simulations were performed on computational resources provided by the University of Chicago Research Computing Center, the NASA High-End Computing Program through the NASA Advanced Supercomputing Division at Ames Research Center, and XSEDE TACC (TG-AST180008).
  
  
  
  § SIMULATION DETAILS
  
  In addition to the the benchmark M = 20 simulation discussed in Section <ref>, five simulations were performed to control for Mach number and box width. The parameters of each simulation are given in Table 1 and from left to right are: Simulation ID, Mach number (M), upstream ion plasma beta (β_i), simulation speed of light (c), grid size (Δ x), time step (Δ t), simulation length (L_x) and simulation width (L_y). Each simulation uses 16 macro-particles per ^2. 
  The reduction in the box length for  M80 was to compensate for the added computational cost of the reduced time step.
  
  
      
  .7
       
   Sim     M     β_i     c (0)     Δ x ()     Δ t () 
           L_x  (d_i)     L_y  (d_i) 
  
   M10     10     2     50     0.5     0.005     10^5     200
  
   M20     20     2     100     0.5     0.0025     10^5     200
  
   M20w     20     2     100     0.5     0.0025     10^5     1000
  
   M40     40     2     200     0.5     0.00125     10^5     200
  
   M80     80     2     400     0.5     0.000625     5× 10^4     200
  
    
  
  
  
  
  
  
  
  
  § THE SOLUTION TO THE CR-MODIFIED JUMP CONDITIONS
  
  The total compression ratio can be found by combining equations <ref>, <ref>, and <ref>. The equations can be rewritten into a single quartic equation in terms of X=√() where the coefficients depend on M_s, M_A and the post-shock pressures ξ_c and ξ_B, namely:
  
      c_1    X^4 + c_2 X^3 + c_3 X^2 + c_4 = 0;
  
  with the coefficients:
  
      c_1    = 1 + η_g/γ_g M_s^2 + 3/M_A^2
      
          c_2    = -√(2ξ_B)(η_cξ_c + 4ξ_B)
      
          c_3    = η_g(ξ_c + ξ_B - 1 - 1/γ_g M_s^2 - 1/2M_A^2) - η_cξ_c - 6ξ_B
      
          c_4    = η_g - 1.
  
  Since γ_g=5/3, η_g =5, while for γ_c between 5/3 and 4/3, η_c varies between 5 and 8;
  in our benchmark simulation, γ_c≃ 1.5 and η_c≃ 6 (see Appendix <ref>), which is also used for the predictions for different Mach number simulations in Figure <ref>.
  Note that it is not possible to introduce an effective adiabatic index for the wave/magnetic fields because the relationship between pressure and energy density is non trivial even when assuming Alfvén waves (Equation <ref>).
  
  Equation <ref> has 4 roots, but only 1 of which is physically relevant. 
  Since ξ_c,ξ_B≪ 1, two of the roots are negative and can be neglected. 
  Of the two positive remaining roots, only one corresponds to an increase in density, temperature and entropy at the shock and thus is the physical solution. 
  Note that for non-null ξ_c and ξ_B, the trivial solution X=1 disappears.
  
  
  
  § THE CR ADIABATIC INDEX
  
  To correctly predict the hydrodynamic modifications, an effective adiabatic index, γ_c, must be determined for the transrelativistic CR distribution.
  γ_c is measured as the ratio of the enthalpy density to the internal energy density:
  
      γ_c = 1 + ∫_p_inj^∞ vp^3/3f dp/∫_p_inj^∞ m_ic^2(Γ - 1)p^2f dp
  
  where Γ is the Lorentz factor.
  For non-relativistic particles, mc^2(Γ - 1) ≈ pv/2 and  γ_c = 5/3, while for relativistic distributions, mc^2(Γ - 1) ≈ pc yielding γ_c = 4/3. 
  In general, the exact value for the γ_c will depend on the CR distribution function, but will be bounded by these two values.
  The top panels of Figure <ref> show the upper and lower bounds for the solution for  based on Equation <ref>. 
  In both the non-relativistic (left) and relativistic (right) cases, we obtain ≳ 4 even for modest values of ξ_B and ξ_c; in general,  is larger in the relativistic case.
  To further illustrate the effect of the CR equation of state, the predicted compression ratio for a fixed value of ξ_c=0.1 and varying ξ_B are shown in the bottom panel of Figure <ref>;
  varying the CR adiabatic index modifies  by about 10% at most, an effect less important than the one induced by the CR drift in the precursor, which is controlled by ξ_B.
  
  
  
  
  