Numerical Modeling of Time Dependent Diffusive Shock Acceleration

Motivated by cosmic ray (CR) re-acceleration at a potential Galactic Wind Termination Shock (GWTS), we present a numerical model for time-dependent Diffusive Shock Acceleration (DSA). We use the stochastic differential equation solver (DiffusionSDE) of the cosmic ray propagation framework CRPropa3.2 with two modifications: An importance sampling module is introduced to improve statistics at high energies in order to keep the simulation time short. An adaptive time step is implemented in the DiffusionSDE module. This ensures to efficiently meet constraints on the time and diffusion step, which is crucial to obtain the correct shock spectra. The time evolution of the spectrum at a one-dimensional planar shock is verified against the solution obtained by the grid-based solver VLUGR3 for both energy-independent and energy-dependent diffusion. We show that the injection of pre-accelerated particles can lead to a broken power law spectrum in momentum if the incoming spectrum of CRs is harder than the re-accelerated spectrum. If the injected spectrum is steeper, the shock spectrum dominates at all energies. We finally apply the developed model to the GWTS by considering a spherically symmetric shock, a spiral Galactic magnetic field, and anisotropic diffusion. The time-dependent spectrum at the shock is modeled as a basis for further studies.


Introduction
The cosmic-ray energy spectrum and composition is by now well-studied from GeV up to ZeV energies.The origin of these highly energetic, charged particles is still largely unknown.
Particularly unresolved is the question of the energy range between 10 15 eV and 3 × 10 18 eV.This energy range is defined by two kinks in the spectrum, the "knee" and the "ankle".Up to ∼ 10 15 eV it is believed that CRs are accelerated in the Galaxy, most likely at Supernova Remnants (SNRs) (see e.g.[1]).Above the ankle (∼ 10 18 eV) [1] CRs have gyroradii larger than the Galaxy and are clearly of extra-galactic origin.The breaks in the spectrum may indicate changes in the contributions to the spectrum: The spectral softening at the knee could be due to the maximal energy reached by Galactic sources.It could also result from a change in the energy-dependence of the residence time of CRs.The break at ≈ 200 − 300 GV especially at lower energies ( [1] and references therein) are also attributed to a change in the transport properties.
In general, the observed power-law spectrum indicates that CRs are accelerated by stochastic processes.There are two possible scenarios of how to accelerate particles stochastically, both developed early-on by Fermi.Fermi Second-order Acceleration [2] was suggested first as acceleration on moving magnetized clouds in the Galaxy.Fermi First-order Acceleration [3][4][5][6][7] is more efficient: CRs are scattered on both sides of a shock and may cross the shock front repeatedly.Each time they cross the shock, they are accelerated depending on their energy, ∆E ∝ E. The energy gain can be derived by Lorentz boosts between the upstream and downstream reference frames.For non-relativistic, stationary shocks and neglecting non-linear interaction between the CRs and the background plasma the slope of the stationary shock spectrum is purely dependent on the shock's compression ratio.Even for non-linear diffusion coefficients the spectral form is unchanged [8].
There are some limitations to DSA, one of which is that particles must be able to reach the shock from downstream to upstream again, see, e.g.[9].In order to cross the shock, CRs must have been accelerated to energies that give them a gyroradius larger than the shock width, which is also known as the injection problem.Also, CRs have to be contained in the acceleration region to experience ample acceleration.This implies a diffusive particle transport with diffusion low enough that CRs are not escaping too quickly.Containment can also be implied by the geometry of the system, e.g. in stellar clusters particles can only escape downstream.
An upper bound of the maximum energy produced by a source can be estimated by the Hillas criterion and depends on the particles' rigidity and the size of the acceleration region [10].Often magnetic field amplification due to CRs determines the maximum energy reached by DSA at the shock.There are other limiting factors like the lifetime of the accelerator or loss processes like synchrotron radiation.
In the transition region from Galactic to extragalactic CRs the Galactic Wind Termination Shock (GWTS) might contribute to the sources of CRs via re-acceleration [11,12].A possible Galactic wind can be driven by, e.g., radiation or the CRs themselves (e.g.[13][14][15][16][17]).It may form a shock when the supersonic wind slows down due to interaction with the Intergalactic Medium (IGM).Especially for starburst galaxies there is evidence for supersonic outflowing winds (e.g.[18,19]).
CRs with their origin in the Galactic disk are advected outwards and are accelerated to higher energies at the GWTS.Either they leave the Galaxy or a fraction of the high-energy CRs may be able to propagate back to the Galaxy against the Galactic wind.The idea of particle acceleration at the GWTS was already discussed e.g. by [20], [11], [12], and most recently by [21].With the assumption that about two solar masses per year are advected outwards, the GWTS cannot be supported by the wind longer than 100 Myr [22].The finite lifetime of the shock may also impact the shock spectrum and its contribution to the observed cosmic-ray spectrum on Earth.
Not only the origin of CR but also their transport properties change in the transition region.Particles undergo a random walk due to scattering on magnetic field turbulence.The interaction of CRs with the magnetic field turbulence depends on the gyroradius and the turbulence spectrum of the magnetic field.CRs with energies above ∼ 10 18 eV engage less in interaction with the turbulent Galactic magnetic field and can only be described by the diffusion approximation at late times, and with that, on larger scales [23].The transition time from ballistic to diffusive transport generally grows with energy.In energy, it is still unclear where exactly the transitions between diffusive and ballistic transport as well as Galactic and extra-galactic CRs lie.
Methods for the simulations of CR transport have to reflect that transition: For highenergy extra-galactic CRs the equation of motion is integrated.In this ballistic regime, individual particle trajectories in arbitrary magnetic field configurations can be simulated precisely.However, for Galactic CRs with energies below 10 16 eV this approach becomes computationally costly as most of the computation is spend on resolving the particles' gyration.Therefore, for Galactic CRs the transport equation is commonly used to describe diffusive transport in space and momentum.
The simplified, assuming isotropic (in momentum space) CR distributions, Fokker-Planck equation, often called Parker transport equation (e.g.[24]) describes the time evolution of the cosmic-ray differential number density N = p 2 f (⃗ x, p, t) in space ⃗ x and momentum p = |⃗ p|, with f being the CR distribution function.Individual terms describe advection of CR with a background flow ⃗ u, spatial diffusion described by the diffusion tensor κ, momentum diffusion described by the coefficient D, adiabatic energy changes and cosmic-ray sources S.
In the test-particle picture, the average change in the particle's momentum can be calculated from changing the local fluid frames when the shock is crossed.Together with the probability to escape downstream and never return to the shock, the well-known power-law in momentum can be derived (e.g.[5] or microscopic derivation in [25]).In the diffusive picture, the ensemble-averaged distribution function f or density N , are described by the transport equation in the frame in which the shock front is stationary.Here, diffusive shock acceleration results from adiabatic heating and the interplay between advection and spatial diffusion (e.g.[4], [26] or macroscopic approach in [25]).
The spatial diffusion tensor is often defined in a local coordinate frame, describing diffusion parallel and perpendicular to the local magnetic field.The momentum diffusion coefficient describes acceleration of CR in the presence of magnetic field turbulence (Second-order Fermi acceleration).The knowledge of both, spatial and momentum diffusion tensors, is essential for the complete description of CR transport.
There is not yet a closed description of the diffusion tensor.Depending on the assumed magnetic field, coherent background and turbulent component, and the CRs properties analytical or numerical approximations might be available (e.g.[27][28][29][30][31][32][33]).On the other, hand it can also be used as a free parameter and fitted to match, e.g., the observed primary and secondary ratios in a model of Galactic cosmic-ray transport.Since all of these approaches come with their own caveats, we will use normalized diffusion coefficients κ∥ = 1 throughout this work.
Considering the GWTS mentioned at the beginning, re-accelerated CRs are expected to be in the transition from diffusive to ballistic transport.We expect a finite lifetime for the GWTS, so that the transition time between the transport regimes is important.Close to the shock due to strong turbulence induced by the CRs themselves the diffusive approach is valid for higher energies than in the interstellar or IGM at a fixed time [5,34].Particles escaping the acceleration region, however, cannot necessarily be described by the diffusion approximation depending on their energy.In order to assess the contribution of the GWTS to the cosmic-ray spectrum they must be propagated back to the Galaxy using either the ballistic or diffusive approach depending on their energy.Thus, the spectrum observed on Earth is further modulated by the energy-dependence and spatial variability of the diffusive transport.
The CR propagation framework CRPropa3.2[35] offers the possibility to simulate acceleration and propagation in both regimes.Also, the structure of CRPropa3.2,which is described in more depth in Section 2, makes it relatively easy to define arbitrary magnetic field and shock configurations in up to three dimensions.Together with the GWTS scenario, this motivates an in-depth study of DSA using the diffusion approach of CRPropa3.2[36].
In test scenarios we explore different questions related to the re-acceleration of CR at the GWTS: How does the spectrum change when a finite shock lifetime is considered?What are the effects on the spectrum when CRs are already pre-accelerated to a power-law?To what extent does a finite acceleration region affect the maximal energy that can be reached?
First, we give a short overview of CRPropa3.2and its diffusion approach which is based on Stochastic Differential Equations (SDEs) in Section 2. A new module is presented that was specifically implemented to enhance statistics when simulating DSA with CRPropa3.2.
In Section 3 we show that DSA simulations with CRPropa3.2result in the expected spectral slope for acceleration at one-dimensional planar shocks.We validate the obtained stationary spectra with predictions from theory and other ensemble-averaged approaches to DSA.The time evolution of the spectrum at the shock is compared to simulations that integrate the transport equation using the finite difference code VLUGR3 [8,37,38].We clarify constraints for simulating DSA based on SDEs.Based on those findings, we consider a finite shock lifetime, energy-dependent diffusion and the injection of pre-accelerated spectra.The effects of each modification on the spectrum at the shock or the downstream boundary are analysed separately.
In Section 4 we take anisotropic diffusion into account and consider oblique shocks.Finally, we present the time-dependent spectrum at a simple model for the GWTS: A spherical symmetric shock and spiral magnetic field.

Methods
In this Section we briefly introduce SDEs and how they relate to the transport equation (1.1).We explain how they are used in CRPropa3.2 to solve the transport equation.For the simulation of DSA a new module, CandidateSplitting, was implemented to enhance statistics at high energies.

Stochastic Differential Equations
One way to simulate DSA is the solution of the transport equation (1.1), a partial differential equation in reduced phase-space.The solution directly results in the distribution function of CRs.We followed this approach using the well-established VLUGR3 [37,38] code to cross check our results.Another way is to make use of the connection between Fokker-Planck type equations and SDEs.In general, Fokker-Planck type equations describe the transition probability from one state to a set of other states.SDEs describe dynamical systems which are subject to noise.
In the following we discuss briefly the set of SDEs implemented in CRPropa3.2.For an in-depth explanation the reader is referred to [36] or e.g.[39].
In the absence of second-order Fermi acceleration processes, i.e.D = 0, and spatially constant diffusion the system of SDEs equivalent to the transport equation (1.1) is given by ) where d⃗ ω t = √ dt ⃗ η is a Wiener process with η i being random numbers from a unit normal distribution.The tensor B corresponds to spatial diffusion.In the local frame of the magnetic field line the diffusion tensor becomes diagonal.Drift terms due to curvature of the magnetic field lines are described by off-diagonal elements and are neglected in the following.With this assumption B is given by B ij = δ ij 2κ ij .Equation (2.1) describes the spatial displacement due to advection and stochastic fluctuations.Equation (2.2) describes the adiabatic energy change.Without diffusion in momentum space it is not an SDE but an ordinary differential equation.
SDEs can also be written in integral form by making use of the Itô integral [40].In this form they can be approximated numerically (see e.g.[41,42]).The SDE approach has the advantage that it is easy to implement and to extend to arbitrary geometries.We describe the first-order Euler-Maruyama scheme implemented in CRPropa3.2for numerical solution of the SDE (Eq.(2.1)) in the following Section.

Simulation with CRPropa3.2
CRPropa3.2 has a modular structure with the Candidate class as central element.When the EoM is used, the candidate simply corresponds to the propagated particle.When the diffusive transport is used, the candidate corresponds to a phase-space element that is propagated.The phase-space element is also called pseudo-particle in the following.The candidate module holds all information about the (pseudo-)particle that can be used and altered by other modules in the ModuleList that are successively called each simulation step.The user can add different kinds of modules to the ModuleList: Propagation, Acceleration, Interaction, Boundaries and Observers.Other modules handle the simulation environment (MagneticField, AdvectionField), the injection of candidates (Source) and the output.The modular structure makes it easy to set up simulations for various scenarios and to add new modules.
To model DSA, we use the DiffusionSDE module [36].Pseudo-particles are propagated parallel and perpendicular to a magnetic field, depending on the diffusion tensor which is constant in space and can depend on energy.The diffusion coefficients for the normal and binormal direction of the magnetic background field line can be assumed to be the same for large curvature radii of the magnetic field lines κ n = κ b ≡ κ ⊥ = ϵκ ∥ .In Section 3.4 energy-dependent diffusion is considered.Other than that, diffusion is considered to be energy-independent.The SDE 2.1 is integrated using the Euler-Maruyama Scheme (e.g.[41]) with time step ∆t and random variables ∆⃗ ω n = ⃗ ω t n+1 − ⃗ ω tn drawn from a normal distribution with an expected value of zero and variance of ∆t.SDE methods have no constraints on the used time step.However, the choice of time step can be crucial to obtain correct results.See [43] for a discussion.
Eq. (2.4) in general is defined in the lab frame.Since the diffusion coefficient in Eq. (2.3) is defined in the local magnetic field frame, the diffusive step √ 2κ • ∆⃗ ω n is calculated in the orthonormal base of the magnetic field by integrating along the magnetic field line for parallel diffusion, calculating the perpendicular diffusive step and transforming back to the lab frame in each time step.
To ensure one-dimensional propagation in Section 3 the magnetic field is set constant in x-direction and diffusion is allowed only parallel to the magnetic field, i.e. ϵ = 0.In Section 4 perpendicular diffusion is taken into account as well as a more complex magnetic field.
Acceleration is handled by the AdiabaticCooling module according to the momentum equation (Eq.(2.2)) based on the divergence of the shock profile specified in the Advection-Field class.The advection field and magnetic fields used for simulations are specified in the following sections.
Free-escape boundaries L ± can be specified in the upstream and downstream region.Observer modules detect pseudo-particles when they cross those boundaries or the TimeEvo-lutionObserver is used to detect pseudo-particles' positions and energies at given times t i during the simulation in the acceleration region.

Analysis
The energy spectrum J(E) = dN/dE of particles in the acceleration region is approximated by the histogram ∆N/∆E, with ∆N being the number of pseudo-particles in each energy bin ∆E.In the high-relativistic limit, assuming E = p/c, J ∝ p 2 f (p, t).Since we expect a power-law spectrum, energy is binned with equal distance in logarithmic space.The error for each bin is then given by ∆J = J/ √ ∆N , any additional errors (e.g. from summing over time "snaps" as explained in Section 3.1) are not included.The spectrum can be weighted by E 2 to highlight the spectral slope s = −2 predicted for acceleration at strong shocks.

Candidate Splitting
With the assumptions of non-relativistic shocks and neglecting non-linear interactions, DSA produces a power-law spectrum in energy and the spectral slope depends only on the shock compression ratio.When mono-energetic candidates are injected at the shock and accelerated until they leave the acceleration region, the number of candidates decreases with increasing energy.Especially for less efficient acceleration at shocks with a low compression ratio, statistics at high energies may be so bad that it is difficult to evaluate whether the stationary solution is already reached or to determine the resulting spectral slope.
The number of injected candidates can be increased in order to get sufficient statistics at high energies, which comes at large computational costs.Another way to overcome this problem is to split candidates in n split copies once they cross energy boundaries in log-space.For DSA the optimal splitting number depends on the expected spectral slope to compensate for the loss of candidates at higher energies.Each splitted candidate is assigned with a weight w = 1/n split during the simulation.In the later analysis candidates are weighted accordingly to obtain the correct spectra.In order to determine the error of the spectrum for each energy bin all candidates can be used, which reduces the uncertainty at high energies.
In Fig. 1 resulting shock spectra simulated with and without candidate splitting are compared at the same simulation time (t = 400 t 0 ).The injected number of candidates is N 0 = 10 6 for the simulation without and N 0 = 10 4 with candidate splitting.For the latter, the  Left: Spectra obtained with and without using the CandidateSplitting module at simulation time t = 400 t 0 evaluated at the shock, x = [0, 2] x 0 .The shock compression ratio is q = 4 and both simulations reproduce the expected spectral slope s = −2.The simulation with candidate splitting reaches higher energies up to 10 7 E 0 with smaller errors.Right: Weighted spectra (top) and relative error (J −2 − J)/J where J −2 is the predicted spectrum (bottom).
actual number of simulated candidates increases during the simulation.Still, computational time goes down by a factor of 251 when candidate splitting is used, since fewer candidates are simulated in low-energy bins.For the set-up shown in Fig. 1, at t = 400 t 0 , the number of candidates is about 2.1N 0 .
Note that the depicted errors are given by ∆J = J/ √ ∆N .To approximate the stationary solution of the energy spectrum at the shock, energy of pseudo-particles stored during simulation are summed over time as explained in Section 3.1.Such time related errors are not shown in Fig. 1 and following figures.
Splitting of candidates at energy bins in logarithmic space is one way of importance sampling used also in other approaches to DSA [44].Another way is the introduction of artificial drift terms pushing particles to the shock [43,45].Such drift terms need to be extracted in the analysis.The candidate splitting method is, however, easier to apply to more complex configurations of the magnetic field and advection field.
Candidate splitting is implemented as an independent module in CRPropa3.2.It is possible to define individual energy bins and number of splits.For use in diffusive shock acceleration only the minimal and maximal energy as well as the expected spectral index need to be specified.

One-dimensional Diffusive Shock Acceleration
In the ensemble averaged picture, energy gain at shocks is described by Eq. (2.2).To apply this adiabatic description a finite shock width is considered since the advective speed u(x) must be continuously differentiable.Deviations from a sharp shock transition can occur if a population of energetic particles is interacting with the thermal background plasma.In One-dimensional advection field with shock at x = 0 and a shock compression ratio q = 4.The smaller the shock width compared to the integration step length, the better the ideal shock solution is approximated.such case the latter can experience a decelerating force due to the pressure gradient resulting from the (energy density) distribution of the suprathermal particle population and either a so-called subshock forms or the shock weakens into a smooth transition (see, e.g., [46][47][48][49][50][51]). The one-dimensional shock profile is implemented with upstream and downstream velocity, u 1 and u 2 = u 1 /q, and shock width L sh .In the following, all units are normalized so that with x 0 /v 0 = t 0 and E 0 being the energy injected at the shock.Thus, different physical scenarios can be modeled with the same simulation by adjusting the time scale.For the CRPropa3.2simulations we used x 0 = 1 km and v 0 = 1 m/s, so that t 0 = 1000 s and κ 0 = x 0 v 0 = 1000 m 2 /s.The velocity profile u(x) is also used in other studies of linear and non-linear DSA [8,44,52].Figure 2 illustrates u(x) for a shock compression ratio q = u 1 /u 2 = 4 and different shock widths.The narrower the shock width, the better the approximation of the ideal shock.However, from the numerical perspective the region where ∂u/∂x ̸ = 0 must be large enough compared to the integration step length for candidates to experience acceleration.Otherwise energy gains may be underestimated (see Section 3.2 for a detailed analysis).

Stationary Solution
Simulating with CRPropa 3.2, pseudo-particles are injected at t = 0 at the shock and propagated until a maximum simulation time or until they leave the acceleration region downstream through free-escape boundaries.
In order to obtain a stationary solution to compare the results with shock spectra predicted from theory, s = 2 − 3q/(q − 1) (e.g.[25]), continuous injection upstream must be assumed.The inherent structure of CRPropa 3.2 does not allow for continuous injection during the simulation itself.To construct the stationary solution at the shock afterwards CR-Propa's TimeEvolutionObserver module is used.Positions and energies of pseudo-particles are stored at times t 0 + n∆T , which are called time snaps and do not necessarily have to be equal to the simulation time step ∆t.With the assumption of continuous injection, the solution at time t is constructed by summing over all time snaps to approximate the timeintegrated solution.This approach is valid given that the solution does not change too much during one time interval ∆T .The time intervals ∆T can be chosen linearly or to increase with time, when the solution undergoes smaller changes.The resulting time evolution of the spectrum at the shock and number density integrated over momentum are shown in Fig. 3 compared to the solutions obtained by the integration of the Fokker-Planck equation (1.1) using the same parameters.For the integration of the partial differential equations VLUGR3 [37,38] is used.The spectra are weighted by energy so that the predicted slope lies horizontally in the plot.The resulting spectra and number densities are in good agreement.Differences may result from binning in energy and space for the CRPropa3.2simulation and from resolution in energy and time for the VLUGR3 simulation.The time evolution can be interpreted in the following ways: The shock is active from time t = 0 on and accelerated particles entering the shock region with energy E 0 .The longer the shock is active the higher the maximum energy the particles can reach and the cut-off of the spectrum moves to higher and higher energies.Having an infinite acceleration region, the stationary solution of a planar one-dimensional shock with compression q = 4 is a power-law with slope s = −2.With increasing simulation time, the approximation of the stationary solution gets better and the solution relaxes to the power-law for ever higher energies.In Fig. 3 the stationary solution at t = 200 is reached for energies E < 10 3 E 0 .Theoretically there is no limit for the maximal energy, due to the infinitely large shock front every n-th crossing there will always be a nonzero number of particles to cross the shock an n + 1-th time for t → ∞.Considering a spherical geometry as discussed in Section 4.3, the maximal energy that can be reached depends on the diffusion coefficient and shock radius.
In order to evaluate the approximation of the stationary solution the spectrum is fitted at t = 500 (not shown in Fig. 3) for energies below 10 4 E 0 .Up to that energy candidates are splitted by the CandidateSplitting module.The resulting spectral slope of −2.054 ± 0.003 matches the expected value s = −2 nicely.
The time evolution of the number density shows how particles gradually migrate into the downstream region due to the advection field.At t = 200 the stationary solution is reached close to the shock, x < 30, at later times also further into the downstream region.The freeescape boundary at x = 100 does not influence the number density profiles at the shock.In principle, free escape boundaries are only needed if the acceleration region has a finite size.In Fig. 18 we show how the resulting spectrum is affected by free-escape boundaries.Only a fraction of particles makes it into the upstream region x < 0 against the advection flow.Acceleration is efficient for a constant, low diffusion coefficient as assumed here, however, particles are unlikely to escape upstream.
The analytical solution by [53] is only for the special case of v 2 /κ = const.across the shock, resulting in a drop of spatial diffusion by 1/16 for a strong shock.Forman and Drury [54] give an approximate solution for momentum dependent κ, which is exact if v 2 /κ = const..Such spatial dependence of the diffusion coefficient is not considered here, since the resulting drift term which adds to the advective step in Eq. (2.1) is currently not implemented in CRPropa3.2.For SDEs the strong additional drift term that is expected at the shock may also lead to numerical difficulties [44].The integration of the transport equation using VLUGR3 is verified against existing analytical solutions in earlier works [8] and by comparisons of the acceleration time scale [25,54].In Appendix D the approximate time-dependent solution from [54] is compared to the solution obtained with VLUGR3.

Constraints
In order to reproduce the predicted spectrum the simulation set-up must fulfill constraints on the integration step length and shock width.For SDEs Achterberg and Krülls [52] and later Achterberg and Schure [44] present a thorough analysis of the choice of shock width and step length and the resulting spectral slope.The simulations resulting in Fig. 3 take their findings into account.In the following we show that DSA simulations with CRPropa3.2are subject to the same constraints.These constraints are not purely numerical but can also be motivated from physics.
Considering a diffusion coefficient constant in space and energy, the advective step length, diffusive step length and shock width must fulfill the inequality to obtain the correct spectral index [52].
The first inequality ensures that pseudo-particles experience the gradient of the advection field and therefore gain energy.Numerically, the time step must be chosen sufficiently small to resolve the shock region.On the other hand, the diffusive step length √ κ∆t -a measure for the stochastic step -must be larger than the shock width to increase the likelihood of pseudo-particles to cross the shock front multiple times in finite simulation time.Depending on the chosen shock width, diffusion coefficient and advection field it may not be possible to fulfill the inequality.In that case, acceleration can be under-or even overestimated.In Appendix A the resulting spectral slopes for various shock widths and time steps are shown, with the physical parameters, advective speed and diffusion coefficient, being held constant.
Constraints on the time step come from the numerical method, however, not all combinations of u, L sh and κ can be physically motivated.For instance, if diffusion is low compared to the advective speed there will be negligible acceleration.
Obeying Eq. (3.3) alone is not sufficient to reproduce the ideal shock spectra.If the shock width is too large, pseudo-particles experience a smoothly changing advection field instead of a discrete shock.For diffusion independent on space and energy Krülls and Achterberg [52] calculated the expected spectral slopes for such smooth velocity gradients.In Fig. 17 we show analogous to [52] that the expected spectral slopes for such smooth shock waves are met if the time step is sufficiently small.Also, when a finite acceleration region is considered, diffusion must not only be large enough so that particles are able to diffusive back to the shock but also low enough so that they are contained in the acceleration region and do not escape too quickly.The effect on a finite acceleration region on the resulting spectral index is shown in Fig. 18.
In the following we show the effect of the chosen time step and shock width as well as the diffusion coefficient.We aim to reproduce the predicted spectral slope s = 2 − 3q/(q − 1) for compression ratios q = 2 and q = 4.

Advective Step
Acceleration at shocks of different widths sh is simulated with decreasing time steps ∆t.The normalized diffusion coefficient is κ = 1 and the upstream velocity ũ1 = 1.The shocks compression ratio is q = 2 leading to a spectral slope s = −4 in order to compare with the predictions from [52].The slope is determined by linear fits in log-space up to the maximum energy of the CandidateSplitting module.Time is chosen such that the cut-off of the spectrum has minimal impact on the slope approximating the stationary solution.The results are shown in Fig. 4.
The expected spectral slopes depending on the shock widths are indicated by lines.With increasing resolution the simulated spectrum well approximates the predicted stationary solution.For sufficiently small shock widths, here Lsh = 0.008, a spectral slope close to that of an ideal shock is obtained.The larger the shock width, the less efficient the acceleration which leads to steeper spectra.Depending on the shock width, with decreasing time step, the simulation approaches the predicted stationary solution from steeper ( Lsh = 0.128) or flatter ( Lsh = 0.008) spectra.Thus, the compression ratio of the shock is over-or underestimated.
We show a more detailed figure with entangled information about the shock width and time step in the Appendix in an analogous way to [52].

Diffusive Step
A similar analysis can be done for the dependence of the diffusion coefficient in relation to the shock width.The larger the diffusive step, the better becomes the approximation of an ideal shock.If the diffusive step however is too large, pseudo-particles may miss the shock region and acceleration is underestimated as a consequence.
For a constant shock width of Lsh = 0.004 the diffusion coefficient was varied over several simulation runs.The resulting spectral slope was fitted in the energy range [E 0 , 10 3 E 0 ] given by the CandidateSplitting.The results are shown in Fig. 4 in comparison to the results from [44].With large diffusive step length the spectrum may also be affected by free-escape boundaries.For the simulations boundaries were chosen such that they do not influence the time evolution of the spectrum for the simulated diffusion coefficients.Likely, that is why we do not find the spectral slope to decrease again for ϵ < 0.04.In the Appendix we show how the free-escape boundaries influence the spectrum at the shock depending on the diffusion coefficient.The compression ratio is q = 2, thus for an ideal shock we expect s = −4.With larger shock widths, the shock produces steeper spectra as calculated by [52].Slopes of the predicted stationary solutions are indicated by the colored lines for the respective shock width.Right: Resulting spectral slope depending on the ratio of upstream velocity and shock width to diffusion coefficient.Data from [44] is shown for comparison.Here ũ1 = 1, q = 4 and Lsh = 0.004 with ∆ t = 0.001 is used while κ is varied.The free escape boundaries are L± = 10.With higher diffusion coefficient the spectral slope approaches the predicted value s = −2.We already show the results for a transversal shock wave and taking perpendicular diffusion into account, which is discussed in Section 4.1.
We conclude that for the simulation of ideal shocks, the shock width must be small compared to the advective field and the time step should be chosen such that the advective step is at most as large as one fourth of the shock width.The diffusive step must be larger than the shock width, but not too large to still contain particles in the acceleration region.Those constraints may be difficult to fulfill when modeling a physical scenario and may require very small time steps close to the shock.

Finite Lifetime
We model two scenarios that lead to a shock spectrum that is not stationary: Either the source of particles reaching the shock is only active for a short period of time or the shock itself has a finite lifetime.Considering the GWTS, the CR flux that reaches the shock from the Galactic disk can be assumed to be stationary.The GWTS itself however cannot be maintained for an infinite time since too much mass would be advected out of the Galaxy.
To model a burst-like particle source, candidates are injected at t = 0 and propagated until they leave the acceleration region through the free-escape boundaries L ± .The time evolution depending on the source duration can be derived similar to the stationary spectrum in Section 3.1.Energy and position of candidates are stored in time intervals ∆T and the resulting spectrum is approximated by summing of those time snapshots from [t − ∆T src , t] with ∆T src being the source lifetime.For more details we point to [12]. Figure 5 shows the time dependent spectrum at the shock for two different source duration ∆T src .Considering a finite lifetime of the shock, we already know the time dependent spectrum at the shock from Section 3. In order to assess the contribution to the overall cosmic-ray spectrum, not only the spectrum at the shock, but also the spectrum of particles escaping downstream, or even upstream is of interest.Fig. 6 (left) shows the time evolution of the downstream spectrum, x = 50, assuming continuous injection at the shock.Compared to Fig. 3 the spectrum evolves later, since particles need time to get from the shock to the downstream position.Considering now a free-escape boundary at L+ = 100, Fig. 6 (right) shows the timedependent spectrum of particles escaping with a lifetimes of the shock ∆ Tsh = 100.For the first time steps, the spectrum still evolves up to the quasi-stationary solution ( t = 450).Particles that stay long in the acceleration region, defined by the free escape boundaries L ± , reach higher energies than particles that escape quickly.Thus, at late times ( t = 620, 860), the spectrum gets flatter since low energy particles already escaped.
With a constant, low diffusion coefficient κ = 1 acceleration is efficient but particles do not escape upstream.In the following Section we introduce energy-dependent diffusion: With increasing energy, diffusion is higher and particles are more likely to escape upstream.A finite lifetime ∆ Tscr = 100 of the shock is assumed.The spectrum of particles escaping through the downstream boundary is shown.

Energy-dependent Diffusion
So far, a diffusion coefficient constant in energy was considered.A more realistic description of diffusive motion takes energy-dependence into account.Assuming that the diffusion coefficient is proportional to (E/E 0 ) α and α > 0, the diffusive step becomes larger with increasing energy.At first glance it is then easy to satisfy the inequality in Eq. (3.3).However, the analysis in Section 3.2 revealed that the diffusive step should not become too large and free-escape boundaries should be dropped or must be set far away from the shock, otherwise high-energy particles escape too quickly.
In order to keep the diffusive time step within a reasonable range, an energy-dependent adaptive time step is implemented in the DiffusionSDE module.Within a range of specified time steps [dt min , dt max ], it is chosen such that the advective step is smaller than 1/4 L sh and the diffusive time step smaller than 100 L sh .This also leads to better performance in case of the energy independent diffusion coefficient (α = 0), since the time step in the downstream region can be chosen higher than in the upstream region by a factor of the shock compression ratio q.
Figure 7 shows the time evolution of the spectrum at the shock and number density in the acceleration region.Again the compression ratio is q = 4 and the predicted spectral slope s = −2.The diffusion coefficient is κ(E) = κ E 0 (E/E 0 ) α with α = 1 and κE 0 = 1.The results obtain with the diffusion approach applied in CRPropa3.2are compared to the simulations using VLUGR3.The same times as in Fig. 3 are used.
With α > 0 the diffusion coefficient κ(E) is equal or greater than in the energyindependent case shown in Section 3.1.With that the acceleration time scale gets larger.The cut-off of the energy spectrum at the shock is at lower energies compared to Fig. 3.  an energy distribution of a power law [11].Depending on the acceleration process a range of spectral slopes is possible.
We injected different spectra at the shock and investigated the impact on the resulting re-accelerated spectra.According to [25] steeper spectra than the one produced by the shock converge to the shock spectrum.When, on the other hand, the shock produces a steeper spectrum than the one continuously injected, the injected spectrum does not change.This results in broken power-law spectra when the injected spectrum has a finite cut-off energy as shown in Fig. 11.
In Fig. 10a a flatter spectrum than produced by the shock is continuously injected up to the maximum energy 10 2 E 0 .Up to that energy, the injected spectrum prevails.At the shock, particles are accelerated independent on their energy and the time evolution of the spectrum at the shock becomes visible for E > 10 2 E 0 .In Fig. 10b a steeper spectrum than produced by the shock is continuously injected up to the maximum energy 10 2 E 0 .The spectrum converges to the one produced by the shock over time, since all particles independent on their energy experience acceleration.Thus, at the shock particles are re-accelerated to higher energies.
In case of flatter injected spectra, broken power-law spectra may emerge.In case of steeper spectra, the spectrum converges to the shock spectrum.If the injected spectrum is the same as produced by the shock it remains the same.This is summarized in Fig. 11: The stationary spectrum up to 10 4 E 0 at the shock with different injected spectra from 10 −1 E 0 − 10 2 E 0 is shown.It is more likely that the injected spectrum has an exponential to decrease with r −2 downstream, r > r sh , modeled by Analogous to Eq. (3.1), a finite shock width L sh is considered.Similar profiles were also used in other studies of spherical DSA [55] and applied to the GWTS [12,20].After acceleration, the wind expands with constant speed u(r).Particles experience adiabatic cooling upstream, r < R sh .Downstream, r > R sh the wind is subsonic and decreases with u(r) ∝ r −2 without cooling.The advection field is shown in Fig. 13 for different shock widths.Figure 13: Radial shock profile.The advective field is constant upstream u(r) = u 1 , drops by 1/q at the shock and decreases downstream with r −2 for r > R sh .
For the magnetic field, analogous to [12], an Archimedean spiral is considered.It is given by with the constant wind velocity v w and Ωr 0 sin θ being the rotational velocity at a reference radius r 0 at latitude θ.The magnetic field lines for different values of Ω are illustrated in Fig. 14.For Ω = 0 the magnetic field is parallel to the shock normal and the advection field.For Ω ̸ = 0 the magnetic field is not parallel to the wind profile and the break of the magnetic field at the shock is neglected.For high values of Ω the magnetic field is almost perpendicular to the shock normal and resembles the transversal shock wave investigated in Section 4.1.In Fig. 15 spectrum and number density are shown for Ω = 10.Here, perpendicular diffusion was also taken into account to achieve sufficient acceleration.In both simulations the stationary shock spectrum matches the predicted spectral slope s = −2 for a strong shock.Due to the decelerating flow for r > R sh the number density profile differs from the planar scenario.
Figure 16 shows the time evolution of the spectrum at the shock and the number density in the acceleration region for Ω = 1.Different values for the perpendicular diffusion are spectral slope, shock width and advective step, as well as diffusive step.Likewise we found that depending on the shock width for too small time steps, the compression ratio of the shock is over-or underestimated (see Fig. 4 and Fig. 17).Presumably, this is due to the pseudo-particles only meeting the shock region by chance: The smaller the shock width the stronger the velocity gradient they encounter.This leads to harder spectra at least for those pseudo-particles experiencing acceleration.Section 3.2 explains in detail how to set up DSA simulations using SDEs in general and CRPropa3.2specifically.
Based on the findings in Section 3.2 for DSA an alternative adaptive time step was implemented in the DiffusionSDE module.This ensures that within a user specified range the largest possible time step based on the inequalities (Eq.(3.3)) for advective and diffusive step is used.Another extension to CRPropa3.2 is the CandidateSplitting module.High energies are only reached by a small fraction of candidates.For better statistics at high energies, candidates crossing energy boundaries are split into n copies depending on the expected spectral slope.This way of importance sampling significantly reduces computation time.
We approached more physical scenarios by considering energy-dependent diffusion, preaccelerated spectra and shocks with finite lifetime.Depending on the life-time of the shock, the spectrum seen by a downstream observer differs from the stationary spectrum.Energy or spatial dependence of the diffusion coefficient would further modulate the downstream spectrum over time.The effects on the shock spectrum discussed in sections 3.3 to 3.6 were deliberately considered separately and are the basis for future combined studies, potentially adding spatial dependent diffusion.
In Section 4 we took anisotropic diffusion into account and showed that planar oblique shocks produce decent results, when diffusion parallel to the shock normal is high enough.A prototype for simulating the GWTS was presented in Section 4.3: A radial shock profile and a Archimedean spiral magnetic field.The time-dependent spectrum obtained at the GWTS as well as at the upstream and downstream escape boundaries can be used in further studies, investigating the portion of particles that is able to propagate back to the Galaxy (as done for a stationary shock spectrum by [12]).Only high-energy CRs are able to propagate against the outflowing wind, so that energy-dependent diffusion must be considered.We already saw in Section 3.4 that energy-dependent diffusion leads to a higher fraction of particles being able to diffusive back into the upstream region against the wind.
For particles leaving the acceleration region the diffusion approximation might not be valid, depending on the escape time, magnetic field of the IGM and the particles' energy.At this point we stress again the advantage of CRPropa3.2:Follow up simulations can be done using both diffusive and ballistic propagation depending on the candidates energy within the same framework.
Not only energy-dependent diffusion but also a spatial dependence might impact the acceleration process as well as the cosmic-ray transport in the Galactic magnetic field.A spatially varying diffusion tensor, however, induces an additional drift term to Eq. (2.1) which is not yet implemented in CRPropa3.2.Also, this term adds to the left-hand side of the inequality in Eq. (3.3) discussed in Section 3.2.This makes it more difficult to choose a valid time step for a given shock width.Achterberg and Schure [44] even found that a second-order scheme might be necessary for handling the strong drift terms that would occur at the shock when the diffusion coefficient drops over the shock.
The GWTS is most likely not a perfect sphere but like the termination shock of the heliosphere a spheroid.Also, it does not necessarily enclose the complete Galaxy.With a more variable wind profile CRs may even propagate back to the Galaxy more easily.A more realistic wind profile presumably changes the time evolution of the spectrum produced by the shock itself.When a finite shock lifetime is considered this may have an impact on the contribution to the shock spectrum at Earth.
With the radial wind profile and Archimedean spiral magnetic field we modeled a quasi three-dimensional shock configuration.For the Galaxy, there exist more realistic magnetic field models, like [57,58] or [59].The magnetic field structure certainly has a great impact on the possibility of CRs travelling back to the Galaxy.Merten et al. [12] showed that the length of the magnetic field lines and the amount of diffusion parallel to the magnetic field lines is critical for a contribution of CRs re-accelerated at the GWTS to the spectrum.The magnetic field also has an impact on the arrival direction of CRs.With CRPropa3.2 the magnetic field can easily be exchanged and the impact on the spectrum and arrival direction on Earth can be studied for various magnetic fields.

A Advective Step
In Section 3.2 the constraints of the SDE approach are discussed.Figure 17 is a detailed version of Fig. 4 where shock width and time step are entangled.Analogously to [52], for different simulation time steps ∆t (or advective steps ∆x adv,1 , with ũ1 = 1) the shock width L sh is varied.Going along the dotted lines, the shock width increases from left to right.Depending on the shock width, the expected spectral slope according to [52] is indicated.In general, the smaller the time step, the better the predicted spectral slope is met.

B Diffusive Step
One advantage of the SDE approach is that there is no need for boundary conditions.However, free escape boundaries can impact the acceleration process and with that the spectral slope.Free escape boundaries are also set in the work of [44]. Figure 18 shows how free escape boundaries impact the resulting spectra, depending on κ analogous to Fig. 4. The higher the diffusion coefficient, the better would be the approximation of the ideal shock spectrum.The smaller the acceleration region, defined by L ± , the more particles leave through the free escape boundaries which leads to steeper spectra.

C VLUGR3 Specification
We compared the SDE approach with a finite-difference method using VLUGR3 [37,38].The computational domains were chosen to be x ∈ [−75x 0 , 75x 0 ] and s = ln p p 0 ∈ [−7.5, 15.0].The resolution was chosen to be 600 × 450 grid points.The time steps were self-determined by VLUGR3 with a given minimal time step of δ min = 10 −9 t 0 .The δ−functions for the monoenergetic injection at the shock were approximated by Gaussian functions of the form δ approx (y) = 1 √ πdy exp − y 2 dy 2 .For both δ approx (p) and δ approx (x) it was dp = dx = 0.01.At last the boundary conditions on the computational domain were given as ∂f ∂x = 0 at x ± 75x 0 and ∂f ∂s = at s = −7.5 v s = 15.0,so vanishing gradients on all four boundaries.

D Approximation of Time-dependent Solution
In Section 3.5 we compare the mean acceleration time (e.g.[25,54]) with the mean acceleration time calculated from CRPropa3.2simulations.We use the approximation given by [54] to also compare the shape of the spectrum over time.The approximation is exact in case of κ/u 2 = const., which is the solution derived by [53].Currently, this cannot be handled by CR-Propa3.2and we only compare the solution obtained with the grid-based method VLUGR3.Figure 19 shows the time evolution of spectra spatial dependence of the diffusion coefficient and spatially constant diffusion, with both energy-independent, α = 0, and energy-dependent diffusion, α = 1.The approximation and VLUGR3 solution for constant diffusion is in good agreement.With energy-dependent diffusion our results diverge from the approximation.

Figure 1 :
Figure1: Left: Spectra obtained with and without using the CandidateSplitting module at simulation time t = 400 t 0 evaluated at the shock, x = [0, 2] x 0 .The shock compression ratio is q = 4 and both simulations reproduce the expected spectral slope s = −2.The simulation with candidate splitting reaches higher energies up to 10 7 E 0 with smaller errors.Right: Weighted spectra (top) and relative error (J −2 − J)/J where J −2 is the predicted spectrum (bottom).

Figure 2 :
Figure2: One-dimensional advection field with shock at x = 0 and a shock compression ratio q = 4.The smaller the shock width compared to the integration step length, the better the ideal shock solution is approximated.

Figure 3 :
Figure 3: Time evolution of spectrum at the shock, x = [0, 2] (left) and number density, n = N dp, in the acceleration region, here defined by the free escape boundaries L ±(right).Particles with energy E 0 are continuously injected at the shock from time t = 0. Upstream speed is ũ1 = 1 with a compression of q = 4 at the shock.The diffusion coefficient is constant in space and energy, κ = 1.Free escape boundaries are at L± = ±100.Spectrum and number density resulting from simulation with CRPropa3.2(dots) are compared to the solutions obtained by integrating the transport equation with VLUGR3 (lines).

Figure 4 :
Figure4: Left: Resulting spectral slope depending on the time step ∆t in relation to the shock width L sh .The upstream speed ũ1 = 1.The compression ratio is q = 2, thus for an ideal shock we expect s = −4.With larger shock widths, the shock produces steeper spectra as calculated by[52].Slopes of the predicted stationary solutions are indicated by the colored lines for the respective shock width.Right: Resulting spectral slope depending on the ratio of upstream velocity and shock width to diffusion coefficient.Data from[44] is shown for comparison.Here ũ1 = 1, q = 4 and Lsh = 0.004 with ∆ t = 0.001 is used while κ is varied.The free escape boundaries are L± = 10.With higher diffusion coefficient the spectral slope approaches the predicted value s = −2.We already show the results for a transversal shock wave and taking perpendicular diffusion into account, which is discussed in Section 4.1.

Figure 5 :
Figure 5: Candidates are injected at t = 0 and simulated until they leave the acceleration region through free-escape boundaries L± = ±100.A finite lifetime ∆ Tscr = 20 (left) and ∆ Tscr = 80 (right) of the particle source is assumed.The spectrum at the shock x = [0, 2] is shown.

Figure 6 :
Figure 6: Left: Assuming continuous injection of candidates, the time evolution of the downstream spectrum at x = 50 is shown.Right: Candidates are injected at t = 0 and simulated until they leave the acceleration region through free-escape boundaries L± = ±100.A finite lifetime ∆ Tscr = 100 of the shock is assumed.The spectrum of particles escaping through the downstream boundary is shown.

Figure 9 :
Figure 9: Top: Mean acceleration time up to momentum p for diffusion coefficients with energy-dependence α = 0 (solid line), α = 1 (dashed line) and α = 2 (dash-dotted line).Results obtained from CRPropa simulations are plotted for comparison.Error bars show the error of the mean time.To calculate the mean acceleration time up to t = 200 t 0 , simulations were run until T = 500 t 0 .Bottom: Relative deviation of the simulated acceleration time to the prediction.With higher momentum the deviation gets larger due to finite simulation time, especially for energy-dependent diffusion.

Figure 17 :
Figure 17: Spectral slope of the stationary solution at the shock depending on the time step (indicated by color) and the ratio of upstream advective step and shock width n.The expected spectral index for finite shock widths is indicated by horizontal lines (data from[52]).

Figure 18 :
Figure18: Slope of the stationary spectra at the shock depending on the diffusion coefficient κ for different free-escape boundaries L ± .With the diffusion coefficient too large compared to the acceleration region, acceleration is less efficient and the spectral slope decreases.