© The Authors 2026

1 Introduction

Woltjer (1972) originally subdivided the evolution of supernova remnants (SNRs) into three phases, followed by a merging with the environment. Phase I lasts until the swept-up mass becomes comparable to the ejecta mass. Phase II then keeps going up to an age comparable to a timescale for the radiative losses close to the shock front. Phase III proceeds with the dynamics dominated by the radiative losses. It then appeared that this picture should be supplemented by two intermediate phases (between I and II, and between II and III), which have their own distinctive physical features.

In an updated scheme, but still idealized (e.g., valid in the case of a homogeneous ambient medium), the adiabatic shock in a SNR evolves through three stages: (1) ejecta-dominated or free expansion (Chevalier 1982a,b), (2) energy conversion, from mostly kinetic toward a fixed balance with the thermal (Petruk et al. 2021; Tang & Chevalier 2017), (3) frequently called “adiabatic” or Sedov-Taylor, with a constant ratio of kinetic to thermal energy (Sedov 1946; Taylor 1950). After the radiative losses of plasma become important, a radiative shock develops further through (Cioffi et al. 1988; Slavin & Cox 1992; Blondin et al. 1998; Petruk et al. 2016, 2018) (4) a partially radiative or post-adiabatic stage, and then (5) a fully radiative one, which is concluded by (6) a dissipation phase. The structure of the flow undergoes essential transformations through the fourth stage, while it remains as simple as a “snowplow“ over the fifth one.

We often use the words stage and phase as synonyms. The term era, however, implies a broader segment of evolution. For instance, the adiabatic era encompasses the entire evolution of the SNR from the explosion up to the point where radiative losses begin to affect the dynamics (i.e., stages 1–3). From that time onward the radiative era extends until the complete dissipation of the SNR into the surrounding environment (stages 4–6).

In the present paper, we focus on the radiative era and, more precisely, on its initial stage: the partially radiative stage. The transition of an SNR into this stage has been studied in a number of papers (see Petruk 2005, for a review). Cioffi et al. (1988); Blondin et al. (1998) focused on the hydrodynamics (HD) of the SNRs, while Slavin & Cox (1992); Petruk et al. (2016, 2018) studied their magneto-hydrodynamical (MHD) properties. The beginning and end of the post-adiabatic phase (stage 4 in the sequence listed above) are marked approximately by the transition time t_tr and the shell formation time t_sf. The first, t_tr, points to the end of the adiabatic regime, namely when the radiative losses become prominent, and the expansion parameter m (R ∝ t^m) deviates from the value 2/5, which remains constant during the Sedov stage. The second time, t_sf, corresponds to the beginning of the formation of a thin fully radiative shell outlining the SNR. The adiabatic regime of expansion is valid up to t_tr, while the fully radiative regime is valid after t_sf. An analytical solution for the radiative shock (Pasko & Silich 1986; Bandiera & Petruk 2004) can be used after t_sf.

One of the important conclusions from the HD studies is that the shock undergoing radiative losses is more compressible than the adiabatic shock. In the purely hydrodynamic simulations, the post-shock density can rise to hundred times the pre-shock value. Such a large increase in density can lead to prominent nonthermal emission, due to an increase in injected CR density, even if the efficiency of cosmic ray (CR) acceleration decreases (Lee et al. 2015; Diesing et al. 2024; Diesing & Gupta 2025). However, the tangential component of the magnetic field (MF) can significantly alter the spatial structure of the flow downstream and radically limit the density increase downstream (Petruk et al. 2016, 2018). This happens because the tangential MF is compressed around the partially radiative shock, together with the matter, which causes the magnetic pressure to grow to a dynamically important level. As a consequence, the magnetic pressure essentially limits the compression of the radiative plasma compared to the case when MF is negligible. In this way, the nonthermal evidence of a radiative shell, as predicted by purely hydrodynamic considerations, will vanish (Diesing & Gupta 2025). In contrast, the plasma compression is not influenced by the radial MF component.

In addition to the increased compression, the post-adiabatic flow features another important property. There are rapid and substantial changes in the spatial structure of the plasma velocity u, which happen very close to the shock, where the diffusive shock acceleration of particles is provided by the gradient of this velocity, du/dx.

Therefore, in the present paper, we address the CR acceleration at the post-adiabatic shocks by expanding the approach developed by Petruk & Kuzyo (2024). The main idea in that study was that the nonuniformity of the downstream flow on the length scales involved in the particle acceleration could affect the CR spectrum. This is analogous to the nonlinear acceleration, where the flow velocity is not spatially uniform upstream due to the back-reaction of CRs. In our case, the difference is in the physical reason for this flow nonuniformity: it is a purely hydrodynamical effect taking place downstream of the shock.

The radiative losses of the plasma are proportional to n²Λ(T, τ) where n is the density and Λ(T, τ) the cooling coefficient, which depends on the temperature T and the ionization timescale τ. Generally, radiative losses become prominent when the shock speed V slows down to values that provide a temperature decrease T ∝ V², which pushes Λ(T) to its maximum. However, the post-adiabatic evolution is not a feature of only old SNRs. It is also natural for the SNRs interacting with dense ambient material, such as molecular or atomic clouds, since radiative losses are proportional to n². Interacting SNRs are expected to be sources of hadronic gamma-rays. It is worth studying, therefore, what kind of changes in the spectrum of particles accelerated at the post-adiabatic shock one could expect.

For the present paper, we adopted 1D MHD numerical simulations. As become clear at the end of the study, CR acceleration at post-adiabatic shocks merits further investigation in higher dimensions. Therefore, the present paper should be considered a first-cut assessment of the problem.

2 Motivations

2.1 MHD features of the post-adiabatic SNRs

In our previous papers (Petruk et al. 2016, 2018), we studied the behavior of the magnetic field in SNRs during the transition phase from adiabatic to fully radiative evolution as well as its role in the modification of the dynamics of post-adiabatic flows. In particular, we explored the evolution of radial profiles of density, pressure, and temperature, as well as of the radial and tangential components of MF. In the present paper, we focus instead on the flow velocity profiles.

We have performed numerical simulations of 1D MHD models in spherical coordinates which reflect the evolution of a sector in SNR, with a small solid angle. The initial conditions correspond to the Sedov–Taylor stage. They are provided by setting up a strong point explosion with energy E₀ = 10⁵¹ erg, released isotropically, in the uniform medium with hydrogen number density n₁ = 1 cm⁻³ and magnetic field with strength B₁. We consider two models: model 1 corresponds to a parallel shock (i.e., ambient MF B₁ parallel to the shock velocity V) with B₁ = 5 µG, while model 2 assumes a perpendicular shock with B₁ = 10 µG. The shock modification due to the back-reaction of CRs was not considered.

The simulations were performed using the PLUTO MHD code (Mignone et al. 2007) on a static uniform grid divided into 100 000 computational cells. The grid size was set to 32 pc. In the numerical setup, we used the Characteristic Tracing time integrator and the HLL Riemann solver. The radiative losses of plasma were treated under the assumption of ionization equilibrium. The setup is similar to that adopted by Petruk et al. (2016, 2018). We refer to these papers for a general overview of the post-adiabatic shock properties in the presence of an ambient magnetic field with different orientations to the shock velocity.

Fig. 1 shows, for reference, the time dependence of the shock radius R, the shock speed V, and the expansion parameter m. This parameter is useful in determining the evolutionary stage of a supernova remnant. Numerically, m = 0.4 during the Sedov stage, while in the course of the fully radiative phase, its time evolution (black dashed line in Fig. 1) is given by an analytical solution (Bandiera & Petruk 2004, and references therein).

In this figure, the two reference times mark approximate limits of the post-adiabatic stage (for a review of different reference times related to cooling of SNRs, see Petruk 2005). The first one is the transition time, when the flow starts to deviate from the Sedov (1959) solution due to the radiative losses of plasma, which become prominent around this time. An analytical estimate for it is given by the expression (Cox 1972; Blondin et al. 1998) $t_{\text{tr}}=2.83\times {10}^4{E}_{51}^{4/17}{n}_1^{-9/17}$ yr where E₅₁ is E₀ in 10⁵¹ erg. The thin cold dense shell is formed around the shell formation time (Cox & Anderson 1982; Cioffi et al. 1988), which is t_sf = 1.83 t_tr for our models. From this time on, the shock is fully radiative, i.e., all thermal energy of the newly shocked ISM material is quickly radiated away. Between these two milestones, the shock is partially radiative, with radiative losses becoming increasingly more important in the SNR’s dynamics. Radiative losses considerably affect the structure of the flow, particularly the flow speed.

Animations demonstrating changes in the flow structure in this stage are presented in the Appendix A.1. The downstream profiles of the flow speed in the shock reference frame u₂, magnetic field strength B₂, and plasma density n₂ in the two models for a few moments of time are shown in Fig. 2. The spatial distributions in this figure are shown versus lg(x) to demonstrate that changes during the post-adiabatic era also affect very small distances from the forward shock. Our simulations resolve the structures down to 10⁻⁵ of the shock radius for t ≥ t_tr. The flow speed in the reference frame of the shock is u₂(x) = V − v₂(R − x) where v₂(r) is the flow speed in the laboratory frame, r and x = R − r are the distances from the center and from the shock, respectively.

We see, in agreement with the results presented by Petruk et al. (2016, 2018), that the MF rapidly decreases downstream of the parallel shock, while it increases behind the perpendicular shock (middle panels in Fig. 2). Actually, the spatial structure of plasma is not affected by the radial MF because plasma flows along the MF lines without any resistance; this model behaves in the same way as a nonmagnetic one. Instead, the tangential MF lines limit the plasma compressibility, which increases because of radiative losses. As a result, the plasma density in model 2 is lower than in model 1 (lower panels in Fig. 2).

The most important HD parameter for the present paper is the flow velocity. It varies considerably in space and time when radiative losses come into play (the upper panels on Fig. 2 and animations in the Appendix A.1). Contrary to the previous evolutionary stages, changes in the flow velocity happen on the length scales less than a few percent of the shock radius, i.e., on the scales involved in the cosmic ray acceleration.

The post-adiabatic flow has a downstream region, between the coordinates x_A and x_B, where the flow velocity is inverted, u₂ < 0. This happens because the pressure drops right after the shock due to radiative losses of plasma there, while the deeper interior remains hot. In the presence of the tangential component of MF, the magnetic pressure increases in this region due to plasma compression, and the magnetic pressure diminishes the overall effect, that is, the flow speed (Fig. 2) and the size of the region with negative u₂ (Fig. 3).

The region with the inverted flow velocity appears first around time 1.49t_tr and extends almost to the end of the post-adiabatic phase, up to 1.75t_tr = 0.97t_sf. Thus, the highest boost in the CR acceleration efficiency should happen around this period.

In summary, there are several prominent features in the shock and flow velocity that can affect the particle acceleration at the post-adiabatic shocks.

• The shock slows down (u₁ diminishes) faster than during the previous adiabatic phase.

• The ratio u₁/u₂(x) increases in a certain range of x as radiative losses start to influence the shock.

• This ratio may become quite large as u₂ approaches zero.

• During some time interval, there is a region where u₂ becomes negative, and particles are not advected away from the forward shock, but instead are pushed back toward it.

• The spatial extent of this region could reach about 10% of the SNR’s radius.

• Changes in the spatial structure of u₂ occur even at distances ~10⁻³R, that is, so close to the shock that they even influence the acceleration of the radio-emitting electrons (as we show below).

• The MF component, orthogonal to the flow velocity, can affect the velocity profile.

Only the first property leads to a reduction in the acceleration efficiency; all others increase it. The shape of the CR spectrum and the maximum energy they can attain are determined by a combination of these trends.

Fig. 1 Top: Evolution of the radius R (dashed lines) and speed V (solid lines) of the forward shock. Bottom. Time dependence of the expansion parameter m. The three vertical lines mark the times ttr (left), 1.49ttr (middle), and tsf (right).

Fig. 2 Spatial distributions of the flow speed u2, magnetic field strength B2, and number density n2 downstream of the shock at several moments of time. The left column is for model 1 and the right column for model 2. The coordinate x marks the distance from the shock, which is used for considerations of the CR acceleration. We note that the horizontal axes are in the log scale. This is necessary to describe equally well the variations on both small and large scales relative to the shock position, which are probed by particles with momenta spanning several decades. The shock is located on the left of each plot, at x/R = 0. The SNR center is on the right, at the point with lg(x/R) = 0.

Fig. 3 Distances xA and xB from the shock for the two models. The flow speed u2 < 0 between the two solid lines. The vertical dotted lines mark the times 1.49ttr and tsf. The flow speed u2 > 0 in the whole domain before 1.49ttr.

2.2 Observational evidence

Simulations show that the flow velocity profile changes continuously after the end of the adiabatic era, as well as the conditions for particle acceleration around the forward shock, which is progressively affected by the increasing radiative losses. Several questions arise, including whether there is observational evidence supporting this.

We also need to understand the kind of dependence that we should look. The ratio of the pre- to the post-shock velocities σ = u₁/u₂ increases during some time after t_tr. This occurs in regions very close to the shock and could affect the acceleration of the radio-emitting electrons. If so, then the radio index α = (s − 1)/2 should reflect changes in the electron spectral index s = (σ + 2)/(σ − 1). That is, the absolute value of α should decrease as the SNR progresses into the post-adiabatic phase. With this trend in mind, we may consider the following two observational hints.

First, we consider a sample of SNRs that have radio and GeV γ-ray emission detected, as well as estimations for their age. Zeng et al. (2019) have analyzed the multi-wavelength spectra of 35 such SNRs under a uniform approach. The authors provide estimates for a spectral index and SNR age. These are SNRs likely interacting with dense clouds, so they are rather in the post-adiabatic stage, with the age ≳t_tr. A broadband (from radio through X-rays to γ-rays) spectral fit provides the value of the particle spectral index. The electron spectral indices of these SNRs estimated in this reference are shown in Fig. 4 as a function of the SNR age. A systematic decline in the index with age is apparent. It means that, independently of local conditions, the ratio u₁/u₂ probed by the radio-emitting electrons increases over time for older SNRs, which could suffer (to a different extent) from radiative losses. Thus, this trend agrees with our hypothesis.

Another possibility is to look for a single SNR where a fraction of a shock interacts with a dense environment. If different parts of the shock enter regions with different densities, then we may expect a lower value of α at locations where the density is higher. There, the shock is affected by radiative losses to a larger extent because the losses are proportional to n². Impressively, the correlation between the radio spectral index and the density is apparent in the supernova remnant Kes 73 (G27.4+0.0). The left plot in Fig. 5 shows maps of the radio spectral index of Kes 73 obtained from the radio observations on the Green Bank Telescope and the Very Large Array (Ingallinera et al. 2014). The CO intensity map in the region around Kes 73 is shown in the right plot. We used CO High-Resolution Survey data (Dempsey et al. 2013)¹ to produce this map. It is important to note that this SNR is young and, as a whole, is not in the partially radiative phase, let alone in the fully radiative phase (Borkowski & Reynolds 2017). We believe that only a small fraction of the forward shock is encountering a molecular cloud, and radiative effects are initiating their development in these isolated regions. Some aspects of the hypothesis that the variation of the radio spectral index in Kes 73 could be due to the shock interaction with dense cloud are discussed in Appendix B.

Fig. 4 Electron spectral index vs. age for a sample of SNRs with the radio and GeV γ-ray emission detected, from the data given by Zeng et al. (2019). The red dashed line represents a linear fit between the index and logarithm of age, with the slope k = −0.37 ± 0.05. The red area shows the 1σ uncertainty. We excluded MSH 15-56 from the sample because its radio emission is likely dominated by a pulsar. There are a few SNRs with two values of α estimated in this reference. We calculated the average α from the two values if the uncertainty is due to various estimates of the ambient density (Tycho, G166.0+4.3, S147, CTB 73B). We use only one value for W51C, RX J0852-4622, namely that which corresponds to the same spectral model as used for other SNRs.

3 CR acceleration at a post-adiabatic shock

3.1 General considerations

We consider two reference frames. The first is the observer’s frame, which is at rest relative to the explosion center. In this frame, the shock velocity V in the 1D sector of the SNR is positive. The second is the shock frame, often used in the description of CR acceleration. In this frame, the shock is at rest, and the positive direction of the coordinate x is defined along the flow velocity of the upstream medium.

The question arises of what could happen with cosmic ray acceleration at the shock when, for some time, material in some region downstream has, in the frame of the observer, a radial velocity higher than the shock velocity itself. Under such conditions, this translates into negative velocities of this material in the reference frame comoving with the shock.

Such a negative velocity in the shock reference frame means that the plasma moves toward the shock in the downstream region, not out of it (as commonly assumed). Such a situation may naturally occur in SNRs right after the end of the Sedov phase when the radiative losses of plasma in the vicinity of the shock become prominent.

As a preliminary point, we have to keep in mind that in the standard model for diffusive particle acceleration, in the shock reference frame, the flow velocity is always taken to be positive, namely the medium moves toward the shock surface in the upstream flow, while it moves away from it downstream. This has a clear consequence in the evolution of particles, as depicted, for instance, in Bell’s model for particle acceleration in shocks. Particles can be accelerated as long as they are hopping back and forth across the shock. The ingredients of the model are the following (Jones & Ellison 1991):

1. the average increase in momentum experienced by a particle over a single cycle; it is ∆p/p ≃ 4/3 (u₁ − u₂)/c, where (u₁ − u₂) is the velocity difference between the upstream (labeled with 1) and the downstream (labeled with 2) flow;

2. the probability that a particle in the downstream region can reach the shock for another cycle, traveling against the general downstream flow.

In this model, quantities that can change and that depend both on the ambient magnetic field and the particle momentum are: (i) the maximum distance from the shock a particle can reach on the shock upstream x_p ∼ D/u, where D is the diffusion coefficient; the same expression is also valid for the maximum distance it can reach downstream, under the condition that it still has some chance to return to the shock again; say D = η r_gc/3, where r_g = r_Lγ = p c/(e B) is the radius of gyration (r_L being the Larmor radius) and η ≥ 1; (ii) the residence time on each side is nearly the previous quantity divided by c, so the time required for the particle to complete one full cycle is t_cycle = (4/c)(D₁/u₁ + D₂/u₂).

This framework requires that the particles eventually flow away in the downstream flow; otherwise a runaway situation develops in which no particle can escape, but they all are accelerated forever. The standard model of diffusive particles acceleration assumes steady-state conditions; while a situation with a negative u₂ and a shock strong enough to efficiently accelerate particles cannot continue forever, otherwise we end up with an infinite density in the vicinity of the shock and infinite energy in the accelerated particles.

We consider a strong nonrelativistic shock discontinuity. In units of the shock speed, the immediate downstream velocity after the shock is u₂(immediate) = u₁/σ with the compression factor σ = 4. We then assume a downstream velocity profile in which, at a given time, the downstream velocities turn negative from the distance x_A to the distance x_B from the shock. We recall that, in the shock frame, the shock is at x = 0, x_A and x_B are positive and x_A < x_B (Fig. 6). The diffusion length x_p2 could be less than x_A, in between these two coordinates, or larger than x_B, depending on momentum.

For the particles, as soon as their diffusion length x_p < x_A, essentially nothing changes from the common picture of acceleration because during their acceleration, they do not notice a velocity inversion further on. In particular, only a fraction of particles can return to the shock against the flow, specifically those with an x-component of their velocity w at a given point x satisfying w_x < −u₂(x). Instead, particles in the region of the negative u₂ will be easily directed toward the shock by the flow. Only particles with w_x > −u₂(x) can escape this region toward the far downstream.

The idea is that, when the accelerated particles gain the momentum p_A corresponding to the diffusion length x_p = x_A, they would remain trapped between the shock and the region with negative u, until they reach the value p_B relevant to x_p = x_B that allows them to be advected downstream or otherwise, if the hydrodynamic evolution proceeded faster than the acceleration process, until the region of negative u eventually disappears. Since generally x_B is a nonnegligible fraction of the SNR radius (Fig. 3), this means that particles, if trapped long enough, could be accelerated up to very high energies.

3.2 Diffusion length

The question is how to compute more precisely x_p, considering the spatial variation of hydrodynamic quantities, such as u and B.

In the case with constant profiles of D(x) and u(x) > 0, the diffusion length x_p of particles with a given p is such that ux_p/D = 1. One could impose a generalization to nonconstant profiles by using the simplest dimensional criterion as $${\displaystyle {\int }_0^{x_{\text{P}}}\frac{\left|u(x^{\prime })\right|}{D(x^{\prime })}d{x}^{\prime }=1}$$(1)

where the absolute value of velocity is taken because x_p, dx′, and D(x′) are always positive while u(x′) can change the sign.

We consider the diffusion coefficient D = ηcr₀ (r_g/r₀)^1/κ/3, where r₀ is a reference length, and κ = 1, 2, 3 for the Bohm, Kraichnan, and Kolmogorov diffusion, respectively. By substituting equation (1) for x_p with this diffusion coefficient, we arrive at the function $\overline{p}(x_{\text{p}})$, which relates the momentum p to a given diffusion length x_p: $$\overline{p}(x_{\text{p}})={\left(\frac{3}{\eta c{r}_0}\right)}^{\kappa }\frac{e{r}_0}{c}[{\displaystyle {\int }_0^{x_{\text{P}}}\left|u(x)\right||B(x){|}^{1/\kappa }dx{]}^{\kappa }}.$$(2)

It is important that the derivative is always $d\overline{p}/d{x}_{\text{p}}\ge 0$. Therefore, there is no multiplicity in x_p for any momentum at any time.

Figure 7 shows the momentum p as a function of the downstream diffusion distance x_p in our MHD models. It is calculated as a solution of the equation $p=\overline{p}(x_{\text{P}})$ with the function $\overline{p}(x)$ given by Eq. (2). Particles with (x, p) below the lines cannot return to the shock for acceleration.

For the sake of comparison, we also plot the dashed lines that correspond to the common approach, with a uniform flow downstream. We see from this figure (cf. Fig. 2) that the nonuniformity of u₂(x) and B₂(x) may affect particles with quite small diffusion length, down to ∼10⁻³ of the shock radius, for the Bohm diffusion. These are the particles with momenta as low as ∼10 GeV in model 1 (yellow line in Fig. 7), i.e., those emitting radio waves. If the Bohm factor η is considerably larger than unity or the Kolmogorov diffusion operates (as observed, Shimoda et al. 2018; Petruk & Kuzyo 2025) then the diffusion coefficient and therefore x_p are larger for the same momenta. Then, the radio-emitting electrons could probe larger distances from the shock, at least some time during the post-adiabatic phase.

The flow velocity profile changes over time with a stronger spatial variation around 1.7t_tr (yellow lines in Fig. 7) when it affects particles with the smallest momenta. At different times and different distances, particles probe different compression factors u₁/u₂(x_p) and get different ∆p that should be reflected in their momentum distribution and in the radio spectrum they emit.

Fig. 5 Distribution of the radio spectral index over Kes 73 measured between 1.4 and 5.0 GHz (left) and the 12CO (J = 3 → 2) intensity map (Dempsey et al. 2013) integrated over the velocity range between 94 and 95 km/s, in units of K km/s (right). The contours correspond to the values of the radio index 0.8, 0.6, and 0.4. Liu et al. (2017) also present 12CO (J = 1 → 0) intensity maps around Kes 73 at various velocities. Those in the range 86–97 km/s correlate with the radio index map as well.

Fig. 6 Spatial structure of the flow velocity (blue line) in the reference frame of the shock, which is located at x = 0. The diffusion length of a particle with momentum p is xp1 upstream and xp2 downstream. These are the largest distances the particle can diffuse and return to the shock to be further accelerated. The shock compression factor σ = u1/u2 is the ratio of the immediate pre- to post-shock velocities. Instead, σp(p) = up1(p)/up2(p) is an effective compression factor “seen” by a relativistic particle with momentum p; it varies with momentum. The shock here is not CR-modified, and therefore up1(p) = u1 upstream, for any momentum. For the downstream region, one may use an approximation up2(p) ≈ u2(xp2) (see Sect. 3.4).

Fig. 7 Minimum momentum p of particles that may reach the distance xp downstream and return to the shock to continue acceleration. It is calculated from the equation $p=\overline{p}(x_{\text{p}})$ for Bohm diffusion. The colored lines correspond to a few time moments (see inset). The dashed lines show p calculated in a common approach, from xp = D2(p)/u2 with uniform u2 and uniform B2 = 5 µG. The top plot is for model 1 and the bottom plot is for model 2.

3.3 Increment in momentum

Another aspect of particle acceleration at the post-adiabatic shocks is related to the momentum gain per cycle. We consider the interaction of a particle (mass m, speed w ∼ c) with scattering centers (mass M ≫ m with velocities u₁ > 0 upstream and u₂ > 0 downstream, spatially constant in their domains and with absolute values u₁ ≪ w, u₂ ≪ w). The initial velocity of a particle is w_o = −w < 0 toward a scattering center in the upstream flow. After interaction there, its velocity is w′ = w + µu₁ (to the order u₁/w) where 0 ≤ µ ≤ 1 is the cosine of the angle between the particle and the flow velocities. Then the particle crosses the shock, interacts with a scattering center in the downstream flow, and its velocity becomes w″ = −(w′ − µu₂). Therefore, the increment in momentum for a given particle in one acceleration cycle is ∆p/p = (w″ − w_o)/w_o = µ(u₁ − u₂)/w. After averaging over the particle flux proportional to wµ, we arrive at a known formula ∆p/p ≃ 4/3 (u₁ − u₂)/c. Instead, if u₂ < 0, then w″ = −(w′ + µ|u₂|) and ∆p/p = µ(u₁ + |u₂|)/w; and the flux average yields ∆p/p ≃ 4/3 (u₁ + |u₂|)/c.

During one cycle, particles with momentum p experience many collisions, but the whole process proceeds in a way that the overall momentum increase during a cycle is determined primarily by the outermost interactions, at the locations x_p1 and x_p2 where the flow velocities are u_p1 and u_p2 (see Fig. 6). It seems that the velocity gains and losses in intermediate interactions between these two points, on average, nearly cancel themselves out. Such a property is demonstrated by a solution for the CRs-modified shock where the spectral index at the momentum p may be estimated from the effective compression u₁(x_p1)/u₂. Therefore, by introducing σ_p = u_p1/|u_p2|, the general formula becomes $$\frac{\text{Δ}p}{p}=\frac{4}{3}\frac{u_{\text{p1}}}{v}\frac{{\sigma }_{\text{p}}-\text{sgn}(u_{\text{p2}})}{{\sigma }_{\text{p}}},$$(3)

where sgn(x) is the sign function. We can see the two effects from this expression: the first is due to the negative u₂ and the second is due to the change in the effective compression σ_p, which may be quite high if u_p2 approaches zero. The factor (σ_p − sgn(u_p2))/σ_p increases linearly with the decrease in the ratio u_p2/u_p1 (e.g., from 0.75 for u_p2/u_p1 = 0.25 to 1.25 for u_p2/u_p1 = −0.25). This means that if there is a downstream region of the negative flow velocity, the increase in momentum per acceleration cycle could be about two times the classic value (for uniform and positive u₂(x)).

In contrast, the shock decelerates considerably over time during the post-adiabatic phase. Therefore, the temporal variation of ∆p/p is determined by changes in the HD structure of the flow as well as in the shock speed.

3.4 Spectrum of accelerated particles

We now consider the particle distribution function. The solution of the steady-state kinetic equation for CRs accelerating in the spatially variable flow downstream has been found by Petruk & Kuzyo (2024). The distribution function is $$f_{\text{o}}(p)=\frac{\xi n_1}{4\pi p_{\text{o}}^3}\frac{3u_1}{u_{\text{p1}}-u_{\text{p2}}}\exp \left[-{\displaystyle {\int }_{p_{\text{o}}}^p\frac{3u_{p^{\prime }1}}{u_{p^{\prime }1}-u_{p^{\prime }2}}\frac{d{p}^{\prime }}{p^{\prime }}}\right],$$(4)

where ξ is the injection efficiency. The spatial variation of the flow velocity in the downstream region is accounted for through u_p2 and through u_p1 in upstream locations: $$u_{\text{p1}}(p)=u_1-{\displaystyle {\int }_{-\infty }^{-0}\frac{f(x,p)}{f_{\text{o}}(p)}\frac{du}{dx}dx},$$(5) $$u_{\text{p2}}(p)=u_2+{\displaystyle {\int }_{+0}^{+\infty }\frac{f(x,p)}{f_{\text{o}}(p)}\frac{du}{dx}dx}.$$(6)

We consider a (reasonable) approximation that the distribution function for particles undergoing acceleration is f (p, x) = f_o(p)ℋ(x_p − x), where f_o(p) is the distribution at the shock and ℋ(x) is the Heaviside step-function, which means that particles with the momentum p which continue acceleration are within their diffusion length x_p from the shock. Then u_p1 = u₁(x_p1) and u_p2 = u₂(x_p2).

The local (in the momentum space) index of the particle distribution function is given by $$s\equiv -\frac{d\ln f_{\text{o}}}{d\ln p}=\frac{3u_{\text{p1}}}{u_{\text{p1}}-u_{\text{p2}}}+\frac{d\ln (u_{\text{p1}}-u_{\text{p2}})}{d\ln p}.$$(7)

4 Toy model A: Case of a constant gradient of the flow velocity

We examine properties of the particle distribution in a toy model that results in a simple analytic solution. Namely, a shock propagating with a constant speed and a downstream flow with a constant gradient du/dx. The CR back-reaction is negligible, and therefore the flow velocity upstream u₁ is spatially constant and u_p1 = u₁. The velocity u_p2 = u(x_p). The flow velocity at a point x downstream is u(x) = u₂ + [du/dx]x; u₁ > 0 and u₂ = u₁/4 > 0 are immediate pre- and post-shock values. By substituting Eq. (1) with this u(x) and considering a spatially constant diffusion coefficient, we derive a quadratic equation for x_p with a real positive solution $$\frac{x_{\text{p}}}{R}=\frac{1}{\mu }\left[\sqrt{1+2\mu \frac{{\tilde{x}}_{\text{p}}}{R}}-1\right],\mu \equiv \frac{\left[du/dx\right]}{u_2/R}>-\frac{R}{2{\tilde{x}}_{\text{p}}}$$(8)

where ${\tilde{x}}_{\text{p}}=D/u_2$. As an example, the parameter µ ≈ 1 for Sedov shock (Taylor 1950). The ratio ${\tilde{x}}_{\text{p}}/R\ll 1$ and we may use the first-order approximation $x_{\text{p}}\simeq {\tilde{x}}_{\text{p}}\propto p^{-1/\kappa }$.

The particle distribution at the shock for this problem can be obtained in analytic form from the general solution (4): $$f_{\text{o}}(p)=\frac{\xi n_1}{4\pi p_0^3}\frac{3u_1}{u_1-u_2}{\left[1-\frac{[du/dx]x_{\text{p}}}{u_1-u_2}\right]}^{\delta }{\left(\frac{p}{p_0}\right)}^{-3u_1/(u_1-u_2)},$$(9)

where δ = [(3κ − 1) u₁ + u₂]/(u₁ − u₂). This is essentially the classic spectrum f_o ∝ p⁻⁴ for momenta p, which corresponds to the diffusion lengths $$x_{\text{p}}\ll \frac{(u_1-u_2)}{\left|du/dx\right|}\iff \frac{x_{\text{p}}}{R}\ll \frac{3}{\left|\mu \right|}.$$(10)

At momenta relevant to x_p/R ≳ 0.3/|µ|, the shape of the particle spectrum depends on the sign of the velocity gradient. When du/dx > 0, the local spectral index s increases with momentum p, leading to a kind of high-energy cutoff in the spectrum. Conversely, if du/dx < 0, the spectral index s decreases with p, resulting in a harder spectrum. For example, in the case of the Sedov-Taylor shock, where du/dx > 0, the spectrum steepens at very high p (figure 3 in Petruk & Kuzyo 2024).

5 Toy model B: Maximum momentum

5.1 Numerical solutions of the time-dependent kinetic equation

In the present section, we consider another toy model to illustrate and further explore the effect of the radial profile of the flow speed downstream on the particle acceleration.

Specifically, we use RATPaC (Radiation Acceleration Transport Parallel Code) (Telezhinsky et al. 2012, 2013; Brose et al. 2020) to derive a numerical solution of the time-dependent kinetic equation for particle acceleration with a given flow velocity profile. Namely, we consider an unmodified shock with the flow speed in the laboratory frame as $$v(r)=\{\begin{array}{ll}{v}_{\text{s}},\hfill & \text{for\hspace{0.17em}}0\le r\le r_{\text{B}},r_{\text{A}}\le r\le R\hfill \\ v_{\text{AB}},\hfill & \text{for\hspace{0.17em}}{r}_{\text{B}}<r<r_{\text{A}}\hfill \\ 0,\hfill & \text{for\hspace{0.17em}}r>R\hfill \end{array},$$(11)

where v_s = 3V/4 is the common post-shock value of the flow speed. We consider the following parameters v_AB/V = 3, r_B/R = 0.97, r_A/R = 0.99. This converts to the shock frame as u = V − v, x = R − r, i.e., $$u(x)=\{\begin{array}{ll}{u}_1,\hfill & \text{for\hspace{0.17em}}x<0\hfill \\ u_{\text{AB}},\hfill & \text{for\hspace{0.17em}}{x}_{\text{A}}<x<x_{\text{B}}\hfill \\ u_2,\hfill & \text{for\hspace{0.17em}}0\le x\le x_{\text{A}},x_{\text{B}}\le x\le R\hfill \end{array},$$(12)

where u₁ = V, u₂ = u₁/4 with the values u_AB = −2u₁, x_A/R = 0.01, x_B/R = 0.03.

The solution for such a problem is compared to the solution for the uniform velocity profile v(r) = v_s for 0 ≤ r ≤ R. The shock speed V = 500 km/s and the shock compression factor σ = 4 are constant over time, while the radius is set to vary as R = Vt. The particles are continuously injected with a constant injection efficiency ξ at a few MeV. The diffusion coefficient is Kolmogorov-like D = 2.5 · 10⁹(w/3)(pc/qB)^1/3 cm²/s where w is the particle speed. The strength of the magnetic field is B₁ = 5 µG and B₂ = σ_B B₁ with ${\sigma }_{\text{B}}=\sqrt{11}$. The diffusion coefficient changes accordingly when a particle crosses the shock. Simulations start at t₀ = 20 yr.

Fig. 8 shows distributions of particles as a function of momentum. We emphasize that the model considered exaggerates the effect of the nonuniformity of the downstream flow. In particular, by keeping x_A/R and x_B/R constant while the shock expands, we increase the region of negative velocity over time. The particle spectra derived for the toy model are nonphysical, with unrealistic energies transferred to CRs. A particular reason is that we set V(t) = const in the simulations. This is done to isolate the effect of the flow velocity profile on the particle acceleration. In reality, the speed of the post-adiabatic shock decreases efficiently with time.

There are some features in the particle spectra in model (11) noticeable in Fig. 8.

• Compared to the case of the uniform flow (represented by dashed lines), one can see that the significantly larger energy is transferred by the shock into the high-energy CRs (solid lines), and this energy grows with time.

• The maximum momenta of particles are larger (solid lines) than in the classic case (dashed lines).

• After the first spectral break in the spectrum (around 0.1m_pc), the spectra are nearly power laws up to some momentum and then, at a higher momentum, the cutoff happens around p_max.

• We have also calculated a model with a wider region of negative u₂, with x_B = 0.05R. The maximum momentum is higher if the region with negative u₂ is wider.

We compare these numerical results with analytical estimates by considering the case with x_B = 0.03R. Particles scattering back to the shock from distances x ≤ x_p are not affected by the negative u₂. Their diffusion length is x_p = D₂(p)/u₂. From this relation, we may estimate the momentum p_A corresponding to the diffusion length x_A = 0.01R. It is p_A = 0.19m_pc. Indeed, the spectra in Fig. 8 (solid lines) deviate from a simple power law with s = 4 around these momenta.

The situation is different for particles with p ≳ p_A because their x_p > x_A. The diffusion length calculated with Eq. (1) for t = 420 yr is shown in Fig. 9 (top plot). The momentum p_B corresponding to x_B = 0.03R is p_B = 915m_pc, as from the analytical formula (1). From Fig. 8, we see that the numerical spectrum for t = 420 yr indeed deviates from the hard power law (i.e., with a spectral index s < 4) around this momentum.

Particles with momenta p_A < p < p_B are scattered back to the shock from smaller distances x_p from the shock (orange line) compared to the case of the uniform u₂ (blue line). For example, particles with momentum p_B return to the shock from a distance 0.17R in the model with the uniform flow but from 0.03R in the model with the flow as in Eq. (12). Clearly, the time needed for one acceleration cycle is shorter in the latter case, and particles can be accelerated to higher p_max during the same period of time. We consider this further in Sect. 5.2.

Under the approximation f (x, p) = f_o(p)ℋ(x_p − x), the integral (4) for the analytical solution may be separated into three integrals, in the domains (p_o, p_A), (p_A, p_B), (p_B, p_max). The flow velocity is constant in each of these domains. Therefore, the distribution function f_o(p) should be given by a sequence of power laws with spectral indices s = 4 until p_A, s = 1 until p_B, and s = 4 until p_max. The numerical spectra in Fig. 8 follow this sequence with the addition of the high-energy cutoff, which smooths the last power-law part of the analytical spectrum. Some additional details on the shape of the CR spectrum in a numerical solution are given in Appendix C.

In summary, if there is a region x > x_A with a ratio u₂/u₁ different from the value relevant for compression at a strong adiabatic shock 0.25, then the spectrum of accelerated particles with p ≳ p_A has a spectral index s = 3/(1 − u₂/u₁) that is different from the classic s = 4. In particular, if u₂/u₁ > 0.3, the index of the distribution function N(p) ≡ p² f_o(p) ∝ p^−Γ is Γ ≳ 2.3 as needed to fit gamma-ray spectra in some middle-aged SNRs interacting with dense molecular clouds, for example IC443 or W44 (e.g., Ackermann et al. 2013). The situation u₂/u₁ > 0.3 occurs if the shock suddenly interacts with a dense material. Then, the shock speed rapidly decreases while the internal material still has a speed corresponding to the pre-interaction times.

Fig. 8 Volume-integrated spectra for particles inside SNRs at different times (marked in the plot legend, in years). The solid lines represent the model with parameters uAB = −2u1, xA/R = 0.01, xB/R = 0.03. The spectra shown by the dashed lines are calculated for the uniform flow speed.

Fig. 9 Diffusion length xp2(p) (top) and mean acceleration time ⟨tacc(p)⟩ (bottom) calculated in the analytical approach (i.e., it is time-independent). The model with xA = 0.01R, xB = 0.03R, uAB = −2u1 is shown with the orange line; it is calculated with Eqs. (1) and (15). The model with uniform u2 is shown with the blue line, as from xp = D2/u2 (then simply xp ∝ p1/3) and (14) for acceleration time. The green dashed line shows xp1(p) upstream, and the gray dashed line corresponds to t = 420 yr.

5.2 Mean acceleration time

Here we discuss how we can estimate the acceleration time and the maximum momentum if the flow velocity is not spatially constant. It is well known that the average acceleration time for the case of spatially constant diffusion coefficient D and flow velocities u₁, u₂ is (Drury 1983) $$\langle t_{\text{acc}}(p)\rangle =\frac{3}{u_1-u_2}{\displaystyle {\int }_{p_{\text{o}}}^p\left(\frac{D_1(p^{\prime })}{u_1}+\frac{D_2(p^{\prime })}{u_2}\right)\frac{d{p}^{\prime }}{p^{\prime }}.}$$(13)

With D₁, D₂, u₁, u₂ spatially constant in their domains and the diffusion coefficient D = ηcr₀ (r_g/r₀)^1/κ/3 used in Sect 5.1, Eq. (13) transforms to $$\langle t_{\text{acc}}(p)\rangle =\frac{3{\beta }_{\text{a}}{D}_1(p)}{V^2},{\beta }_{\text{a}}=\frac{\kappa \sigma \left(1+\sigma /{\sigma }_{\text{B}}^{1/\kappa }\right)}{\sigma -1}$$(14)

for p ≫ p_o, where β_a accounts for the fractions of time particles spend upstream and downstream. It is β_a = 14.7 for σ = 4, ${\sigma }_{\text{B}}=\sqrt{11}$ and Kolmogorov diffusion. This dependence is shown in Fig. 9 with the blue line. The maximum momentum which may be reached at t = 420_yr is p_max = 55m_pc. This corresponds to the maximum momentum in the spectrum from the bottom right plot in Fig. 8.

However, the maximum momentum in the spectrum for t = 420 yr is considerably higher in the region where the negative u₂ is accounted for (Fig. 8 bottom left plot). Therefore, the formula (13) should be modified.

The two ratios under the integral in Eq. (13) are actually estimates for x_p(p). Therefore, we suggest an ansatz for the acceleration time that incorporates the spatial variations of the flow speed and magnetic field $$\langle t_{\text{acc}}(p)\rangle \simeq \frac{3}{u_1(x_{\text{p1}})-u_2(x_{\text{p2}})}{\displaystyle {\int }_{p_{\text{o}}}^p\left(x_{\text{p1}}(p^{\prime })+x_{\text{p2}}(p^{\prime })\right)\frac{d{p}^{\prime }}{p^{\prime }}},$$(15)

with x_p(p) given by the equation (1). The orange line on the bottom plot in Fig. 9 shows the acceleration time calculated with the expression (15). The acceleration time is shorter for particles with momenta p > p_A if there is a region with an inverted flow velocity downstream.

Since acceleration is faster, the maximum momentum is limited by the available time ⟨t_acc(p_max)⟩ = t may be as high as 915m_Hc in 420 yr at the shock with v_AB = 3V, while it could reach just 47m_Hc at the shock with uniform u₂ (Fig. 9).

To demonstrate the reliability of the expression (15), we compare the maximum momenta it gives with those derived from the numerical solution in Fig. 10 (orange line and dots). Impressively, we can see an almost perfect correspondence. In parallel, the blue line is calculated with Eq. (14) and it corresponds to the numerically derived p_max in the model with uniform u₂. The numerical solution clearly demonstrates the effect shown by analytical formulae, namely, the higher p_max if there is a region with a negative flow velocity downstream.

Finally, the orange line in Fig. 10 perfectly matches p_max ∝ t ³, as expected. The dependence p_max ∝ t ^κ (with κ = 3 for the Kolmogorov diffusion coefficient) follows from a simple relation t_acc ∝ D(p) ∝ p^1/κ.

Fig. 10 Maximum momenta pmax for two numerical solutions (dots) and from analytical expressions (solid lines). The blue line is from Eq. (14) and the dots are from the numerical solution with uniform u2 (dashed lines in Fig. 8). The orange line is calculated with Eq. (15) with xp1 and xp2 from (1) and the dots are from the numerical solution with xB = 0.03R, uAB = −2u1 (solid lines on Fig. 8). The value of pmax is measured at the level where the distribution function drops below the power-law extrapolation in e times for the orange dots and in 10 times for the blue dots.

6 Evolution of spectra in the MHD models

In the present section, we return to our MHD models, namely model 1 and model 2.

Fig. 11 shows the mean acceleration time ⟨t_acc(p)⟩ calculated for our model 1 in the two approaches. The solid lines correspond to the detailed calculations with Eqs. (15) and (1) which account for the spatial profiles of MHD parameters. The dashed lines represent calculations with Eq. (14) assuming the spatially uniform flow. A given line corresponds to a system configuration (e.g., the shock speed, profile of the flow velocity) at the time indicated in the legend. The line shows the time a particle needs to reach a given momentum. For SNR in the Sedov phase (blue lines), the differences due to assumed uniform and nonuniform fluid velocity (cf. the solid and dashed lines) become apparent only for large momenta, p/m_pc ≳ 10⁴. When SNR is in the post-adiabatic phase (green lines), the differences between the two approaches appear for lower momenta, p/m_pc ≳ 10^2.6. These differences affect the low momenta as p/m_pc ≳ 30 for later evolutionary epochs (the yellow and the red lines). There are two effects apparent from this plot. On the one hand, by comparing the lines with each other, we see that as the system enters the radiative phase, the acceleration time for most momenta first increases (blue, green, and yellow lines) and then decreases (red lines). This behavior is primarily caused by the temporal deceleration and subsequent acceleration of the shock, since the same trend is observed in both approaches, for the sequences of solid and the sequence of dashed lines. On the other hand, by comparing a solid with a corresponding dashed line, the flow nonuniformity downstream of the post-adiabatic shock shortens the acceleration time compared to the case of a uniform flow. We consider, for example, the yellow solid and dashed lines and compare the points of their intersection with the yellow dotted horizontal line. This leads to a conclusion that the system with the flow nonuniformity (solid line) is able to accelerate particles to the momentum p/m_pc ≃ 10^3.5 while the system with a uniformity to p/m_pc ≃ 10^2.5 only.

Fig. 12 shows the maximum momenta derived as a solution of equation t = ⟨t_acc(p_max)⟩, on the flow structures relevant for a given time t. Although both approaches (with uniform and nonuniform flow) typically result in similar values of p_max before ≈ 1.5t_tr, the differences may reach an order of magnitude during the post-adiabatic phase, due to impact from the region with negative u₂.

The evolution of the particle distribution function f_o(p) for the two models is shown in Fig. 13 and in animations in Appendix A.2. Similarly to a nonlinear acceleration regime, the spectrum becomes concave, though we neglect the back-reaction of CRs on the shock, and u₁ is uniform. The spectra are so hard because of the quite strong gradients of u₂, large differences between u₁ and u₂, and the presence of the region with u₂ directed toward the shock (Fig. 2), as considered in detail in Sections 3 and 4.

From spectra f_o(p) at different times, we may calculate the evolution of the spectral index of the radio emission generated by electrons accelerated at the forward shock. Fig. 14 shows the evolution of the radio spectral index for the post-adiabatic shock for model 1. It is calculated as described in Appendix D, considering the emission from a thin shell with varying B₂(x) and n₂(x). From Fig. 14 we can see that the radio spectral index decreases starting from time about 1.6t_tr. This is in perfect agreement with expectations and observational hints presented in Sect. 2.2.

In the case of Bohm diffusion (top plot), the higher frequencies respond earlier to radiative cooling because of the higher energies of electrons emitting at these frequencies, so they can diffuse to larger distances from the shock (cf. Fig. 2). In a similar fashion, the sensitivity of the radio index to the structure of the flow is higher for the Kolmogorov diffusion coefficient (bottom plot) because the electrons in this case may probe the deeper regions downstream of the shock.

The radio spectral index at these frequencies is close to 0.5 all the time in model 2, either with the Bohm or Kolmogorov diffusion coefficient. The reason is that in this model, particles become sensitive to the flow structure only at higher momenta (Figs. 7 and 13). By increasing the Kolmogorov diffusion coefficient, the behavior of the radio spectral index in model 2 can be made similar to that in model 1.

Fig. 11 Mean acceleration time ⟨tacc(p)⟩ for our numerical model 1 and Bohm diffusion coefficient. The solid lines show calculations with Eqs. (15) and (1) for xp. The colored lines are calculated with the spatial distributions of MHD parameters relevant to the time moments t shown in the legend. The dashed lines of the same color correspond to the same time moments and uniform flow, Eq. (14) with βa = 6.67. The dotted horizontal lines are added to guide the eye; they correspond to the time moments shown in the legend.

Fig. 12 Time-limited maximum momentum for our MHD models and Bohm diffusion coefficient. As in Fig. 11, the solid lines represent calculations with the spatial structure of the flow, while the dashed lines are for the uniform flow. The three vertical lines mark the times ttr (left), 1.49ttr (when the region with negative u2 appears), and tsf (right).

Fig. 13 Spectra of accelerated particles for models 1 and 2 with Bohm diffusion obtained with equations from Sect. 3.4. The colored lines correspond to the different time moments (see inset).

Fig. 14 Evolution of the radio index for model 1 at a few frequencies shown in the legend, in GHz. Bohm (top) and Kolmogorov (bottom) diffusion are assumed. We set DKolm(p∗) = DBohm(p∗) at p∗ = 55mpc for this plot. The vertical lines mark the time t/ttr = 1.49 (then the negative u2 first appears) and tsf.

7 Discussion and conclusions

The main findings of the present paper can be summarized as follows.

• MHD simulations of SNRs demonstrate essential changes in the flow structure behind the shock on the post-adiabatic stage.

• The changes also affect the profiles of the flow velocity on the length scales of the CR acceleration.

• As a result, the acceleration efficiency increases. This happens due to three reasons: a higher u₁/u₂, an increased ∆p per acceleration cycle, and a possible existence of a region with u₂ < 0.

• Together, these factors make the particle spectrum harder.

• It could also result in a decrease in the radio spectral index as the shock progresses into the post-adiabatic regime.

• MF becomes dynamically important in the post-adiabatic shocks; it affects the profiles of the flow velocity and, therefore, could influence the particle spectrum.

We also report the observed data to support these conclusions: the correlation of the radio index in Kes 73 with the interstellar density and the decrease in the radio index with age in interacting SNRs.

However, the application of our theoretical considerations to 1D numerical simulations of the post-adiabatic SNRs seems to exaggerate the overall effect. There are a few limitations that should be taken into account.

In particular, fluctuations in the shock speed V and the parameter m (Fig. 1) seem to appear because the numerical MHD simulations are 1D. In 3D, the flow has more degrees of freedom, and we may expect, in particular, that the first minimum in V and m around time 1.7t_tr could not be so deep. One might also expect that the spatial profiles of the flow velocity in more realistic 3D simulations could have softer gradients and contrasts. However, simulations of post-adiabatic flows, in order to be coupled to the CR acceleration, require a resolution down to ∼10⁻⁵R around the shock. Such simulations are very demanding in three dimensions, considering the spatial and temporal ranges that must be covered.

Another point is that we used the steady-state solution for particle acceleration on a stationary MHD background. This implicitly assumes that the acceleration time is shorter than the hydrodynamic timescale for spatial changes in the flow. However, once the radiative losses of plasma appear, the structure of the flow changes rapidly (animations in Appendix A.1). Therefore, though the stationary solution could be reasonable for the radio-emitting electrons, the spectra of CRs at the highest energies (emitting X- and gamma-rays) could be altered in a time-dependent treatment which considers nonstationary acceleration of particles on the nonstationary profiles of the flow velocity. The numerical treatment of such a problem, together with high-resolution 3D MHD simulations, is almost unfeasible at the present time.

Post-adiabatic shocks, despite their rather low shock speed, could be efficient CR accelerators. Our simulations show that a large amount of energy is transferred into the high-energy CRs (Fig. 13). Numerically, the energy in CRs (as from an instant f_o(p)) reaches 50% of the shock kinetic energy about time 1.6t_tr and rapidly grows from that time on. If the energy in CRs becomes this high, then their back-reaction on the flow upstream should also be taken into account.

Furthermore, such high efficiency in transferring kinetic energy into the CRs creates a new channel of energy losses of plasma, in addition to the radiative losses. This should affect both the MHD structure of the flow and the shock speed. Thus, the shock decelerates in the post-adiabatic regime not only due to the expansion and losses of thermal radiation but also due to the efficient acceleration of CRs. The influence of CRs on the post-adiabatic plasma has recently been considered by Diesing & Gupta (2025), by introducing a free parameter which sets the ratio of CR pressure to the thermal pressure. The authors have shown that CRs can strongly affect the plasma dynamics at the post-adiabatic phase.

The radiative energy losses, as well as efficient transfer of the shock energy into CRs, should lead to a rapid, successive drop in acceleration efficiency and the disappearance of a post-adiabatic SNR as a nonthermal source. In fact, there is evidence of a fade in radio emission in SNRs shortly after the end of the Sedov stage, which comes from an analysis of the Σ − D relation (Bandiera & Petruk 2010).

When this paper was under review, a publication has appeared at arXiv exploring the CR spectra at shocks with radiative losses (Cristofari 2025). The author has also demonstrated that the shock during the early stage of the radiative phase could be a significantly efficient CR accelerator.

A more refined future treatment of the problem should involve multi-dimensional MHD simulations coupled with particle acceleration.

Data availability

Four movies associated to Appendix A are only available at https://www.aanda.org

Acknowledgements

We thank the referee, Pat Slane, for constructive comments. We acknowledge Yang Chen for providing us with the ¹²CO (J = 1 → 0) maps. We acknowledge the support from INAF 2023 RS4 Theory grant as well as INAF 2023 RS4 mini-grant. OP and TK thank the Armed Forces of Ukraine for providing security to perform this work. This research used an HPC cluster at DESY, Zeuthen, Germany, and the HPC system MEUSA at INAF-Osservatorio Astronomico di Palermo, Italy.

References

Appendix A Animations

A.1 Evolution of the flow downstream of the post-adiabatic shock

The animations demonstrate the time development of the profiles of density (upper right), flow speed v(r) in the laboratory frame (lower left), and magnetic field (lower right) for two models of SNRs described in Sect. 2.1. Namely, Animation 1 for model 1 with parallel shock and B_o = 5 µG (the associated movie is available online) and Animation 2 for model 2 with perpendicular shock and B_o = 10 µG (the associated movie is available online). We note that the horizontal axis in the plot for the flow speed spans only 15% of the shock radius, to emphasize changes around the region where the shock particle acceleration happens. The upper left plot shows the expansion parameter and the moving pointer, which indicates the time moment for which the hydrodynamic profiles have been shown.

A.2 Evolution of the spectrum of particles accelerated at the post-adiabatic shock

Animation 3 for model 1 (the associated movie is available online) and Animation 4 for model 2 (the associated movie is available online) show the evolution of the flow velocity in the reference frame of the shock u₂(x), the spectrum of accelerated particles f_o(p) and the spectral index s(p) = −d ln f_o/d ln p.

Appendix B Some considerations regarding Kes 73

By assuming the uniform expansion and considering the Sedov shock in a uniform ambient medium, Borkowski & Reynolds (2017) have derived the average shock speed in Kes 73 V ≃ 1400 km/s, the average radius R = 6 pc, the average ambient number density n_o ≃ 2 cm⁻³, the explosion energy 1.9 foe and the age t = 1800 yr. As we see (Fig. 14), the radio spectral index is expected to begin decreasing due to radiative losses before 1.6t_tr. We need to determine what ambient density n_mc would allow for a part of the SNR shock to be at the stage corresponding to 1.6t_tr with $t_{\text{tr}}=2.83\times {10}^4{E}_{51}^{4/17}{n}_1^{-9/17}\text{yr}$. Writing t = 1.6t_tr and substituting the explosion energy into the equation for the transition time, we derive the density n_mc ≃ 580 cm⁻³. This value could naturally be in a molecular cloud near Kes 73. Indeed, the column density of material SNR interacts with on the West is $N_{{\text{H}}_2}\simeq {10}^{22}{\text{\hspace{0.17em}cm}}^{-2}$ (Liu et al. 2017). The size of a cloud with such $N_{{\text{H}}_2}$ and n_mc is ≃ 5.6 pc along the line of sight. From Fig. 5 the thickness of a cloud filament in the plane of the sky is about the SNR radius, that perfectly matches the estimated size of the cloud along the line of sight. With such density, the shock speed in the interaction region drops from the average 1400 km/s to the local value 1400 (n_o/n_mc)^1/2 ≈ 82 km/s. Then the post-shock plasma temperature is about $T=14V_{\text{km/s}}^2\simeq 9\times {10}^4\text{\hspace{0.17em}K}$, which corresponds to the maximum of the cooling curve Λ(T). And finally, we determine the area of the shock surface interacting with such dense material. It is the ratio of the area of the spherical cap to the area of a sphere, $$\text{Ω}=\frac{2\pi R^2}{4\pi R^2}\left(1-\sqrt{1-{\left(\frac{\varrho }{R}\right)}^2}\right)\simeq 0.03,$$(B.1)

where ϱ is the average radius of the red area on the spectral index map. It appears that the interacting region is a quite small fraction of the SNR surface. In addition, this region is projected on the SNR interior. Therefore, it does not significantly affect the estimates of the average values derived by Borkowski & Reynolds (2017). All these considerations support our idea that the radio spectral index variation across Kes 73 could mark the onset of radiative losses in the interacting part of the SNR. Nevertheless, this idea remains a hypothesis until a dedicated study. In particular, such a study would help discriminate it from other potential effects. For example, the shock compression can be sufficiently high in radiative shocks to produce prominent free-free emission in the radio band, which would flatten the spectrum.

Appendix C Shape of the numerical CR spectrum

The numerical solution for the momenta p ≳ p_A, as measured from Fig. 8, is nearly a power law p⁴ f_∗(p) ∝ p^a with a from 1.7 for t = 120 yr to 1.2 for t = 1000 yr. This corresponds to the spectral index s_∗ = 4 − a (for the function f_∗(p) ∝ p^−s∗) in the range 2.3 ÷ 2.8. The analytical solution gives s = 3/(1 − u_AB/u₁) = 1 between p_a and p_B. Thus, we have s_∗ − s ≈ 1.3÷1.8, and we need to know whether the numerical and analytical solutions correspond to each other.

The CR spectra from the numerical solution (Fig. 8) are integrated over the entire volume. One should take into account differences between the analytical (Sect. 3.4) and the numerical (Sect. 5.1) treatment of the CR spectra. The latter is a nonstationary solution for expanding spherical shock, and the former, instead, is a steady-state solution for a plane-parallel shock on an immediate MHD background.

The volume-integrated spectrum represents particles injected at different times. The injection efficiency ξ (a fraction of the density) does not change over time in both approaches. However, in the numerical solution, the number of injected particles increases over time due to the increase in the spherical shock surface ∝ (R(t)/R₀)². The number density of particles injected at some time moment t_i decreases over time t as (R(t)/R(t_i))⁻³ because the volume increases. So, the number density of particles accelerated until the time t to some momentum p is smaller than the number density of particles freshly injected at this time t. The shock surface at the time t_i when particles were injected was smaller, and the volume has increased since then. As a consequence, the slope of the numerical spectrum is softer: there are more particles with low momenta just because of the expanding spatial domain. To compare with the analytical solution, we should apply a correction to the numerical solution by making the early injection artificially larger: ξ(t_i) → ξ · [R(t)/R(t_i)]⁵ where t > t_i.

The application of this correction is not straightforward. The CR spectrum at a given time t covers a wide range of momenta. Particles that have some momentum p were injected at different times t_i. However, if we assume that all particles with momentum p spend the same time for acceleration to this momentum, then we may write an approximate relation $$f_{\text{o}}(p)\simeq f_{*}(p,t)\cdot {\left[R(t)/R(t_{\text{i}}(p))\right]}^5,$$(C.1)

where f_o(p) and f_∗(p) are the analytical and numerical distribution functions, R(t) = Vt, t_i(p) = t − t_acc(p). This transforms to a relation for the spectral indices $$s(p)\simeq s_{*}(p,t)-\text{Δ}s(p,t),\text{Δ}s=\frac{5}{3}{\left[\frac{t}{t_{\text{acc}}(p)}-1\right]}^{-1},$$(C.2)

where t_acc ∝ p^1/3 is used, s and s_∗ are indices for f_o(p) and f_∗(p) respectively. We note that t ≥ t_acc and, therefore, ∆s ≥ 0 as needed. If t_acc ≪ t, then ∆s → 0. On the other hand, if t_acc → t, then ∆s could be quite large. For intermediate ratios t/t_acc ∼ 2, the difference ∆s ∼ 5/3 is close to the differences s_∗ − s we observed in our results.

Appendix D Calculation of the radio spectral index

Synchrotron emissivity is $$P(\nu )=\underset{0}{\overset{\infty }{{\displaystyle \int }}}N(p)\epsilon (p,\nu )dp,$$(D.1)

where p is the electron momentum, N(p) = 4πp² f_o(p) is the electron momentum distribution, ν is the synchrotron frequency, ϵ is the single-electron emissivity $$ϵ=c_1\mu BF(\nu /{\nu }_c),{\nu }_c=c_2\mu B{p}^2,$$(D.2)

where $c_1=\sqrt{3}{e}^3/(m_{\text{e}}{c}^2),c_2=3e/(4\pi m_{\text{e}}^3{c}^3)$. We take the geometrical factor µ = 1, i.e., MF is perpendicular to the line of sight. We use the delta-function approximation for the special function $$F(x)=\frac{8\pi }{9\sqrt{3}}\delta (x-0.29).$$(D.3)

In this approximation $$P(\nu )=c_3{\nu }^{1/2}{B}^{1/2}N(p_{\text{eff}}),p_{\text{eff}}={\left(\frac{\nu }{0.29c_2\mu B}\right)}^{1/2}$$(D.4)

where c₃ is a constant, p_eff is the momentum of the electron that emits most of its energy at frequency ν in the magnetic field B. For example, this approximation yields the well-known dependence P ∝ B^(s+1)/2ν^−(s−1)/2 for N ∝ p^−s.

Most of the emission comes from a shell of the thickness about 0.1R, where MHD parameters vary in space. The flux density is therefore proportional to an integral of emissivity P over this shell, $$S(\nu )\propto {\nu }^{1/2}{\displaystyle \underset{0}{\overset{0.1R}{\int }}{B}_2{(x)}^{1/2}{n}_2(x)f_{\text{o}}(p_{\text{eff}}(x,\nu ))p_{\text{eff}}{(x,\nu )}^{2}dx},$$(D.5)

where p_eff depends on x through B₂(x). The spatial profile of advected particles downstream is approximated as f (x, p) ≈ n₂(x) f_o(p)/n₂(0) (Orlando et al. 2007, equation A.6).

The radio spectral index at the frequency v₀ is determined from fluxes at two frequencies: $$\alpha ({\nu }_0)=-\frac{\lg [S({\nu }_1)/S({\nu }_2)]}{\lg ({\nu }_1/{\nu }_2)}.$$(D.6)

We choose these frequencies as $$\frac{{\nu }_2}{{\nu }_0}=\frac{{\nu }_0}{{\nu }_1}=\mathrm{1.05.}$$(D.7)

We note that a simpler approach, namely, α = (s − 3)/2 where the index s of the electron momentum distribution is measured at momentum p_eff(0, ν₀), results in a similar temporal variation of the radio spectral index as in a more complex and accurate approach described above.

Appendix E Erratum

We noticed misprints in the text of Appendix C in the paper Petruk & Kuzyo (2024).

In the sentence between equations (C.1) and (C.2), there should be “ f (x, p) = 0 at x = −∞ and ∂f /∂x = 0 at x = ±∞” instead of “ f (x, p) = 0 at x = ±∞ and ∂f /∂x = 0 at x = −∞”.

Expression (C.7) should read as follows: $${\left[D\frac{\partial f}{\partial x}\right]}_{+\infty }-u_3{f}_3\equiv -u_{\text{x}}{f}_{\text{o}}.$$(E.1)

All Figures

	Fig. 1 Top: Evolution of the radius R (dashed lines) and speed V (solid lines) of the forward shock. Bottom. Time dependence of the expansion parameter m. The three vertical lines mark the times ttr (left), 1.49ttr (middle), and tsf (right).
In the text

Fig. 2 Spatial distributions of the flow speed u2, magnetic field strength B2, and number density n2 downstream of the shock at several moments of time. The left column is for model 1 and the right column for model 2. The coordinate x marks the distance from the shock, which is used for considerations of the CR acceleration. We note that the horizontal axes are in the log scale. This is necessary to describe equally well the variations on both small and large scales relative to the shock position, which are probed by particles with momenta spanning several decades. The shock is located on the left of each plot, at x/R = 0. The SNR center is on the right, at the point with lg(x/R) = 0.

In the text

	Fig. 3 Distances xA and xB from the shock for the two models. The flow speed u2 < 0 between the two solid lines. The vertical dotted lines mark the times 1.49ttr and tsf. The flow speed u2 > 0 in the whole domain before 1.49ttr.
In the text

Fig. 4 Electron spectral index vs. age for a sample of SNRs with the radio and GeV γ-ray emission detected, from the data given by Zeng et al. (2019). The red dashed line represents a linear fit between the index and logarithm of age, with the slope k = −0.37 ± 0.05. The red area shows the 1σ uncertainty. We excluded MSH 15-56 from the sample because its radio emission is likely dominated by a pulsar. There are a few SNRs with two values of α estimated in this reference. We calculated the average α from the two values if the uncertainty is due to various estimates of the ambient density (Tycho, G166.0+4.3, S147, CTB 73B). We use only one value for W51C, RX J0852-4622, namely that which corresponds to the same spectral model as used for other SNRs.

In the text

Fig. 5 Distribution of the radio spectral index over Kes 73 measured between 1.4 and 5.0 GHz (left) and the 12CO (J = 3 → 2) intensity map (Dempsey et al. 2013) integrated over the velocity range between 94 and 95 km/s, in units of K km/s (right). The contours correspond to the values of the radio index 0.8, 0.6, and 0.4. Liu et al. (2017) also present 12CO (J = 1 → 0) intensity maps around Kes 73 at various velocities. Those in the range 86–97 km/s correlate with the radio index map as well.

In the text

Fig. 6 Spatial structure of the flow velocity (blue line) in the reference frame of the shock, which is located at x = 0. The diffusion length of a particle with momentum p is xp1 upstream and xp2 downstream. These are the largest distances the particle can diffuse and return to the shock to be further accelerated. The shock compression factor σ = u1/u2 is the ratio of the immediate pre- to post-shock velocities. Instead, σp(p) = up1(p)/up2(p) is an effective compression factor “seen” by a relativistic particle with momentum p; it varies with momentum. The shock here is not CR-modified, and therefore up1(p) = u1 upstream, for any momentum. For the downstream region, one may use an approximation up2(p) ≈ u2(xp2) (see Sect. 3.4).

In the text

Fig. 7 Minimum momentum p of particles that may reach the distance xp downstream and return to the shock to continue acceleration. It is calculated from the equation $p=\overline{p}(x_{\text{p}})$ for Bohm diffusion. The colored lines correspond to a few time moments (see inset). The dashed lines show p calculated in a common approach, from xp = D2(p)/u2 with uniform u2 and uniform B2 = 5 µG. The top plot is for model 1 and the bottom plot is for model 2.

In the text

	Fig. 8 Volume-integrated spectra for particles inside SNRs at different times (marked in the plot legend, in years). The solid lines represent the model with parameters uAB = −2u1, xA/R = 0.01, xB/R = 0.03. The spectra shown by the dashed lines are calculated for the uniform flow speed.
In the text

Fig. 9 Diffusion length xp2(p) (top) and mean acceleration time ⟨tacc(p)⟩ (bottom) calculated in the analytical approach (i.e., it is time-independent). The model with xA = 0.01R, xB = 0.03R, uAB = −2u1 is shown with the orange line; it is calculated with Eqs. (1) and (15). The model with uniform u2 is shown with the blue line, as from xp = D2/u2 (then simply xp ∝ p1/3) and (14) for acceleration time. The green dashed line shows xp1(p) upstream, and the gray dashed line corresponds to t = 420 yr.

In the text

Fig. 10 Maximum momenta pmax for two numerical solutions (dots) and from analytical expressions (solid lines). The blue line is from Eq. (14) and the dots are from the numerical solution with uniform u2 (dashed lines in Fig. 8). The orange line is calculated with Eq. (15) with xp1 and xp2 from (1) and the dots are from the numerical solution with xB = 0.03R, uAB = −2u1 (solid lines on Fig. 8). The value of pmax is measured at the level where the distribution function drops below the power-law extrapolation in e times for the orange dots and in 10 times for the blue dots.

In the text

Fig. 11 Mean acceleration time ⟨tacc(p)⟩ for our numerical model 1 and Bohm diffusion coefficient. The solid lines show calculations with Eqs. (15) and (1) for xp. The colored lines are calculated with the spatial distributions of MHD parameters relevant to the time moments t shown in the legend. The dashed lines of the same color correspond to the same time moments and uniform flow, Eq. (14) with βa = 6.67. The dotted horizontal lines are added to guide the eye; they correspond to the time moments shown in the legend.

In the text

	Fig. 12 Time-limited maximum momentum for our MHD models and Bohm diffusion coefficient. As in Fig. 11, the solid lines represent calculations with the spatial structure of the flow, while the dashed lines are for the uniform flow. The three vertical lines mark the times ttr (left), 1.49ttr (when the region with negative u2 appears), and tsf (right).
In the text

	Fig. 13 Spectra of accelerated particles for models 1 and 2 with Bohm diffusion obtained with equations from Sect. 3.4. The colored lines correspond to the different time moments (see inset).
In the text

	Fig. 14 Evolution of the radio index for model 1 at a few frequencies shown in the legend, in GHz. Bohm (top) and Kolmogorov (bottom) diffusion are assumed. We set DKolm(p∗) = DBohm(p∗) at p∗ = 55mpc for this plot. The vertical lines mark the time t/ttr = 1.49 (then the negative u2 first appears) and tsf.
In the text