11institutetext: INAF - Osservatorio Astronomico di Palermo, Piazza del
Parlamento 1, 90134 Palermo, Italy
11email<EMAIL_ADDRESS>22institutetext: Institute for Applied Problems
in Mechanics and Mathematics, National Academy of Sciences of Ukraine, Naukova
St. 3-b, 79060 Lviv, Ukraine

# Individual particle approach to the diffusive shock acceleration. Effect of
the non-uniform flow velocity downstream of the shock

O. Petruk 1122 T. Kuzyo 22

(Received …; accepted …)

###### Abstract

Aims. The momentum distribution of particles accelerated at strong non-
relativistic shocks may be influenced by the spatial distribution of the flow
speed around the shock. This phenomenon becomes evident in the cosmic-ray
modified shock, where the particle spectrum itself determines the flow
velocity profile upstream. However, what if the flow speed is not uniform
downstream as well? Hydrodynamics indicates that its spatial variation over
the length scales involved in the acceleration of particles in supernova
remnants (SNRs) could be noticeable.

Methods. In the present paper, we address this issue, initially following
Bell’s approach to particle acceleration and then by solving the kinetic
equation. We obtained an analytical solution for the momentum distribution of
particles accelerated at the cosmic-ray modified shock with spatially variable
flow speed downstream.

Results. We parameterized the downstream speed profile to illustrate its
effect on two model cases, the test particle and non-linear acceleration at
the shock.The resulting particle spectrum is generally softer in Sedov SNRs
because the flow speed distribution reduces the overall shock compression
accessible to particles with higher momenta. On the other hand, the flow
structure in young SNRs could lead to harder spectra. The diffusive properties
of particles play a crucial role as they determine the distance from the
shock, and, as a consequence, the flow speed that particles encounter
downstream. We discuss the effect of the plasma velocity gradient to be
(partially) responsible for the evolution of the radio index and for the high-
energy break visible in gamma rays from some SNRs. We expect that the effect
from the gradient of the flow velocity downstream could be prominent in
regions of SNRs with higher diffusion coefficient and lower magnetic field,
i.e. where acceleration of particles is not very efficient.

###### Key Words.:

acceleration of particles – shock waves – ISM: supernova remnants

##  1 Introduction

The theory of diffusive shock acceleration (DSA) was introduced in four
independent papers and formulated initially in two alternative approaches. The
first approach was through the kinetic equation (Krymskii, 1977; Axford et
al., 1977; Blandford & Ostriker, 1978), while the second was by using the
microscopic description (Bell, 1978a). Although Bell’s approach is quite
intuitive, most of the results related to DSA have been derived by considering
the diffusion-convection equation.

The approach of Bell (1978a) deserves more attention as it may be used to
derive the known solutions of the kinetic equation and help with new problems.
To illustrate the possibilities of the Bell approach, we obtain in Sect. 2 the
solution of Blasi (2002) for the non-linear acceleration (NLA) problem.

It has been known for a long time that the spectrum of accelerated particles
is affected by the non-uniformity of the plasma flow (e.g. Bogdan & Volk,
1983). Actually, the NLA solution reflects the non-uniform velocity
distribution upstream of the shock, where the particles of higher energies
diffuse farther upstream and are scattered back to the shock by scattering
centres whose speed differs from the speed of the centres affecting particles
of the smaller momenta which are able to penetrate to shorter distances in the
pre-shock region. This leads to the formation of a concave particle spectrum
(e.g. Malkov & Drury, 2001; Blasi, 2002) instead of a simple power-law derived
from the linear approach where the plasma speed $u$ is considered uniform both
upstream and downstream of the shock (e.g. Drury, 1983).

The common boundary condition for DSA problems is a zero gradient of the flow
velocity behind the shock. However, if $du/dx\neq 0$ in the upstream medium
affects the particle spectrum shape, what could be the effect if we also
consider a non-zero gradient downstream?

The diffusion length of particles accelerated at the shock may reach $\sim
0.1$ of the shock radius. At this distance downstream in real objects, e.g. in
supernova remnants (SNRs), the flow velocity differs from the immediate post-
shock value. Hence, we may expect a ‘downstream’ effect similar to the
‘upstream’ one from the non-linear approach. Moreover, it is known from self-
similar solutions (e.g. Chevalier, 1982; Sedov, 1959) and hydrodynamic
simulations (e.g. recent Patnaude et al. (2015); Orlando et al. (2015); Gabler
et al. (2021) for the case of the core-collapse SNe and Ferrand et al. (2019,
2021) for the type Ia SNe) that the structure of the flow undergoes essential
changes during the early evolution of SNRs. Therefore, we may expect that the
shape of the spectrum of accelerated particles should vary with time due to
the variable non-uniformity of the downstream flow.

In this paper, we address this effect in the framework of the Bell approach to
particle acceleration at shocks. Sect. 3 presents detailed motivations by
considering the spatial profiles of the plasma velocity behind the SNR shocks
from numerical simulations. Appendix A includes an animation of the flow speed
evolution. In Sect. 4, we use the generalized description of the Bell approach
from Sect. 2 to DSA with the non-uniform downstream flow. We consider two
practical cases: test-particle and non-linear acceleration (NLA) of particles.
Appendix B provides relevant details of derivations. In Sect. 5, we discuss
possible scenarios where the effect could be visible. Our conclusions are
presented in Sect. 6. Appendix C derives the same solution as in Sect. 4 by
solving the kinetic equation.

##  2 NLA in the framework of the Bell approach

Let us consider the shock reference frame with the shock at the coordinate
$x=0$. In this frame, the flow velocity of the shock $\mathbf{u}$ is in the
positive $x$ direction. Particles are accelerated by repeating the
acceleration cycle multiple times, i.e. crossing the shock from downstream
($x>0$, index ‘2’) to upstream ($x<0$, index ‘1’) and then again to
downstream. Approach of Bell (1978a) consists of finding the increase of the
particle momentum per cycle, the probability to participate in the next
acceleration cycle and then considering many repetitions. We refer the reader
to Sect. 3.3 by Jones & Ellison (1991) for the extended formulation of the
approach. In the present section, we generalize it to the modified shock.

The main idea behind the DSA is that the spatial structure of the flow
determines the change in the particle momentum $\dot{p}\equiv dp/dt$:

$\dot{p}=-\frac{p}{3}\frac{du}{dx}.$ (1)

The average momentum gain is derived by integration of this expression,

$\displaystyle\frac{\Delta
p}{p}=-\frac{1}{3}\int_{x_{a}}^{x_{b}}\displaystyle\frac{du}{dx}\frac{dx}{v_{\mathrm{x}}}$
(2)

where $v_{\mathrm{x}}=dx/dt$ is the projection of the particle velocity
$\mathbf{v}$ on the $x$-axis, and by averaging over the flux. The classic
formula

$\displaystyle\frac{\Delta
p}{p}=\displaystyle\frac{4}{v}\frac{u_{1}-u_{2}}{3}$ (3)

comes from (2) for any $x_{\mathrm{a}}<0$ and $x_{\mathrm{b}}>0$ if we assume
that the upstream $u_{1}(x)$ and the downstream $u_{2}(x)$ velocities of the
flow are uniform in their domains.

Under the NLA conditions, the pressure of accelerated particles modifies the
structure of the flow upstream of the shock, creating some profile
$u_{\mathrm{1}}(x)$. The pressure of these particles downstream is negligible
compared to the thermal pressure, so the uniformity of
$u_{\mathrm{2}}(x)=\mathrm{const}$ remains unaffected.

Let us consider the problem in a simplified and more intuitive way first. We
assume that particles with the momentum $p$ can move away from the shock to a
maximum distance $|x_{\mathrm{p}}|$ upstream where the flow speed is
$u_{\mathrm{p1}}=u_{\mathrm{1}}(x_{\mathrm{p}})$ (Blasi, 2002). Changing
$x_{a}$ to $x_{\mathrm{p}}$ in (2), we arrive at the formula (more accurate
derivation, refer to Sect. 4.1)

$\displaystyle\frac{\Delta
p}{p}=\displaystyle\frac{4}{v}\frac{u_{\mathrm{p1}}-u_{2}}{3}.$ (4)

In order to participate in the next acceleration cycle, the particle should
return to the shock from downstream against the flow direction, i.e. its
velocity component $v_{\mathrm{x}}$ should be $-v_{\mathrm{x}}>u_{2}$.
Therefore, the ratio of particle fluxes in the positive and the negative $x$
directions gives the probability of continuing acceleration for particles with
an isotropic distribution of v (Jones & Ellison, 1991)

$P=\displaystyle\left(\frac{1-u_{2}/v}{1+u_{2}/v}\right)^{2}.$ (5)

The probability to participate at least in $N$ (or larger) number of cycles is

$P=\displaystyle\prod_{i=0}^{N-1}\left(\frac{1-u_{2}/v_{i}}{1+u_{2}/v_{i}}\right)^{2}.$
(6)

Decomposition of the logarithm of $P$ to the first order on $u_{2}/v_{i}$
results in

$\ln P=\displaystyle\sum_{i=0}^{N-1}\frac{4u_{2}}{v_{i}}.$ (7)

The original approach generalised the equation for the momentum increase (4)
to $N$ cycles, and its logarithm is decomposed to the first order in
$u_{2}/v_{i}$. We, instead, express $4/v_{i}$ from the equation (4) and
substitute (7) with it. Then, considering that the ratio $\Delta
p_{i}/p_{i}\sim u_{2}/c\ll 1$ is small for relativistic particles and the
number $N$ is high, we make a transition from the summation to the
integration: $\sum\Delta p_{i}/p_{i}\rightarrow\int dp^{\prime}/p^{\prime}$.
After that, we obtain the probability for a particle to be accelerated to
momentum $p$ (or higher) from the injection momentum $p_{\mathrm{o}}$ on the
steady background of the nonuniform profile of the flow velocity upstream

$P(>\\!p)=\exp\left(-\int_{p_{\mathrm{o}}}^{p}\frac{3u_{2}}{u_{\mathrm{p^{\prime}1}}-u_{2}}\frac{dp^{\prime}}{p^{\prime}}\right).$
(8)

This probability and the definition of the isotropic distribution function
$f_{\mathrm{o}}(p)$ (Jones & Ellison, 1991)

$f_{\mathrm{o}}(p)=-\displaystyle\frac{n}{4\pi p^{2}}\frac{dP(>\\!p)}{dp}$ (9)

where index ‘o’ refers to the shock location ($x=0$), $n$ is the number
density of particles participating in the acceleration process by returning
from downstream, results in the same steady-state distribution function of
accelerated particles at the modified shock as derived by Blasi (2002) from
the alternative approach of solving the kinetic equation,

$f_{\mathrm{o}}(p)=\displaystyle\frac{\eta n_{1}}{4\pi
p_{\mathrm{o}}^{3}}\frac{3u_{1}}{u_{\mathrm{p1}}-u_{2}}\exp\left[-\int_{p_{\mathrm{0}}}^{p}\frac{3u_{\mathrm{p^{\prime}1}}}{u_{\mathrm{p^{\prime}1}}-u_{2}}\frac{dp^{\prime}}{p^{\prime}}\right]$
(10)

where $\eta$ is the injection efficiency, $n_{1}$ and $u_{1}$ should be taken
at the same point, e.g. in the immediate pre-shock location (at $x=-0$). It is
clear that the shape of the particle spectrum depends on the effective
compression factor $\sigma_{\mathrm{eff}}(p)=u_{\mathrm{p1}}/u_{2}$ which is
‘seen’ by particles with momentum $p$. The function (10) transforms to the
power-law test-particle version by assuming that the flow velocity upstream is
uniform, meaning that the same convergence of the flow for particles of any
momentum: $\sigma_{\mathrm{eff}}(p)=u_{1}/u_{2}$. The nonlinearity of this
solution lies in the fact that $\sigma_{\mathrm{eff}}(p)$ itself depends on
$f(p)$, which requires another equation for these two functions. Blasi (2002)
developed a way to relate $\sigma_{\mathrm{eff}}(p)$ to $f(p)$ by considering
the hydrodynamics of the flow. In their notation,
$\sigma_{\mathrm{eff}}(p)=R_{\mathrm{tot}}U(p)$ where
$R_{\mathrm{tot}}=u_{0}/u_{2}$ represents the total compression factor of the
shock and $u_{0}$ is the flow speed far upstream.

Figure 1: Profiles of the flow velocity $v(r)$ in the observer reference frame
(left) and $u(x)$ in the shock reference frame (right) for the three models of
SNRs. Upper row: three SNe models at 500 yrs. Lower row: SNIa model at the
ages 20, 200, and 1000 yrs. On the plots to the left, the explosion centre is
at $r=0$ with the zero velocity $v$ there. On the plots to the right, the
horizontal axes display $x=1-r/R$. Therefore, negative values correspond to
the upstream region, whereas positive values correspond to the downstream
region. The centre of the explosion on the right-hand side plots is at $x=1$,
and the velocity $u$ at this point equals the upstream velocity $u_{1}$. We
conducted the numerical simulations for the unmodified shocks, therefore $u$
is spatially constant upstream. The dashed line represents the velocity
profile from the Sedov (1959) solution.

##  3 Profiles of the flow velocity

Studies of the cosmic-ray modified shocks demonstrate that the spatial non-
uniformity of the plasma velocity upstream of the shock affects the shape of
the momentum distribution of accelerated particles. It is natural to ask how
the downstream flow structure modifies the particle spectrum. Numerical
simulations of CR acceleration which accounts for the complex hydrodynamic
structure of the flow inside SNRs demonstrate a variety of shapes of the
volume-averaged CR spectra (e.g. Telezhinsky et al., 2012, 2013; Brose et al.,
2020). In the present paper, we would like to consider the effect of the non-
uniform flow in the vicinity of a single shock where the CR acceleration
happens.

The length scales involved in the particle acceleration at the forward shock
are of the order of the diffusion length which may eventually be up to a few
times $0.1R$ where $R$ is the SNR radius. Most of the mass swept up by the
forward shock is located right behind the adiabatic shock within about $10\%$
of the radius. The ratio of densities is $\rho(0.9R)/\rho(R)=0.49$ for the
Sedov shock in the uniform medium. The variation of the flow velocity in the
laboratory reference frame is not so strong, $v(0.9R)/v(R)=0.82$, but it is
not negligible.

The leading motivation for the present study is the fact that the velocity of
the flow $u(x)$ is not uniform downstream of the shock on the length scale
involved in the particle acceleration. It varies spatially from the value
$u_{2}$ at $x=+0$ in accordance with the profile determined by the
hydrodynamics. However, the length scale is small enough to avoid considering
the effect of the shock curvature on particle acceleration. In fact, the
distance $0.1R$ on the SNR surface corresponds to the sector with a central
angle in $\arctan(0.1)=5.7$ degrees. Therefore, we can consider the one-
dimensional acceleration problem within such a small sector with a plane
shock.

To understand the profiles $u(x)$ we can expect downstream, we suggest
referring to our previous paper (Petruk et al., 2021). There, we have
investigated the MHD evolution of the remnants of three supernova models,
namely, thermonuclear SNIa and two core-collapse SNIc, SNIIP. We have adopted
a uniform ISM to model the type SNIa explosion, whereas the core-collapse SNe
evolved in a CSM bubble with a density shaped by $\rho_{\mathrm{CS}}=Dr^{-2}$
where the value of $D$ depends on the explosion model. Core-collapse SNe
produces a range of explosion conditions, unlike the thermonuclear explosion
of a WD star (SNIa event). To account for this, we have selected the two most
representative cases of BSG star (SNIc) and RSG star (SNIIP) undergoing the
core collapse. The SNIc and SNIIP models differ in the effective explosion
mass and CSM number density of a wind-blown bubble. We have performed all
numerical simulations with PLUTO MHD code (Mignone et al., 2007, 2012). The
spherically symmetric uniform grid spanned 20 000 zones, providing a very
detailed spatial resolution during the entire simulation period. Petruk et al.
(2021) provides more details on the technical model setup as well as initial
conditions for the explosion models.

Fig. 1 shows the spatial distributions of the flow speed for the three models
at the age of $500{\,\rm yrs}$ and for three different SNR ages in the case of
the SNIa model. The plots in the left column correspond to the laboratory
frame, and those in the right column correspond to the shock frame. The
velocity profile $u_{2}(x)$ may rise or decline with the distance from the
shock. The distribution evolves with time, e.g. changing from a decreasing to
an increasing function of $x$ in the SNIa model (Fig. 1 bottom right plot).
Animation in Appendix A shows the detailed evolution of the velocity
distributions up to $10\,000{\,\rm yrs}$ in these three models.

##  4 Particle acceleration on the background of the non-uniform flow
velocities upstream and downstream

###  4.1 Distribution function

Let us consider the particle acceleration on the flow where the velocity
profiles are non-uniform both upstream and downstream of the shock. We start
again from the equation (1). To calculate the increase of momentum per
acceleration cycle, we need to spatially separate the integration

$\Delta p=\int\limits_{-\infty}^{-0}\dot{p}(x){\cal
X}_{1}(x)\frac{dx}{v_{\mathrm{x}}}+\int\limits_{-0}^{+0}\dot{p}(x){\cal
X}_{0}(x)\frac{dx}{v_{\mathrm{x}}}+\int\limits_{+0}^{+\infty}\dot{p}(x){\cal
X}_{2}(x)\frac{dx}{v_{\mathrm{x}}}$ (11)

where the functions ${\cal X}_{i}(x)$ mean of the probability of finding a
cosmic ray with a given momentum at a given coordinate $x$. The momentum of a
particle is modified not only at the shock transition (where the flow speed
changes from $u_{1}(0)$ to $u_{2}(0)$) but also in the convergent (or
divergent) flows upstream and downstream.

The formula (11) considers one-dimensional flow. We will use it along with the
velocity profiles calculated in Sect. 3 using the spherical geometry. There
are two important factors we have to keep in mind. First, the centre of SNR is
at the point $x=R$ in the shock reference frame. Therefore, the upper limit in
the last integral should be $R$. However, for the sake of generalisation, we
use the notation $+\infty$ since the acceleration region is much smaller than
the radius $R$ of SNR. Second, there is a geometrical effect. In the
laboratory frame, the integration over the SNR volume is given by $r^{N}dr$
where $N=0,1,2$ for the plane, cylindrical, and spherical shocks,
respectively. By integrating over $x$ in the shock reference frame, we assume
that this effect is accounted for by the probabilities ${\cal X}_{i}(x)$.

The spatial distribution of particles is given by the functions $f(x,p)$.
Therefore, the probability for the accelerated particle to be at a coordinate
$x$ is

${\cal X}_{i}(x)=\frac{f_{i}(x,p)}{f_{\mathrm{o}}(p)}$ (12)

where $f_{\mathrm{o}}(p)\equiv f(0,p)$ is the distribution function at the
shock. Obviously, ${\cal X}_{0}(x)=1$.

Equation (1) reads $\dot{p}={p}\left(u_{1}(0)-u_{2}(0)\right)\delta(x)/3$ for
the second integral in (11). Then, the expression (11) may be written as

$\frac{\Delta
p}{p}=\frac{1}{3}\frac{u_{\mathrm{p1}}-u_{\mathrm{p2}}}{v_{\mathrm{x}}}$ (13)

with the notations

$u_{\mathrm{p1}}=u_{1}(0)-\int_{-\infty}^{-0}{\cal X}_{1}(x)\frac{du}{dx}dx,$
(14) $u_{\mathrm{p2}}=u_{2}(0)-\int_{+\infty}^{+0}{\cal
X}_{2}(x)\frac{du}{dx}dx,$ (15)

The transition to the classical test-particle case is provided by $du/dx=0$.

The probabilities ${\cal X}_{i}(x)$ are related to particles that participate
in the acceleration around the shock. They do not account for particles
escaping upstream or advected downstream. In particular, by equating the
diffusive and convective flows in the problem of the test-particle
acceleration, we have (Drury, 1983)

${\cal
X}_{1}(x)=\exp\left[\int_{0}^{x}\frac{u_{1}(x^{\prime})dx^{\prime}}{D_{1}(x^{\prime},p)}\right]$
(16)

where $D_{1}$ is the diffusion coefficient. For the NLA problem (Amato &
Blasi, 2005; Caprioli et al., 2010), the expression for this probability is
more complex. However, the same references suggest its accurate approximation:

${\cal
X}_{1}(x)=\exp\left[\frac{s(p)}{3}\frac{\sigma_{\mathrm{s}}-1}{\sigma_{\mathrm{s}}}\int_{0}^{x}\frac{u_{1}(x^{\prime})dx^{\prime}}{D_{1}(x^{\prime},p)}\right]$
(17)

where $\sigma_{\mathrm{s}}=u_{1}(0)/u_{2}(0)$ is the sub-shock compression
ratio. The probability ${\cal
X}_{1}(x)\simeq\exp(u_{1}x/D_{1})\simeq\exp(-x/x_{\mathrm{p}})$ where
$x_{\mathrm{p}}\simeq-D_{1}/u_{1}$ is the coordinate that corresponds to the
diffusion length for particles with momentum $p$. This may roughly be
approximated by the Heaviside step function ${\cal X}_{1}(x)\approx{\cal
H}(x-x_{\mathrm{p}})$. This is why we used $u_{\mathrm{p}}\approx
u_{1}(x_{\mathrm{p}})$ in Sect. 2.

In the test particle problem, a formula similar to equation (16) may be
derived in the same way for the region behind the shock. In the future, a
derivation of $f_{2}(x,p)$ and ${\cal X}_{2}(x)$ in a more general case,
considering NLA effects and $du_{2}/dx\neq 0$ should be taken into account. In
the present paper, we assume that the probability downstream is similar to
that upstream, Which means that ${\cal
X}_{2}(x)\simeq\exp(-x/x_{\mathrm{p}})\approx{\cal H}(x_{\mathrm{p}}-x)$ and
thus $u_{\mathrm{p2}}\approx u_{2}(x_{\mathrm{p}})$. In fact, at a given point
$x$ downstream, only a fraction $\xi$ of particles with momentum $p$ can
return to the shock and continue acceleration. A fraction $1-\xi$ is advected
downstream from this point and never returns to the shock. The fraction $\xi$
decreases with distance from the shock because the particles can return from
about one diffusion length. This distance is $x_{\mathrm{p}}\simeq
D_{2}(p)/u_{2}\propto p^{\alpha}/u_{2}$ with positive $\alpha$, i.e. there are
progressively fewer particles with momentum $p$ able to return being farther
from the shock.

The equation (13) gives the average $\Delta p/p$ for the first half-cycle
(upstream-downstream) of acceleration. When we also consider the second half-
cycle (downstream-upstream) and take the average over the flux, we obtain the
expression for the average momentum increase

$\displaystyle\frac{\Delta
p}{p}=\displaystyle\frac{4}{v}\frac{u_{\mathrm{p1}}-u_{\mathrm{p2}}}{3}.$ (18)

The presented equation accounts for the fact that the plasma velocity profiles
are not constant upstream and downstream, contrary to what is assumed in the
classical picture of the test-particle acceleration at the shock. This
equation shows that the momentum change depends on the difference between the
two utmost values of the flow speed that are ‘reachable’ for accelerating
particles with a given momentum.

In the Sect. 2, the probability for a particle to cross the shock from
downstream to upstream and to participate consequently in the next
acceleration cycle is given by the ratio of the particle fluxes to the
positive and to the negative $x$ directions for the isotropic distribution of
v and uniform profile of the plasma velocity downstream
$u_{2}(x)=\mathrm{const}$ (Jones & Ellison, 1991). In our case, $u_{2}$ varies
with the distance from the shock.

Let us consider a simple and intuitive assumption that particles with a given
momentum $p$ can reach the maximum coordinate $x_{\mathrm{p2}}$ downstream of
the shock and this distance is bigger for particles with larger momenta.
Particles must return to the shock from this point to continue acceleration.
Therefore, a component of a particle velocity $v_{\mathrm{x}}$ must be larger
than the flow speed $u_{\mathrm{p2}}\approx u_{2}(x_{\mathrm{p}})$ there. With
these assumptions, the we can write the expression for the probability from
(5) as

$P=\displaystyle\left(\frac{1-u_{\mathrm{p2}}/v}{1+u_{\mathrm{p2}}/v}\right)^{2}.$
(19)

The probability for $n$ cycles is

$P=\displaystyle\prod_{i=0}^{n-1}\left(\frac{1-u_{\mathrm{p2,\mathit{i}}}/v_{i}}{1+u_{\mathrm{p2,\mathit{i}}}/v_{i}}\right)^{2}.$
(20)

Decomposition of the logarithm to the first order on
$u_{\mathrm{p2,\mathit{i}}}/v_{i}$ gives

$\ln
P=-\displaystyle\sum_{i=0}^{n-1}\frac{4u_{\mathrm{p2,\mathit{i}}}}{v_{i}}$
(21)

Appendix B presents a more general derivation of this expression.

Let us now substitute this equation with $4/v_{i}$ from (18) and convert the
sum $\sum\Delta p_{i}/p_{i}$ to the integral $\int dp^{\prime}/p^{\prime}$. As
a result, we derive the probability of a particle being accelerated to
momentum $p$ (or higher) on a steady background of the nonuniform profiles of
the flow velocity upstream and downstream

$P=\exp\left(-\int_{p_{\mathrm{o}}}^{p}\frac{3u_{\mathrm{p^{\prime}2}}}{u_{\mathrm{p^{\prime}1}}-u_{\mathrm{p^{\prime}2}}}\frac{dp^{\prime}}{p^{\prime}}\right).$
(22)

The isotropic distribution function $f_{\mathrm{o}}(p)$ is defined by equation
(9) with $n$ the number density of all accelerated particles, making it a
fraction $n=\eta n_{0}$ of all the incoming particles. The flux is conserved
$n(x)u(x)=\mathrm{const}$. Therefore, we can take $n$ at any point, e.g. we
use $n=\eta n_{1}u_{1}/u_{\mathrm{p2}}$ to resemble the formula (10). Thus,
the definition (9) finally yields the distribution function

$f_{\mathrm{o}}(p)=\frac{\eta n_{1}}{4\pi
p_{\mathrm{o}}^{3}}\frac{3u_{1}}{u_{\mathrm{p1}}-u_{\mathrm{p2}}}\exp\left[-\int_{p_{\mathrm{o}}}^{p}\frac{3u_{\mathrm{p^{\prime}1}}}{u_{\mathrm{p^{\prime}1}}-u_{\mathrm{p^{\prime}2}}}\frac{dp^{\prime}}{p^{\prime}}\right].$
(23)

Known solutions may be easily derived from the equation (23). Namely,
$u_{2}(x)=\mathrm{const}$ results in the non-linear solution of Blasi (2002).
If $u_{1}(x)=\mathrm{const}$ as well, then we have the common test-particle
power-law spectrum for the momentum distribution of accelerated particles.

For the cosmic-ray modified shock, the way to calculate $u_{\mathrm{p1}}$ with
an integro-differential equation which accounts for the regulation of the flow
hydrodynamics by $f_{\mathrm{o}}(p)$ is developed by Blasi (2002); Blasi et
al. (2005). In these references, the flow is assumed to follow the adiabatic
law in the pre-shock region. Amato & Blasi (2006); Caprioli et al. (2009)
examine the impact of other effects in the precursor, such as turbulent
heating and magnetic field amplification, on the plasma velocity profile
upstream.

The equation (23) describes the momentum distribution of accelerated particles
at the shock and is derived in the frame of the approach of Bell (1978a) to
the particle acceleration. The solution of the same problem derived in the
frame of an alternative approach, namely by solving the kinetic equation for
DSA, is presented in Appendix C.

###  4.2 Spectral index

The ‘local’ spectral index $s(p)$, i.e. the local slope of the particle
distribution $f_{\mathrm{o}}(p)$ at a given momentum $p$, is defined as

$s(p)=-\frac{d\ln f_{\mathrm{o}}}{d\ln p}.$ (24)

Applying it to the function (23), we have

$s(p)=\frac{3u_{\mathrm{p1}}}{u_{\mathrm{p1}}-u_{\mathrm{p2}}}+\frac{d\ln(u_{\mathrm{p1}}-u_{\mathrm{p2}})}{d\ln
p}.$ (25)

This expression represents the spectral index for particles accelerated around
the shock with spatially non-uniform distribution of the flow speed $u(x)$
upstream and downstream.

Figure 2: Schematic representation of the structure of the flow velocity
$u(x)$ (solid line). The shock is located at the coordinate $x=0$. The
upstream profile of $u(x)$ is shown for the case of the cosmic-ray modified
shock.

Let us introduce the effective compression factor that is ‘seen’ by particles
with the momentum $p$ (as shown in Fig. 2):

$\sigma_{\mathrm{eff}}(p)=\frac{u_{\mathrm{p1}}}{u_{\mathrm{p2}}}\approx\frac{u_{1}(x_{\mathrm{p1}})}{u_{2}(x_{\mathrm{p2}})}.$
(26)

The spectral index for the cosmic-ray modified shock (it is given by equation
(8) in Blasi, 2002) comes from (25) by assuming $u_{2}=\mathrm{const}$ and
therefore $\sigma_{\mathrm{eff}}=u_{\mathrm{p1}}/u_{2}$:

$s(p)=\frac{3\sigma_{\mathrm{eff}}}{\sigma_{\mathrm{eff}}-1}\left(1+\frac{1}{3}\frac{p}{\sigma_{\mathrm{eff}}}\frac{d\sigma_{\mathrm{eff}}}{dp}\right).$
(27)

Equation (25) gives the classic test-particle index

$s=\frac{3\sigma_{\mathrm{s}}}{\sigma_{\mathrm{s}}-1}$ (28)

for $u_{1}=\mathrm{const}$ and $u_{2}=\mathrm{const}$, namely, by substitution
$\sigma_{\mathrm{eff}}=\sigma_{\mathrm{s}}=u_{1}/u_{2}$. In the general case
when $u_{1}(x)\neq\mathrm{const}$ and $u_{2}(x)\neq\mathrm{const}$, the
contribution to $s(p)$ from the first term in (25) is the same as for the
unmodified shock (28) but with the effective compression factor
$\sigma_{\mathrm{eff}}$ instead of $\sigma_{\mathrm{s}}$. This contribution is
modified somehow by the second term in (25), which represents the ‘local’ (in
the momentum space) slope of the difference $u_{\mathrm{p1}}-u_{\mathrm{p2}}$.

Let us consider the effects of the spatially nonuniform $u_{2}(x)$ on a few
practical cases.

###  4.3 Benchmark problems

####  4.3.1 Test particles at Sedov shock

We start from the acceleration of the test particles, i.e.
$u_{1}=\textrm{const}$, at adiabatic shock with the Sedov profile for
$u_{2}(x)$.

For simplicity and clarity, we use the approximation by Taylor (1950) for the
distribution of the flow velocity downstream of the Sedov (1959) shock. It is

$v(r)=v_{2}(R)r/R$ (29)

in the laboratory frame, the immediate post-shock value is $v_{2}(R)=3V/4$,
$V$ is the shock speed, $r$ the distance from the explosion center, $R$ the
shock radius. The accuracy of the approximation is better than $4\%$ for
$>0.7R$ (figure 1 in Petruk, 2000). In the shock reference frame, the profile
is

$u_{2}(x)=u_{2}(0)\cdot(1+3x/R)$ (30)

where $u_{2}(0)=u_{1}/\sigma_{\mathrm{s}}$. We calculate $u_{\mathrm{p2}}$
with expression $u_{\mathrm{p2}}\approx u_{2}(x_{\mathrm{p}})$ which follows
from the definition (15) and ${\cal X}_{2}(x)\approx{\cal
H}(x_{\mathrm{p}}-x)$. In order to parametrise the dependence of the diffusion
length $x_{\mathrm{p}}$ on $p$, we take

$x_{\mathrm{p}}(p)\propto\frac{D_{2}(p)}{u_{2}(0)}$ (31)

where $D$ is the spatially constant diffusion coefficient of the form
$D_{2}\propto p^{\alpha}$ with $\alpha\geq 0$. Therefore,

$u_{\mathrm{p2}}=u_{2}(0)\cdot(1+3\bar{x}_{\mathrm{p}}),\qquad\bar{x}_{\mathrm{p}}(p)=\bar{x}_{\mathrm{pm}}\cdot\left({p}/{p_{\mathrm{max}}}\right)^{\alpha}$
(32)

where $\bar{x}_{\mathrm{pm}}\equiv x_{\mathrm{p}}(p_{\mathrm{max}})/R$ is the
diffusion length for particles with the maximum momentum $p_{\mathrm{max}}$
normalized to $R$. In such a model, we can evaluate the spectral index
analytically:

$s(p)=\frac{3\sigma_{\mathrm{s}}-3\alpha\bar{x}_{\mathrm{p}}(p)}{\sigma_{\mathrm{s}}-1-3\bar{x}_{\mathrm{p}}(p)}.$
(33)

Figure 3: The local slopes $s(p)$ and the distribution functions
$f_{\mathrm{o}}(p)/n_{\mathrm{o}}$ for the unmodified Sedov shock. The
parameters are $\bar{x}_{\mathrm{pm}}=0.1$ (blue lines),
$\bar{x}_{\mathrm{pm}}=0.05$ (red lines), $\alpha=1$ (solid lines),
$\alpha=1/3$ (dashed lines), $p_{\mathrm{max}}=10^{4}m_{\mathrm{p}}c$. The
classic value (for $u_{2}=\mathrm{const}$) is $s=4$ and is independent of $p$,
$\alpha$ and $\bar{x}_{\mathrm{pm}}$. Vertical lines: left – the energy of
electrons for which the radio index is calculated in Fig. 6; right – protons
with energy $240{\,\rm GeV}$ (see Sect. 5.2).

To see the sensitivity of the problem to $\bar{x}_{\mathrm{pm}}$, we consider
a few values of this parameter. We take these values for granted in the
present section and explore their role further in Sect. 5.3.

In Fig. 3, we can see the local spectral index and the distribution function
for the model we are considering. It is evident that the effect of the
function $u_{2}(x)$, which increases with $x$, leads to a higher $s$ for the
higher $p$. (Note that in the observer frame, the flow speed $v$ decreases
downstream with distance from the shock.) It is becase the effective
compression factor $\sigma_{\mathrm{eff}}$ decreases with growing $p$ (Fig.
2). Therefore, the momentum distribution function $f_{\mathrm{o}}(p)$ descends
more rapidly than a common power-law (Fig. 3).

It’s worth noting that the shapes of the particle distributions shown in Fig.
3 are only due to the plasma velocity gradient downstream of the shock. We
applied only the sharp high energy cutoff to the spectrum of particles ${\cal
H}(p_{\mathrm{max}}-p)$, where ${\cal H}(x)$ is the Heaviside step function.
For comparison, the spectrum is just the power law up to the maximum momentum
$N(p)=Kp^{-s}{\cal H}(p_{\mathrm{max}}-p)$ in case of the constant flow speed
$u_{2}(x)=\mathrm{const}$.

The classical test-particle index for the case of $u_{2}=\mathrm{const}$ is
$s=4$. From Fig. 3, we observe that the local index $s(p)$ for the Sedov
velocity profile is close to this value for particles with low momenta $p$,
i.e. for the particles which are accelerated from the small distances
$\bar{x}_{\mathrm{p}}\approx 0$, where deviations of the flow velocity from
$u_{2}(0)$ are minor. The index $s(p)$ increases with $p$ for Sedov shock, and
its maximum value happens for particles around $p_{\mathrm{max}}$ which are
able to return to the shock from deeper downstream locations, i.e. from around
$\bar{x}_{\mathrm{pm}}$.

The effect described in the present paper is sensitive to $\alpha$. For a
fixed $\bar{x}_{\mathrm{pm}}$, the effect is more prominent for smaller
$\alpha$ (Fig. 3) because particles with low momenta may return to the shock
from deeper $\bar{x}_{\mathrm{p}}$ downstream regions.

Electrons in the magnetic field $B=10{\,\rm\mu G}$ emitting most of their
energy through the synchrotron radiation at the radio frequency $5{\,\rm GHz}$
have energy $10{\,\rm GeV}$. It is apparent from Fig. 3 that, at this energy,
the spectral index differs by a few percent from the usual value. For example,
the radio index is $0.52$ for $\alpha=1/3$ and $\bar{x}_{\mathrm{pm}}=0.1$
instead of classic $0.5$. The differences are larger for smaller $\alpha$. In
the limit of a momentum-independent $D$ (i.e. $\alpha=0$), the spectral index
is also momentum independent. In such a toy model, its value depends on the
distance $x_{\mathrm{p}}$ which is the same for any $p$, as it follows from
the equation (32). Numerically, for the Sedov profile of $u_{2}(x)$, the
values are $s=4.4$ for $\bar{x}_{\mathrm{pm}}=0.1$ and $s=4.2$ for
$\bar{x}_{\mathrm{pm}}=0.05$. The radio index at $5{\,\rm GHz}$ is $0.7$ and
$0.6$ respectively.

####  4.3.2 Cosmic-ray modified Sedov shock

Let us consider the acceleration of particles on the cosmic-ray modified shock
with the upstream profile $u_{1}(x)\neq\mathrm{const}$ and the downstream
distribution for $u_{2}(x)$ given by equation (30), i.e. like in a Sedov SNR.
The goal is to compare the spectra of particles with the case
$u_{2}(x)=\mathrm{const}$.

The essence of the non-linearity of the cosmic-ray modification of the shock
is in the dependence of the plasma speed profile upstream of the shock
$u_{1}(x)$ on the particle distribution function $f_{\mathrm{o}}(p)$. If we
have $u_{2}(x)=\mathrm{const}$ downstream of the shock, then the NLA solution
of Blasi (2002) holds. In our problem, we should account for the spatial
variation $u_{2}(x)$ as well. We show above that the non-uniform profile
$u_{2}(x)$ affects the shape of the distribution function $f_{\mathrm{o}}(p)$.
Therefore, the profile $u_{1}(x)$ upstream will be modified by cosmic rays
differently compared to the uniform $u_{2}(x)$ distribution. In turn, this
modified upstream profile affects $f_{\mathrm{o}}(p)$, too.

Figure 4: Plots of $s(p)$ and $f_{\mathrm{o}}(p)/n_{\mathrm{o}}$ for the
cosmic-ray modified shock with $u_{2}(x)=\mathrm{const}$ (red line) and with
$u_{\mathrm{p2}}(p)$ given by the equation (32) with $\alpha=1$ (solid blue
line), $\alpha=1/3$ (dashed blue line), $\bar{x}_{\mathrm{pm}}=0.1$,
$p_{\mathrm{max}}=10^{4}m_{\mathrm{p}}c$.

We calculated the NLA solution for such a problem numerically in a self-
consistent way by following the approach described by Blasi (2002); Blasi et
al. (2005). Namely, we found the function $u_{\mathrm{p1}}(p)$ by solving
iteratively the integro-differential equation which depends on
$f_{\mathrm{o}}(p)$, accounts for the adiabatic evolution of the flow in the
shock precursor, for the relation between the injection efficiency and the
compression factor at the sub-shock and – our addition – for the flow velocity
structure downstream. Then the distribution function $f_{\mathrm{o}}(p)$ may
be found from the equation (23). The iterations are repeated until the
boundary condition $u_{\mathrm{p1}}(p_{\mathrm{max}})=u_{0}$ is satisfied. We
considered the following set of parameters: the Mach number ${\cal
M}_{\mathrm{o}}=100$, the flow speed far upstream $u_{0}=6.5\times
10^{8}{\,\rm cm/s}$, the maximum momentum
$p_{\mathrm{max}}=10^{4}m_{\mathrm{p}}c$ and the injection efficiency
$\eta=1.30\times 10^{-4}$. The red line on Fig. 4 shows the solution for
$u_{2}(x)=\mathrm{const}$ while the blue lines correspond to models with the
velocity gradient downstream of the shock.

Qualitatively, we can understand the behavior of the momentum distribution of
accelerated particles from Fig. 2. In this scheme, we show a model of the
cosmic-ray modified shock (it is represented by the velocity profile in the
$x<0$ domain) with increasing function $u_{2}(x)$ downstream (in the $x>0$
domain). If we first consider the uniform $u_{2}(x)$ (dashed grey horizontal
line on Fig. 2), we can observe that particles with the lowest momenta ‘see’
the compression $\sigma_{\mathrm{s}}$. Particles with larger momenta are
affected by continuously higher compression $\sigma_{\mathrm{eff}}(p)$ which
raises due to the velocity profile upstream. It is the main characteristic of
the NLA solution shown by the red line in Fig. 4: the higher the effective
compression, the lower the spectral index $s$ and the harder the distribution
$f_{\mathrm{o}}(p)$. However, the increasing function $u_{2}(x)$ (solid line
downstream in Fig. 2) progressively weakens this trend for higher momenta. As
a result, the spectrum’s shape may be less concave compared to the pure NLA
case or have a complex shape with an inflection point (blue lines in Fig. 4).
In contrast, the decreasing $u_{2}(x)$ causes stronger effective compression
compared to the uniform $u_{2}$ downstream and strengthens the NLA concavity.

The effective compression factors of are typically high in the NLA problem.
They often lead to the steep shapes of the emission spectra which are not
observed in SNRs. To make $\sigma_{\mathrm{eff}}$ lower compared to the pure
NLA case, Caprioli et al. (2009) consider the relative Alfv́enic speed
$v_{\mathrm{A}}$ of the scattering centres (i.e. $u+v_{\mathrm{A}}$ instead of
$u$ in the diffusion-convection equation). Our model with a rising flow speed
in the shock reference frame could be a natural way to lower
$\sigma_{\mathrm{eff}}$ and make the spectrum flatter in circumstances where
the cosmic-ray shock modification is an essential component.

It’s worth noting that the injection efficiency we consider for Fig. 4 is not
high, and the shock modification is not very strong. Indeed, the compression
factor $\sigma_{\mathrm{tot}}\equiv u_{0}/u_{2}(0)$ is $5.434$ for the red
line, $5.359$ for the blue solid line and $4.897$ for the blue dashed line. If
the back-reaction of accelerated particles is higher, e.g.
$\sigma_{\mathrm{tot}}$ is about $10$ or more, then the role of the upstream
profile $u_{1}(x)$ in the shaping $f_{\mathrm{o}}(p)$ dominates the effect
from the downstream distribution $u_{2}(x)$.

##  5 Discussion

We saw in the previous section that the spatial variation of the flow speed
downstream of the shock results in deviation of the momentum distribution of
accelerated particles from the classic power law. The spectrum of emission
produced by these particles should also be affected. Could we see such effect
in some observables? How strong could it be, if any? Let us look for an answer
in the frame of the test-particle acceleration, i.e. to assume in the present
section that $u_{1}=\textrm{const}$.

###  5.1 Evolution of the radio index of young SNRs

Previously, we explored magneto-hydrodynamic properties of young SNRs on
examples of remnants from the three types of supernovae (Petruk et al., 2021).
Sect. 3 shows (see also an animation in Appendix A) that the temporal changes
in the plasma velocity structure downstream of the shock $u_{2}(x,t)$ are
essential during the development of SNR from the early stage (Chevalier, 1982;
Dwarkadas & Chevalier, 1998) to the Sedov (1959) epoch. These changes may
cause time variations in the momentum distribution of particles and,
consequently, in the evolution of the radio index of SNRs.

Fig. 5 shows the evolution of the effective compression factor for particles
with the maximum momentum
$\sigma_{\mathrm{eff}}(p_{\mathrm{max}})=u_{1}/u_{2}(\bar{x}_{\mathrm{pm}})$
with $\bar{x}_{\mathrm{pm}}=0.1$ for three models of SNRs. The flow dynamics
are different between the models. Specifically, the flow speed within
$\bar{x}_{\mathrm{pm}}$ from the shock is almost the same in the SNIIP model,
which is one of the representative cases for the DSA theory that assumes
$u_{2}(x)=\mathrm{const}$. On the other hand, the SNIa model demonstrates
prominent changes in conditions for the shock particle acceleration during its
evolution.

These differences may be present in the evolution of the radio index. To
expose the effect, we approximate the profile of $u_{2}(x,t)$ for the
downstream region close to the shock by the linear function with a variable
coefficient $k$:

$u_{2}(x,t)=u_{2}(0)\cdot\left(1+k(t)x/R\right).$ (34)

This approximation is justifiable by the shape of the $u_{2}(x)$ profiles
within the interval of $x$ from zero to $\sim 0.1R$ (Fig. 1 and animation in
Appendix A). Then the function $u_{\mathrm{p2}}(p)$ is

$u_{\mathrm{p2}}(p)=u_{2}(0)\cdot\left(1+k(t)\bar{x}_{\mathrm{p}}(p)\right)$
(35)

instead of (32) and the local slope of the particle distribution is, instead
of (33),

$s(p,t)=\frac{3\sigma_{\mathrm{s}}-k(t)\alpha\bar{x}_{\mathrm{p}}(p)}{\sigma_{\mathrm{s}}-1-k(t)\bar{x}_{\mathrm{p}}(p)}.$
(36)

Figure 5: The effective compression factor
$u_{1}/u_{2}(\bar{x}_{\mathrm{pm}})$ with $\bar{x}_{\mathrm{pm}}=0.1$ for
different ages. Models of SNRs are the same as in Sect. 3.  Figure 6: The
evolution of the radio index at $5{\,\rm GHz}$ for the three models of SNRs
and electrons with energy $10{\,\rm GeV}$ (so, the magnetic field strength is
$10{\,\rm\mu G}$), $p_{\mathrm{max}}=10^{4}m_{\mathrm{p}}c$,
$\bar{x}_{\mathrm{pm}}=0.1$, $\alpha=1/3$ (solid lines). The dashed line
represents the same model as the solid line but with artificially taken
$\alpha=1/5$.

The evolution of the radio index $(s-3)/2$ calculated with this expression for
$s(p,t)$ is shown in Fig. 6. In the model SNIa displaying the most prominent
changes in the plasma speed’s structure (Fig. 5), the radio index varies from
$0.46$ to $0.54$ during about $1000{\,\rm yrs}$ (red solid line, $s$ changes
from $1.9$ to $2.1$ respectively). The range for the radio index variation
could be even wider if particles probe deeper interiors of SNR during
acceleration (i.e. when ($x_{\mathrm{pm}}>0.1R$ or $\alpha$ smaller in the
assumed dependence $x_{\mathrm{p}}(p)$, Eq. 32), as shown by the dashed line
in Fig. 6. It is also consistent with Fig. 3.

We can conclude that the spectral index in the SNIa model is smaller than the
value $0.5$ up to SNR age $\approx 100{\,\rm yrs}$ and higher at later times,
as shown in Fig. 6. The reason is given by Equation (36), which shows that the
classic value $s=4$ (and the radio index $0.5$) corresponds to $k=0$.
Furthermore, this formula reveals that $s$ could be below $4$ when $k<0$
(which happens when the flow velocity $u_{2}$ decreases with $x$, and thus the
effective compression factor $u_{1}/u_{2}(\bar{x}_{\mathrm{pm}})>4$, as seen
in Fig. 5) or above $4$ when $k>0$.111This is valid for
$\alpha<3\sigma_{\mathrm{s}}/(\sigma_{\mathrm{s}}-1)$ that is always the case.

It’s worth noting that the left vertical line in Fig. 3 corresponds to
electrons that mostly determine the radio index shown in Fig. 6. We see that
the deviation of their spectrum from the power law is minor. However, the
deviation is more significant for particles with higher momenta that emit
X-rays and gamma-rays.

###  5.2 High-energy spectral break

A spectral feature is observed in the gamma-ray spectra of several SNRs. It is
attributed to a change in the slope of the proton spectrum which is typically
referred to as the broken power law. The change in the slope $\Delta s\approx
0.4-1$ and should be added to proton spectrum models around the break momentum
$p_{\mathrm{br}}$, which has values from $p_{\mathrm{br}}=10{\,\rm GeV}/c$ to
$p_{\mathrm{br}}=240{\,\rm GeV}/c$ to fit the observed spectra of some SNRs
(e.g. Ackermann et al., 2013; H. E. S. S. Collaboration et al., 2015, 2018).
Such a feature is also visible in radio (Loru et al., 2019) and microwave
(Planck Collaboration et al., 2016) spectra of SNRs, although it has roots in
the momentum distribution of electrons.

Our results suggest that the particles accelerated at the shock with spatially
variable flow velocity downstream have a momentum distribution with a smooth
change of slope over a wide range of momenta. The radio-to-X ray synchrotron
spectrum is convex for evolved SNRs (when $u_{2}(x)$ is increasing function,
Sect. 3) that could mimic the ‘break’ if it is not very sharp. Our
calculations in Sect. 4.2 show that the slope changes between the particle
energies $1{\,\rm Gev}$ and $300{\,\rm Gev}$ as $\Delta s=0.11$ (from
$s=4.02$) or $\Delta s=0.30$ (from $s=4.3$) for Sedov shock and parameters
$(\bar{x}_{\mathrm{pm}},\alpha)$ from $(0.1,0.3)$ to $(0.2,0.1)$. Comparing
these values of $\Delta s$ with those required by observations, we see that
the model could (at least partially) be responsible for the ‘broken power
law’.

The difference of $s(p)$ from the standard value $s=4$ increases toward the
maximum energy, faster for larger $\bar{x}_{\mathrm{pm}}$ (i.e. how far from
the shock the particles can return back), for the higher value of $k$ (i.e.
how fast the flow speed $u_{2}$ grows with $x$) and for the lower $\alpha$
(i.e. how the diffusion coefficient depends on the particle momentum).

The last point leads to an important conclusion: the shape of the accelerated
particle distribution in the momentum space is sensitive to the dependence of
$D$ on $p$ if the flow speed is not spatially uniform. It is clear from
equation (36): the larger the absolute value of $k$, the higher the
sensitivity.

Figure 7: Distribution functions for the same cases as in Fig. 3. The only
difference is the diffusion coefficient. It is the same as in Fig. 3 up to the
momentum $0.1p_{\mathrm{max}}$ and independent of $p$ for higher momenta.
Above $0.1p_{\mathrm{max}}$, the slopes $s$ for the solid blue and red lines
are $2.33$ and $2.16$ respectively. Note that all the lines will be straight
horizontal up to $p_{\mathrm{max}}$ in case of the uniform flow speed
$u_{2}(x)=\mathrm{const}$.

In our model, the diffusion coefficient is a power-law function of the
particle momentum, $D\propto p^{\alpha}$ across the whole range of $p$. As a
result, the momentum spectrum of particles changes smoothly over the wide
range of $p$. The observed spectral feature may reflect the change of
dependence of the diffusion coefficient on $p$ around the break momentum
$p_{\mathrm{br}}$.

Fig. 7 illustrates this concept, showing the same models as Fig. 3 with the
only difference in the diffusion coefficient for particles with the highest
momenta. Namely, it is

$D(p)\propto\left\\{\begin{array}[]{ll}p^{\alpha},&p<p_{\mathrm{br}}\\\
p_{\mathrm{br}}^{\alpha}=\mathrm{const},&p\geq
p_{\mathrm{br}}\end{array}\right.$ (37)

with $p_{\mathrm{br}}=0.1p_{\mathrm{max}}$. We see a sharp break around $p\sim
10^{3}/mc$ in the solid lines representing the Bohm diffusion up to the break.
The slope change is $\Delta s=0.33$ for $x_{\mathrm{pm}}=0.1$ and $\Delta
s=0.16$ for $x_{\mathrm{pm}}=0.05$, meaning that it depends on the depth of
the ‘zone’ downstream which is sampled by the accelerating particles. It is
interesting to note that the slope $s=2.33$ above $p_{\mathrm{br}}$ around the
lower boundary of the values of $s\approx 2.3-2.6$ is required to fit the
observed gamma-ray spectra of some SNRs at the photon energies $\gtrsim
1{\,\rm TeV}$. The index $s$ may even be higher in our model for appropriate
$x_{\mathrm{pm}}$ and $k$.

We would like to emphasize that $D(p)$ influences $s(p)$ only if the flow
speed is non-uniform, in either downstream or upstream direction. If
$u_{2}(x)=\mathrm{const}$ and $u_{1}(x)=\mathrm{const}$, then particles of all
momenta sample the same difference between $u_{1}$ and $u_{2}$ and, therefore,
$s$ is independent of $p$.

###  5.3 Evolution of the diffusion length

The effect we are considering depends on the spatial extent of the downstream
region probed by particles accelerating at the shock. We expect that the
effect from the non-uniform flow velocity downstream of the shock may be
noticeable in real SNRs if the diffusion length $x_{\mathrm{pm}}$ is higher
than say $5\%$ of the radius $R$.

The Larmor radius of a proton or electron with energy $E=E_{\mathrm{TeV}}\cdot
10^{12}{\,\rm eV}$ in the magnetic field $B=B_{\mathrm{\mu}}\cdot
10^{-6}{\,\rm G}$ is
$r_{\mathrm{L}}=10^{-3}E_{\mathrm{TeV}}/B_{\mathrm{\mu}}{\,\rm pc}$. The
length-scale downstream involved into formation of $f_{\mathrm{o}}(p)$ is of
order of the diffusion length $x_{\mathrm{p}}\simeq D_{2}/u_{2}$. If diffusion
is Bohm-like, then $D_{\mathrm{B}}=\eta r_{\mathrm{L}}c/3$ and

$x_{\mathrm{p}}\simeq\frac{4\eta c}{3V}r_{\mathrm{L}}.$ (38)

The diffusion length could be up to $x_{\mathrm{p}}\sim 10^{5}r_{\mathrm{L}}$
(for $\eta=100$ and a slow shock $V=300{\,\rm km/s}$). Under such extreme
circumstances, even the radio-emitting electrons may effectively feel the non-
uniformity of the flow speed distribution because the diffusion length may be
up to $x_{\mathrm{p}}\sim 0.1{\,\rm pc}$ for $10{\,\rm GeV}$ electrons in
$10{\,\rm\mu G}$ magnetic field. Particles with energy
$E_{\mathrm{max}}=10{\,\rm TeV}$ may have $x_{\mathrm{p}}\sim 0.1{\,\rm pc}$
in more common situations, such as in a magnetic field with $50{\,\rm\mu G}$
strength behind the shock with speed $6000{\,\rm km/s}$ and $\eta=10$.

In order to see the sensitivity of the problem to the parameter
$\bar{x}_{\mathrm{pm}}=x_{\mathrm{p}}(p_{\mathrm{max}})/R$, we adopt its fixed
values for calculations in the previous sections. Let us examine how
$\bar{x}_{\mathrm{pm}}$ evolves in the two alternative scenarios: the
efficient acceleration with the Bohm diffusion and the test-particle one with
the Kolmogorov diffusion coefficient.

The turbulent component of the magnetic field $\delta B$ causes the scattering
of particles. The equation (38) yields

$x_{\mathrm{p,B}}\simeq\frac{\eta\sigma_{\mathrm{tot}}c}{3\sigma_{\mathrm{B}}Ve\delta
B_{1}}E=\frac{\eta\sigma_{\mathrm{tot}}}{30\sigma_{\mathrm{B}}}\frac{E_{\mathrm{TeV}}}{V_{\mathrm{3000}}\delta
B_{\mathrm{\mu}}}\ {\,\rm pc}$ (39)

where $V_{\mathrm{3000}}$ is the shock speed in units $3000{\,\rm km/s}$ and
we used $\delta B_{2}=\sigma_{\mathrm{B}}\delta B_{1}$. It follows that the
normalized diffusion length varies over time due to the expansion of SNR, the
evolution of $p_{\mathrm{max}}$, and the magnetic field:
$\bar{x}_{\mathrm{pm}}\propto p_{\mathrm{max}}/(VBR)$. We know the
dependencies $R(t)$ and $V(t)$ for our SNR models from simulations.

Under the efficient acceleration regime, we consider both the ordered $B_{1}$
and the random $\delta B_{1}$ components. At the young shocks, the Bell non-
resonant instability operates with $\delta B_{1}>B_{1}$. The saturated level
of magnetic energy is about the $V/c$ fraction of energy $E_{\mathrm{cr}}$
stored in CRs (Bell, 2004)

$\frac{\delta
B_{1}^{2}}{4\pi}\simeq\frac{V}{c}\xi_{\mathrm{cr}}\rho_{\mathrm{o}}V^{2}$ (40)

where $\xi_{\mathrm{cr}}$ is the fraction of the kinetic energy transferred to
CRs. When this $\delta B_{1}$ is less than the strength of the ordered field
$B_{1}$, the resonant waves are responsible for the particle scatterings with
(Amato & Blasi, 2006)

$\frac{\delta B_{1}^{2}}{8\pi}\simeq\frac{v_{\mathrm{A}}}{V}P_{\mathrm{cr}}$
(41)

where $v_{\mathrm{A}}=B_{1}/\sqrt{4\pi\rho_{\mathrm{o}}}$ and
$P_{\mathrm{cr}}=(\gamma_{\mathrm{cr}}-1)E_{\mathrm{cr}}$ the CR pressure,
$\gamma_{\mathrm{cr}}=4/3$ and $\delta B_{1}/B_{1}\leq 1$. The highest value
of $\delta B_{1}$ for the resonant waves could be $B_{1}$. We assume $\delta
B_{1}=B_{1}$ in situations where these formulae give, at the same time,
simultaneously for the nonresonant $\delta B_{1}<B_{1}$ and the resonant
$\delta B_{1}>B_{1}$. Numerically, for the plasma with the mean mass of
particle in units of the proton mass $\mu=0.609$,

$\delta B_{1}=\left\\{\begin{array}[]{ll}110\
\xi_{\mathrm{cr}}^{1/2}n_{\mathrm{o}}^{1/2}V_{\mathrm{3000}}^{3/2}\ {\,\rm\mu
G},&\delta B_{1}>B_{1}\\\\[8.0pt] 27\
\xi_{\mathrm{cr}}^{1/2}B_{\mathrm{\mu}}^{1/2}n_{\mathrm{o}}^{1/4}V_{\mathrm{3000}}^{1/2}\
{\,\rm\mu G},&\delta B_{1}<B_{1}\\\\[8.0pt]
B_{1}&\mathrm{otherwise}\end{array}\right.$ (42)

Fig. 8 shows the variation of $\delta B_{1}$ with the SNR age in our models.

Figure 8: Evolution of $\delta B_{1}$ in three models of SNRs.

The mean acceleration time $t_{\mathrm{acc}}(p_{\mathrm{max}})$ could be
estimated from a classic integral (eq. 3.31 in Drury, 1983) assuming $D\propto
p^{\alpha}$, $u_{1}=V$, $u_{2}=V/\sigma_{\mathrm{tot}}$:

$t_{\mathrm{acc}}(p_{\mathrm{max}})\simeq\frac{\beta_{\mathrm{a}}}{\alpha}\frac{D(p_{\mathrm{max}})}{V^{2}},\quad\beta_{\mathrm{a}}=\frac{3\sigma_{\mathrm{tot}}(1+\sigma_{\mathrm{tot}}/\sigma_{\mathrm{B}})}{\sigma_{\mathrm{tot}}-1}$
(43)

where $\beta_{\mathrm{a}}$ accounts for the fractions of times a particle
spends upstream and downstream. As expected, the acceleration time is longer
in $1/\alpha$ times compared to the case of the Bohm diffusion ($\alpha=1$),
e.g. in 3 times for the Kolmogorov diffusion ($\alpha=1/3$).

One of three processes determines the maximum momentum of accelerated
particles $p_{\mathrm{max}}$. The condition
$t_{\mathrm{acc}}(p_{\mathrm{max}})=t$ gives the maximim momentum limited by
the SNR age. We have for the Bohm-like diffusion

$p_{\mathrm{max,age,B}}=\frac{3e}{\eta\beta_{\mathrm{a}}}\left(\frac{V}{c}\right)^{2}\delta
B_{1}t.$ (44)

The synchrotron loss time, applicable for electrons only, is

$t_{\mathrm{rad}}=\frac{1}{A_{\mathrm{p}}\langle B^{2}\rangle p},\quad
A_{\mathrm{p}}=\frac{4}{3}\frac{\sigma_{\mathrm{T}}}{8\pi}\frac{1}{(m_{\mathrm{e}}c)^{2}}=4.7\times
10^{7}{\,\rm cgs}$ (45)

with $\langle B^{2}\rangle=\beta_{\mathrm{b}}B_{\mathrm{rad}}^{2}$ where
$B_{\mathrm{rad}}=\max(\delta B_{1},B_{1})$ and $\beta_{\mathrm{b}}$ accounts
for the times particles spend in the upstream and downstream (Reynolds, 1998)

$\beta_{\mathrm{b}}=\frac{\sigma_{\mathrm{B}}+\sigma_{\mathrm{B}}^{2}\sigma_{\mathrm{s}}}{\sigma_{\mathrm{B}}+\sigma_{\mathrm{s}}}.$
(46)

The equation
$t_{\mathrm{acc}}(p_{\mathrm{max}})=t_{\mathrm{rad}}(p_{\mathrm{max}})$
provides

$p_{\mathrm{max,rad,B}}=\left[\frac{3e}{\eta\beta_{\mathrm{a}}\beta_{\mathrm{b}}A_{\mathrm{p}}}\left(\frac{V}{c}\right)^{2}\frac{\delta
B_{1}}{B_{\mathrm{rad}}^{2}}\right]^{1/2}.$ (47)

The third important time to account is the time-to-grow for the Bell
instability that produces the magnetic field (Bell et al., 2013; Cardillo et
al., 2015). In this case, the maximum momentum is estimated by the equation
(5) in Cristofari et al. (2020):

$p_{\mathrm{max,gr}}\simeq\frac{\xi_{\mathrm{cr}}e\pi^{1/2}}{5\Lambda}\frac{4}{4-\omega}R\rho_{\mathrm{o}}^{1/2}\left(\frac{V}{c}\right)^{2}$
(48)

where $e$ is the electron charge, $\omega=-d\ln n_{\mathrm{o}}(R)/d\ln R$ is
the ambient density steepness (it is denoted by $m$ in Cardillo et al., 2015),
$\Lambda\approx 12$ for $p_{\mathrm{max}}\sim 100{\,\rm TeV}$. This expression
is valid for any type of diffusion and when the non-resonant $\delta
B_{1}>B_{1}$.

Table 1: We consider two scenarios for particle acceleration to calculate the evolution of $p_{\mathrm{max}}$ and $\bar{x}_{\mathrm{pm}}$. The type of particles is different to show the role of the radiative losses in one case and to remove it from analysis in another scenario for simplicity. | efficient | test-particle
---|---|---
| acceleration | regime
particle | proton | electron
diffusion | Bohm | Kolmogorov
$\alpha$ | $1$ | $1/3$
$\eta$ | $1$ | $10$
$\xi_{\mathrm{cr}}$ | 0.15 | –
$E_{\mathrm{o}}$ | – | $80{\,\rm GeV}$
$\delta B_{1}$ | (42) | $B_{1}$
$B_{\mathrm{rad}}$ | $\max(\delta B_{1},B_{1})$ | $B_{1}$
$B_{1}$ | $10{\,\rm\mu G}$ | $10{\,\rm\mu G}$
$\sigma_{\mathrm{B}}$ | $\sqrt{11}$ | 1
$\sigma_{\mathrm{tot}}$ | 5.5 | 4
$p_{\mathrm{max}}$ | min (44), (48) | min (50), (51)
$\beta_{\mathrm{a}}$ | 9.7 | 20
$\beta_{\mathrm{b}}$ | 7.2 | 1

Figure 9: Evolution of the maximum momentum $p_{\mathrm{max}}$ (a) and the
normalized diffusion length $\bar{x}_{\mathrm{pm}}$ (b) for particles with the
maximum momentum in the test-particle acceleration scenario with Kolmogorov
diffusion. Different colours show three of our models of SNRs. Line types
represent the time-limited (dashed lines), loss-limited (dot-dashed lines)
$p_{\mathrm{max}}$ and the minimum of the two (solid lines). Lines for
$\bar{x}_{\mathrm{pm}}$ correspond to the solid lines on the plot for the
maximum momenta. $p_{\mathrm{max,rad}}<p_{\mathrm{max,age}}$ from $t=1.3{\,\rm
kyr}$ for the blue and $t=2.3{\,\rm kyr}$ for the green solid lines.

Figure 10: The same as in Fig. 9 for the scenario of efficient acceleration.
Dot-dashed lines represent $p_{\mathrm{max}}$ due to the Bell instability
which operates up to the time when the non-resonant $\delta B_{1}>B_{1}$.

In an alternative scenario, we assume the test-particle acceleration and the
Kolmogorov diffusion. We set the diffusion coefficient the same as the Bohm-
like one at the particle energy $E_{\mathrm{o}}$, i.e. $D_{\mathrm{K}}=\eta
cE_{\mathrm{o}}^{2/3}E^{1/3}/(3e\delta B)$. Then

$x_{\mathrm{p,K}}\simeq\frac{\eta\sigma_{\mathrm{tot}}c}{3\sigma_{\mathrm{B}}Ve\delta
B_{1}}E_{\mathrm{o}}^{2/3}E^{1/3}=\frac{\eta\sigma_{\mathrm{tot}}}{30\sigma_{\mathrm{B}}}\frac{E_{\mathrm{o,TeV}}^{2/3}E_{\mathrm{TeV}}^{1/3}}{V_{\mathrm{3000}}\delta
B_{\mathrm{\mu}}}\ {\,\rm pc}$ (49)
$p_{\mathrm{max,age,K}}=\left[\frac{ec^{2/3}}{\eta\beta_{\mathrm{a}}E_{\mathrm{o}}^{2/3}}\left(\frac{V}{c}\right)^{2}B_{1}t\right]^{3}.$
(50)
$p_{\mathrm{max,rad,K}}=\left[\frac{ec^{2/3}}{\eta\beta_{\mathrm{a}}\beta_{\mathrm{b}}A_{\mathrm{p}}E_{\mathrm{o}}^{2/3}}\left(\frac{V}{c}\right)^{2}\frac{\delta
B_{1}}{B_{\mathrm{rad}}^{2}}\right]^{3/4}.$ (51)

‘Test-particle’ relativistic protons do not amplify the magnetic field.
Therefore, we do not account for $p_{\mathrm{max,gr}}$ in this scenario. MF
strength for the diffusion coefficient and the radiative losses is the ordered
one at the parallel shock $\delta B_{1}=B_{\mathrm{rad}}=B_{1}=B_{2}$.

The two scenarios for the evolution of the diffusion length
$\bar{x}_{\mathrm{pm}}=x_{\mathrm{pm}}/R$ are summarized in Table 1. Figs. 9,
10 show results of calculations.

In the test-particle case, the maximum momentum in our models is limited
predominantly by the age of the system (Fig. 9a). Interestingly, the three
lines for $\bar{x}_{\mathrm{pm}}$ (Fig. 9b) have shapes similar to the
evolution of the deceleration parameter $m=Vt/R$ (figure 2 in Petruk et al.,
2021). Indeed, if we express $D$ from the relations $t_{\mathrm{acc}}=t$ or
$t_{\mathrm{acc}}=t_{\mathrm{rad}}$ and substitute
$\bar{x}_{\mathrm{pm}}=\sigma_{\mathrm{tot}}D(p_{\mathrm{max}})/(RV)$ with it,
we get

$\bar{x}_{\mathrm{pm}}=\left\\{\begin{array}[]{ll}\displaystyle\frac{\alpha\sigma_{\mathrm{tot}}}{\beta_{\mathrm{a}}}m(t),&t\geq
t_{\mathrm{rad}}\\\\[8.0pt]
\displaystyle\frac{\alpha\sigma_{\mathrm{tot}}}{\beta_{\mathrm{a}}}\frac{t_{\mathrm{rad}}(p_{\mathrm{max}})}{t}m(t),&t<t_{\mathrm{rad}}\\\
\end{array}\right.$ (52)

that is valid for diffusion of any type. If the age of the systems limits
$p_{\mathrm{max}}$, then the correlation $\bar{x}_{\mathrm{pm}}(t)\propto
m(t)$ is clear. When the radiative losses limit the maximum energy of
particles, then $\bar{x}_{\mathrm{pm}}$ scales also with parameters
determining $t_{\mathrm{rad}}$ (Fig. 9b, blue line after $t=1.3{\,\rm kyr}$,
green line after $t=2.3{\,\rm kyr}$). It could be essential in real SNRs when
radiative losses limit $p_{\mathrm{max}}$, e.g. when MF is amplified to high
values. Then $\bar{x}_{\mathrm{pm}}\propto B^{-1/2}t^{8/5}$ for a Sedov SNR as
an example, that might provide a way to determine properties of magnetic field
from the non-thermal spectra of SNRs. When $p_{\mathrm{max}}$ is age-limited,
the normalized diffusion length for particles with maximum momentum is
determined solely by the shock compression and diffusion type (represented by
$\alpha$). In particular, for the Bohm-like (Kolmogorov) diffusion and
unmodified shock, $\bar{x}_{\mathrm{pm}}\approx 20\ (6)\%$ of the SNR radius
in the young SNRs (when $m\approx 1$) or $\bar{x}_{\mathrm{pm}}\approx 8\
(3)\%$ at the Sedov stage (when $m=0.4$).

In the scenario of efficient acceleration, $p_{\mathrm{max}}$ is limited by
the grow-rate of the Bell instability up to $t\approx 3{\,\rm kyr}$ (Fig. 10a,
cf. Fig. 8). Eqs. (39), (40) and (48) yield then

$\bar{x}_{\mathrm{pm}}=0.025\
\frac{\eta\sigma_{\mathrm{tot}}}{\sigma_{\mathrm{B}}}\frac{4}{4-\omega}\frac{\xi_{\mathrm{cr}}^{1/2}}{V_{\mathrm{3000}}^{1/2}}\sim
0.02\eta V_{\mathrm{3000}}^{-1/2}.$ (53)

The dependence on the shock speed is not strong, and the evolution of
$\bar{x}_{\mathrm{pm}}$ is slow. Indeed, we see in Fig. 10b that the region
downstream involved into the acceleration process is about $1\div 3\%$ of the
radius when the Bell instability is effective. However, if the diffusion is
Bohm-like with $\eta=2\div 3$, then $\bar{x}_{\mathrm{pm}}$ could reach $5\%$
of the SNR radius. At times after $\sim 3{\,\rm kyr}$, when the non-resonant
instability is not effective in scattering the particles, then the SNR age
limits the maximum momentum. At this stage, the diffusion length does not
depend on $\eta$ (Eqs. 41 and 44) and may be up to $\sim 0.1R$.

Thus, the downstream structure of the flow might have some influence on the
shape of the non-thermal spectra of SNRs in both scenarios we are considering.
The effect is not strong enough but may modify the spectral indices in
realistic cases under some circumstances.

###  5.4 Limitations

The non-uniformity of the plasma speed behind the shock affects the cosmic-ray
spectrum to a larger extent if particles can penetrate deeply downstream while
accelerating at the shock. If so, they can probe a larger (for decreasing
$u_{2}(x)$) difference between the upstream $u_{\mathrm{p1}}$ and downstream
$u_{\mathrm{p2}}$ flow speeds and get a higher bust in momentum $\Delta
p\propto u_{\mathrm{p1}}-u_{\mathrm{p2}}$. The distance from the shock depends
on the properties of diffusion, $x_{\mathrm{p2}}\simeq D/u_{2}$. Therefore, in
order the effect from the downstream gradient of $u_{2}$ be visible in the
spectrum, the diffusion coefficient should be large. In our estimates we
considered $x_{\mathrm{p1}}$ to be $5\div 10\%$ of the SNR radius $R$ for the
particles with the maximum momentum
$p_{\mathrm{max}}=10^{4}m_{\mathrm{p}}c\approx 10{\,\rm TeV}/c$. This
corresponds to a diffusion coefficient for particles with such momentum $D\sim
0.1Ru_{2}\sim 3\times 10^{25}{\,\rm cm^{2}/s}$ (we used $R\sim 1{\,\rm pc}$
and $u_{2}\sim 10^{3}{\,\rm km/s}$ for this estimate), that is well within
typical values used in the literature. On the contrary, to be capable of
providing the proton acceleration to energies beyond $100{\,\rm TeV}$, the
diffusion coefficient should be the smallest one, namely, the Bohm one
$D=r_{\mathrm{L}}c/3$ in the magnetic field of few $100{\,\rm\mu G}$, that
requires substantial magnetic field amplification. Under such conditions, the
diffusion length for particles with energy $E_{\mathrm{max}}=10{\,\rm TeV}$ is
$D_{2}/u_{2}\simeq r_{\mathrm{L}}c/u_{2}\sim 0.01{\,\rm pc}$ that is just
$\sim 1\%$ of $R$ for particles with the maximum momentum. Of course, the
variation of the plasma speed provided by the SNR hydrodynamics on such length
scales downstream does not affect the shape of the CR spectrum. On the other
hand, the length-scale $x_{\mathrm{p}}\propto D\propto
p_{\mathrm{max}}^{\alpha}$ also increases with rise of the maximum momentum.

In general, we expect that the effect of the downstream velocity gradient
could be prominent in SNRs or their regions where the particle acceleration is
not very efficient, i.e. with a lower magnetic field and higher diffusion
coefficient. Thus, in the present paper, we are not in the conditions to
provide the highest $E_{\mathrm{max}}\sim{\,\rm PeV}$ energies. Instead, the
maximum energy $E_{\mathrm{max}}\approx 10{\,\rm TeV}$ we used in our examples
is a value which may be reached in SNRs with a common interstellar magnetic
field and without the fastest diffusion regime. It seems that the conditions
for the very efficient acceleration (amplified $B\sim\mathrm{few}{\,\rm mG}$,
the smallest Bohm diffusion coefficient), while crucial for the investigation
of PeVatrons, are not widespread in SNRs. Such conditions could be expected
locally, i.e. present only in some small regions of the SNR shock surface. As
illustration of this point, we refer the reader to the figure 4 in Long et al.
(2003). The thinnest radial profile of the X-ray surface brightness in region
A of SN1006 is due to the very fast radiative cooling of electrons in the high
magnetic field. This profile provides evidence of effective acceleration with
magnetic field amplification there. However, two other profiles B and C, being
even nearby, are thicker (i.e. magnetic field is lower in these regions).
Therefore, the high magnetic strength inferred from the thinnest profile may
not be generally attributed to the whole SNR. The less efficient particle
acceleration could be in operation throughout other regions of the SNR shock
surface.

Another point is related to the dependence of the diffusion coefficient on the
particle momentum, i.e. the value of the index $\alpha$. The effect from the
gradient of $u_{2}$ could be more prominent at the smaller particle momenta if
the index $\alpha$ is less than unity (like in the Kolmogorov or Kraichnan
diffusion coefficients) because $x_{\mathrm{p}}(p)\propto
x_{\mathrm{p}}(p_{\mathrm{max}})p^{\alpha}$ (see Figs. 3, 4, 7).

Interestingly, numerical studies of the non-linear particle acceleration
(where the kinetic particle equation was solved, along with the system of the
hydrodynamic equations for the flow) by different groups (Kang & Jones, 1995,
2007; Gieseler et al., 2000; Kang et al., 2002; Berezhko & Völk, 1997, 2000;
Zirakashvili & Aharonian, 2010; Zirakashvili & Ptuskin, 2012) used either the
Bohm or the Bohm-like diffusion coefficient and did not find an appreciable
impact of the downstream velocity gradient on the particle spectrum.

We found one report where the effect appears visible (though it is difficult
to separate from the eventual influence of other factors). Brose et al. (2020)
numerically solved a system of equations for cosmic rays, magnetic turbulence
and plasma hydrodynamics in SNRs. Figure 2 in their calculations compares the
particle spectra for the Bohm diffusion (red lines) and the diffusion
coefficient calculated self-consistently from the magnetic turbulence (green
lines). In the latter case, the coefficient becomes Kolmogorov-like for most
energies (Brose 2023, private communication). The differences between the
spectra shapes for $\alpha=1$ and $\alpha=1/3$ at the ages $10$ and $60{\,\rm
kyr}$ are similar to ours in Fig. 3.

In summary, the effect of the non-uniform $u_{2}(x)$ on the particle spectrum
could be prominent in SNRs where adoption of Kolmogorov ($\alpha=1/3$) or
Kraichnan diffusion ($\alpha=1/2$) in magnetic fields with standard strength
$B$ of few $10{\,\rm\mu G}$ is acceptable, or, in case of a more efficient
acceleration with the Bohm-like diffusion coefficient with $\eta>1$. In
contrast, such values of $\alpha$, $B$ and $\eta$ are not compatible with the
highest maximum energies of cosmic rays observed in several SNRs and are not
consistent with a modern view of anisotropic cascade of MHD turbulence in the
sources of highly efficient particle acceleration based on the Bohm diffusion
and high values of amplified $B$.

##  6 Conclusions

Our research aims to study the possible effect of the spatially variable
distribution of the flow speed downstream of the shock on the momentum
distribution of particles accelerated at the shock.

To achieve this, we generalized the approach of Bell (1978b) which is an
alternative way to describe the particle acceleration at the shock. This
result is of general academic interest of its own.

We applied this generalized individual particle approach to the NLA solution
of Blasi (2002), which assumes $u_{2}(x)=\mathrm{const}$ and then to a more
general problem with $u_{2}(x)\neq\mathrm{const}$.

In order to cross-check our results, we derived in Appendix C the same general
solution for $f_{\mathrm{o}}(p)$ by solving the kinetic equation.

The formula for the general solution (23) is similar to the known NLA
solution. The difference consists in substitution $u_{\mathrm{p2}}(p)$ instead
of the constant $u_{2}$.

In Sect. 3, we present profiles of the flow speed derived in numerical MHD
simulations of remnants evolving after supernova explosions of different
types. We demonstrate that the flow speed may vary prominently on the spatial
scales involved in the acceleration of particles around the SNR shocks.

Then, we calculate the momentum spectra of particles for several practical
cases to demonstrate the effect and its limitations. In these applications, we
treat the problem by taking the formula (35) for the function
$u_{\mathrm{p2}}(p)$. Such a simplified description of $u_{\mathrm{p2}}(p)$
allows us to study the effect of a non-uniformity of the flow speed on the
particle spectrum accelerated in the test-particle limits or the NLA regime
(Sect. 4.3).

In Sect. 5, we suggest two cases where the phenomenon might affect, at least
partially, some observed properties of SNRs.

The overall effect may be visible from the equation (36) written for the test-
particle limit. If the flow speed in the frame of the shock rises with
distance from the shock (i.e. $k>0$ like in evolved SNRs) and particles with
higher $p$ reach larger $x$, then $s$ grows with $p$. This results to a convex
spectrum of particles. If the flow speed declines ($k<0$), then $s$ decreases
and the spectrum is concave. The typical assumption of the uniform flow
downstream ($k=0$) results in the classic formulae.

An important conclusion apparent from the equation (36) is that the diffusion
coefficient (represented by $\alpha$ in this formula) directly affects the
shape of $f_{\mathrm{o}}(p)$ only if the flow speed is not uniform upstream
and/or downstream. Diffusion properties of accelerating particles are
important in shaping the $f_{\mathrm{o}}(p)$ because the diffusion coefficient
determines the diffusive lengths for particles with different momenta. The
momentum increase of a particle is sensitive to the difference of the flow
velocities it probes
$u_{\mathrm{1}}(x_{\mathrm{p1}})-u_{\mathrm{2}}(x_{\mathrm{p2}})$. Therefore,
the same speed profile $u_{2}(x)$ affects the spectrum $f_{\mathrm{p}}(p)$ in
a different way in models with distinct diffusion coefficients.

The effect of interest is rather sensitive to the diffusion properties of
cosmic rays. Particles should diffuse deeply downstream to probe a
hydrodynamic variation of $u_{2}(x)$. Therefore, their diffusion length and
diffusion coefficient should be high. As the acceleration of cosmic rays in
SNRs to energies beyond $100{\,\rm TeV}$ is likely achievable only through the
Bohm diffusion mechanism in an amplified magnetic field (then the diffusion
length is small), the significance of the downstream gradient of plasma
velocity in SNRs, capable of generating PeV cosmic rays, may generally be
disregarded.

###### Acknowledgements.

We acknowledge Elena Amato for useful discussions. OP acknowledges the OAPa
grant number D.D.75/2022 funded by Direzione Scientifica of Istituto Nazionale
di Astrofisica, Italy. The study was partially supported by the INAF mini-
grant 1.05.23.04.04. This project has received funding through the
MSCA4Ukraine project, which is funded by the European Union. Views and
opinions expressed are however those of the author(s) only and do not
necessarily reflect those of the European Union. Neither the European Union
nor the MSCA4Ukraine Consortium as a whole nor any individual member
institutions of the MSCA4Ukraine Consortium can be held responsible for them.
We thank the Armed Forces of Ukraine for providing security to perform this
work. This research used an HPC cluster at DESY, Germany, and the HPC system
MEUSA at INAF-Osservatorio Astronomico di Palermo, Italy.

## References

  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807
  * Amato & Blasi (2005) Amato, E. & Blasi, P. 2005, MNRAS, 364, L76
  * Amato & Blasi (2006) Amato, E. & Blasi, P. 2006, MNRAS, 371, 1251
  * Axford et al. (1977) Axford, W. I., Leer, E., & Skadron, G. 1977, in International Cosmic Ray Conference, Vol. 11, International Cosmic Ray Conference, 132
  * Bell (1978a) Bell, A. R. 1978a, MNRAS, 182, 147
  * Bell (1978b) Bell, A. R. 1978b, MNRAS, 182, 443
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415
  * Berezhko & Völk (1997) Berezhko, E. G. & Völk, H. J. 1997, Astroparticle Physics, 7, 183
  * Berezhko & Völk (2000) Berezhko, E. G. & Völk, H. J. 2000, A&A, 357, 283
  * Blandford & Ostriker (1978) Blandford, R. D. & Ostriker, J. P. 1978, ApJ, 221, L29
  * Blasi (2002) Blasi, P. 2002, Astroparticle Physics, 16, 429
  * Blasi et al. (2005) Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907
  * Bogdan & Volk (1983) Bogdan, T. J. & Volk, H. J. 1983, A&A, 122, 129
  * Brose et al. (2020) Brose, R., Pohl, M., Sushch, I., Petruk, O., & Kuzyo, T. 2020, A&A, 634, A59
  * Caprioli et al. (2010) Caprioli, D., Amato, E., & Blasi, P. 2010, Astroparticle Physics, 33, 307
  * Caprioli et al. (2009) Caprioli, D., Blasi, P., Amato, E., & Vietri, M. 2009, MNRAS, 395, 895
  * Cardillo et al. (2015) Cardillo, M., Amato, E., & Blasi, P. 2015, Astroparticle Physics, 69, 1
  * Chevalier (1982) Chevalier, R. A. 1982, ApJ, 258, 790
  * Cristofari et al. (2020) Cristofari, P., Blasi, P., & Amato, E. 2020, Astroparticle Physics, 123, 102492
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973
  * Dwarkadas & Chevalier (1998) Dwarkadas, V. V. & Chevalier, R. A. 1998, ApJ, 497, 807
  * Ferrand et al. (2019) Ferrand, G., Warren, D. C., Ono, M., et al. 2019, ApJ, 877, 136
  * Ferrand et al. (2021) Ferrand, G., Warren, D. C., Ono, M., et al. 2021, ApJ, 906, 93
  * Gabler et al. (2021) Gabler, M., Wongwathanarat, A., & Janka, H.-T. 2021, MNRAS, 502, 3264
  * Gieseler et al. (2000) Gieseler, U. D. J., Jones, T. W., & Kang, H. 2000, A&A, 364, 911
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A5
  * H. E. S. S. Collaboration et al. (2015) H. E. S. S. Collaboration, Abramowski, A., Aharonian, F., et al. 2015, A&A, 574, A100
  * Jones & Ellison (1991) Jones, F. C. & Ellison, D. C. 1991, Space Sci. Rev., 58, 259
  * Kang & Jones (1995) Kang, H. & Jones, T. W. 1995, ApJ, 447, 944
  * Kang & Jones (2007) Kang, H. & Jones, T. W. 2007, Astroparticle Physics, 28, 232
  * Kang et al. (2002) Kang, H., Jones, T. W., & Gieseler, U. D. J. 2002, ApJ, 579, 337
  * Krymskii (1977) Krymskii, G. 1977, Akademiia Nauk SSSR Doklady, 234, 1306
  * Long et al. (2003) Long, K. S., Reynolds, S. P., Raymond, J. C., et al. 2003, ApJ, 586, 1162
  * Loru et al. (2019) Loru, S., Pellizzoni, A., Egron, E., et al. 2019, MNRAS, 482, 3857
  * Malkov & Drury (2001) Malkov, M. A. & Drury, L. O. 2001, Reports on Progress in Physics, 64, 429
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228
  * Mignone et al. (2012) Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7
  * Orlando et al. (2015) Orlando, S., Miceli, M., Pumo, M. L., & Bocchino, F. 2015, ApJ, 810, 168
  * Patnaude et al. (2015) Patnaude, D. J., Lee, S.-H., Slane, P. O., et al. 2015, ApJ, 803, 101
  * Petruk (2000) Petruk, O. 2000, A&A, 357, 686
  * Petruk et al. (2021) Petruk, O., Kuzyo, T., Orlando, S., Pohl, M., & Brose, R. 2021, MNRAS, 505, 755
  * Planck Collaboration et al. (2016) Planck Collaboration, Arnaud, M., Ashdown, M., et al. 2016, A&A, 586, A134
  * Reynolds (1998) Reynolds, S. P. 1998, ApJ, 493, 375
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic Press Inc., New York)
  * Taylor (1950) Taylor, G. 1950, Proceedings of the Royal Society of London Series A, 201, 159
  * Telezhinsky et al. (2012) Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012, Astroparticle Physics, 35, 300
  * Telezhinsky et al. (2013) Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2013, A&A, 552, A102
  * Zirakashvili & Aharonian (2010) Zirakashvili, V. N. & Aharonian, F. A. 2010, ApJ, 708, 965
  * Zirakashvili & Ptuskin (2012) Zirakashvili, V. N. & Ptuskin, V. S. 2012, Astroparticle Physics, 39, 12

##  Appendix A Evolution of the flow velocity downstream of the forward shock.
Animation

The animation shows the time development of the profiles of the flow speed
$v(r)$ in the observer reference frame (left) and $u(x)$ in the shock
reference frame (right) for the same models of SNRs as described in Sect. 3
and shown in Fig. 1, up to $10\,000{\,\rm yrs}$. The horizontal axes represent
$r$ on the right plot and $x=1-r/R$ on the left plot. Vertical axes are
normalized on $v_{\mathrm{2}}(0)$ (left) and $u_{1}(0)$ (right). The dashed
line on each plot corresponds to the flow speed in the Sedov (1959) solution.

##  Appendix B Probability of return

The probability of return $P$ in one acceleration cycle is given by the ratio
$P={F_{-x}}/{F_{+x}}$ where $F_{-x}$ and $F_{+x}$ are the fluxes of particles
with momentum $p$ in the negative and positive $x$-directions (Jones & Ellison
1991). In our case, the flow speed varies with $x$ and we need to integrate
over all places where the probability ${\cal X}_{2}(p,x)$ to find particles
with this momentum is above zero, i.e.

$\begin{array}[]{ll}F_{-x}&=\displaystyle\int\limits_{0}^{1}d{\cal
X}_{2}\left|\int\limits_{-v}^{-u_{2}(x)}dv_{x}(u_{2}(x)+v_{x})\right|\\\
&\displaystyle=\frac{v^{2}}{2}\int\limits_{+\infty}^{+0}\left(1-\frac{u_{2}(x)}{v}\right)^{2}\frac{d{\cal
X}_{2}}{dx}dx\\\ F_{+x}&=\displaystyle\int\limits_{0}^{1}d{\cal
X}_{2}\int\limits_{-u_{2}(x)}^{v}dv_{x}(u_{2}(x)+v_{x})\\\
&\displaystyle=\frac{v^{2}}{2}\int\limits_{+\infty}^{+0}\left(1+\frac{u_{2}(x)}{v}\right)^{2}\frac{d{\cal
X}_{2}}{dx}dx.\end{array}$ (54)

Transition to the rightmost equalities assumes that ${\cal X}_{2}=1$ at $x=+0$
and ${\cal X}_{2}=0$ for $x=+\infty$. The simplified treatment of the
probability in Sect. 4.1 may be recovered by taking approximately

${\cal X}_{2}(x)={\cal H}(x_{\mathrm{p}}-x)$ (55)

and considering that $d{\cal
H}(x_{\mathrm{p}}-x)/dx=-\delta(x-x_{\mathrm{p}})$. Then

$P=\displaystyle\left(\frac{1-u_{\mathrm{p2}}/v}{1+u_{\mathrm{p2}}/v}\right)^{2}$
(56)

where $u_{\mathrm{p2}}=u_{2}(x_{\mathrm{p}})$ as it follows from the
approximation (55) and the definition

$u_{\mathrm{p2}}\equiv\int_{0}^{1}u_{2}(x)d{\cal
X}_{2}=u_{2}(0)-\int_{+\infty}^{+0}{\cal X}_{2}(x)\frac{du_{2}}{dx}dx.$ (57)

In a more general treatment, we neglect the terms of the order $(u_{2}/v)^{2}$
in (54). Then, $(1\pm u_{2}/v)^{2}\approx 1\pm 2u_{2}/v$ and, with the use of
the definition (57), the probability of return in one cycle becomes

$P=\frac{1-{2u_{\mathrm{p2}}}/{v}}{1+{2u_{\mathrm{p2}}}/{v}}.$ (58)

In many cycles, we have the product of the probabilities for individual
cycles. Decomposition of the logarithm of the product into the series to the
first order in $u_{\mathrm{p2}}/v$ yields the equation (21).

##  Appendix C Solution of the kinetic equation

This Appendix presents an alternative approach to derive the expression (23)
from the kinetic equation.

We consider the standard steady-state equation for the particle distribution
function $f(x,p)$:

$u\frac{\partial{f}}{\partial{x}}=\frac{\partial{}}{\partial{x}}\left[D\frac{\partial{f}}{\partial{x}}\right]+\frac{1}{3}\frac{du}{dx}p\frac{\partial{f}}{\partial{p}}+Q_{\mathrm{o}}(p)\delta(x)$
(59)

together with the typical conditions: the flow speed is along the $x$ axis, in
the positive $x$ direction; the shock is located at the coordinate $x=0$; the
velocity jump $du/dx=(u_{2}-u_{1})\delta(x)$ through the point $x=0$;
$f(-0,p)=f(+0,p)$; $f(x,p)=0$ at $x=\pm\infty$ and $df/dx=0$ at $x=-\infty$;
the momentum part of the injection term

$Q_{\mathrm{o}}(p)=\displaystyle\frac{\eta n_{1}u_{1}}{4\pi
p^{2}}\delta(p-p_{\mathrm{o}})$ (60)

where $p_{\mathrm{o}}$ is the injection momentum,$\eta$ is the injection
efficiency (a fraction of particles which start acceleration),
$p_{\mathrm{o}}$ the injection momentum, $\delta(x)$ is the delta-function,
the indices ‘1’ and ‘2’ mark upstream and downstream locations.

We consider the general situation when the flow speed profiles are not
spatially constant: $u_{1}(x)\neq\mathrm{const}$,
$u_{2}(x)\neq\mathrm{const}$. The approach to solving the equation (59) is
similar to solving the non-linear problem (Blasi 2002; Blasi et al. 2005),
i.e. we have to account for the spatial variations of the flow velocity.

The integration of the main equation (59) from $-\infty$ to $x=-0$ results in

$\left[D\frac{\partial{f}}{\partial{x}}\right]_{1}=u_{\mathrm{p1}}f_{\mathrm{o}}-\frac{p}{3}\frac{\partial{}}{\partial{p}}\int_{-\infty}^{-0}f\frac{du}{dx}dx,$
(61)

where $f_{\mathrm{o}}(p)$ is the distribution function at the shock. The
integration of the equation (59) from $x=-0$ to $x=+0$ and from $x=+0$ to
$+\infty$ yields

$\left[D\frac{\partial{f}}{\partial{x}}\right]_{2}-\left[D\frac{\partial{f}}{\partial{x}}\right]_{1}+\frac{p}{3}\frac{\partial{f_{\mathrm{o}}}}{\partial{p}}(u_{2}-u_{1})+Q_{\mathrm{o}}(p)=0$
(62)
$\left[D\frac{\partial{f}}{\partial{x}}\right]_{2}=\frac{p}{3}\frac{\partial{}}{\partial{p}}\int_{+0}^{+\infty}f\frac{du}{dx}dx-\int_{+0}^{+\infty}u\frac{\partial{f}}{\partial{x}}dx-
u_{\mathrm{x}}f_{\mathrm{o}}.$ (63)

Integration of the last integral in (63) by parts transforms it to

$-\int_{+0}^{+\infty}u\frac{\partial{f}}{\partial{x}}dx=u_{\mathrm{p2}}f_{\mathrm{o}}.$
(64)

We introduced the notations

$\left[D\frac{\partial{f}}{\partial{x}}\right]_{+\infty}=-u_{3}f_{3}\equiv-
u_{\mathrm{x}}f_{\mathrm{o}}$ (65)

where $u_{3}$ and $f_{3}(p)$ are at the downstream infinity, $u_{\mathrm{x}}$
an unknown function, as well as (cf. with equations (14)-(15))

$u_{\mathrm{p1}}(p)=u_{1}-\int_{-\infty}^{-0}\frac{f(x,p)}{f_{\mathrm{o}}(p)}\frac{du}{dx}dx,$
(66)
$u_{\mathrm{p2}}(p)=u_{2}-\int_{+\infty}^{+0}\frac{f(x,p)}{f_{\mathrm{o}}(p)}\frac{du}{dx}dx$
(67)

With all these, the equation (62) may be re-written as the equation for the
distribution function at the shock $f_{\mathrm{o}}(p)$:

$\frac{\partial{f_{\mathrm{o}}}}{\partial{p}}+\frac{f_{\mathrm{o}}}{p}s_{\mathrm{o}}(p)=\frac{3Q_{\mathrm{o}}(p)}{p(u_{\mathrm{p1}}-u_{\mathrm{p2}})},$
(68)
$s_{\mathrm{o}}=\frac{3}{u_{\mathrm{p1}}-u_{\mathrm{p2}}}\left(u_{\mathrm{p1}}-u_{\mathrm{p2}}+u_{\mathrm{x}}+\frac{p}{3}\frac{\partial{(u_{\mathrm{p1}}-u_{\mathrm{p2}})}}{\partial{p}}\right).$
(69)

The solution of this equation is

$f_{\mathrm{o}}(p)=\frac{\eta n_{1}}{4\pi
p_{\mathrm{o}}^{3}}\frac{3u_{1}}{u_{\mathrm{p1}}-u_{\mathrm{p2}}}\exp\left[-\int_{p_{\mathrm{o}}}^{p}\frac{3(u_{\mathrm{p1}}-u_{\mathrm{p2}}+u_{\mathrm{x}})}{u_{\mathrm{p1}}-u_{\mathrm{p2}}}\frac{dp^{\prime}}{p^{\prime}}\right].$
(70)

The unknown function $u_{\mathrm{x}}$ should be
$u_{\mathrm{x}}=u_{\mathrm{p2}}$ in order (i) to match the expression (23)
derived in the frame of the approach of Bell, (ii) to give the Blasi (2002)
solution with $u_{2}=\mathrm{const}$ (cf. eq.(8) in Blasi et al. 2005), (iii)
to yield the classic test-particle solution $f_{\mathrm{o}}\propto
p^{-3u_{1}/(u_{1}-u_{2})}$ with $u_{1}=\mathrm{const}$ and
$u_{2}=\mathrm{const}$ (then $u_{\mathrm{p1}}=u_{1}$ and
$u_{\mathrm{p2}}=u_{2}$).

Strictly speaking, $u_{\mathrm{x}}$ in the solution (70) is
$u_{\mathrm{x}}=u_{3}f_{3}(p)/f_{\mathrm{o}}(p)$ (see eq. 65). The equality
$u_{\mathrm{x}}=u_{\mathrm{p2}}$ is valid in the case (i) above at least to
the level of approximation for the spatial ‘part’ of the distribution function
$f_{2}(x,p)$ adopted in the present paper, i.e. $f_{2}(x,p)\approx
f_{\mathrm{o}}(p){\cal H}(x_{\mathrm{p}}-x)$. Indeed, let us consider some
point $x_{\mathrm{p^{\prime}}}$ downstream. The advected flux of particles
with momentum $p^{\prime}$ from this point toward $x>x_{\mathrm{p^{\prime}}}$
is proportional to $u_{2}(x_{\mathrm{p^{\prime}}})f_{\mathrm{o}}(p^{\prime})$.
The flux conservation gives
$u_{3}f_{3}(p^{\prime})=u_{2}(x_{\mathrm{p^{\prime}}})f_{\mathrm{o}}(p^{\prime})$.
Similarly,
$u_{3}f_{3}(p^{\prime\prime})=u_{2}(x_{\mathrm{p^{\prime\prime}}})f_{\mathrm{o}}(p^{\prime\prime})$
for particles with another momentum $p^{\prime\prime}$, and so on. Therefore,
$u_{3}f_{3}(p)=u_{2}(x_{\mathrm{p}})f_{\mathrm{o}}(p)$ holds for any momentum
$p$ of relativistic particles. From the definition (15) and ${\cal
X}_{2}(x)={\cal H}(x_{\mathrm{p}}-x)$, we have that
$u_{2}(x_{\mathrm{p}})=u_{\mathrm{p2}}$. Thus,
$u_{3}f_{3}(p)/f_{\mathrm{o}}(p)=u_{\mathrm{p2}}$ and, from (65),
$u_{\mathrm{x}}=u_{\mathrm{p2}}$.

It would be desirable in the future to determine the distribution $f_{3}(p)$
from more general considerations. This task is related to the problem of
distinguishing between the distribution functions of two populations of
particles, namely, those which continue to participate in acceleration and
those advected downstream. Nevertheless, the exact evolution of particles
advected toward the downstream infinity seems unimportant to determine the
distribution function at the shock because they are ‘removed’ from the
accelerator. Therefore, one may set up the boundary condition (65) as
$u_{3}f_{3}(p)=u_{\mathrm{p2}}f_{\mathrm{o}}$ a priori.