important to note that, besides only applying to high-energy particles with mean free paths larger than $L$, these equations are qualitatively different from those of [test-particle]{} DSA in an isolated shock: the spectrum $f_i$ below the exponential break is harder and the acceleration time is shorter.

Our model assumes a high level of CR-driven magnetic instabilities and Bohm diffusion for CRs in the close vicinity of the shocks [see, e.g., @bell04; @schureea12; @beov14]. However, once high-energy particles obtain mean free paths on the order or larger than the distance between the shocks, $L$, scattering will become much weaker. A specific feature of the simulation is that the highest energy CRs, with ${p_\mathrm{max}}\geq  p \gsim p_{\star}$, propagate with little scattering. Despite the weak scattering between the shocks, their momenta are still nearly isotropic since they scatter for long periods in the extended regions downstream from the shocks. Here, $p_{\star}$ is the momentum such that the proton mean free path $\Lambda(p_{\star}) \gsim  L$. Since the particle distributions are nearly isotropic even for high-energy particles with $\Lambda(p) \gsim  L$, the kinetic equation reduces to the so-called telegrapher equation [e.g., @earl74] and this allows a smooth transition between the diffusive and the scatter-free propagation 