Cosmic-Ray Energy-Dependent Injection (CREDIT)

• CREDIT is a non-stationary cosmic-ray model where injection time depends on particle energy or rigidity, shaping post-supernova acceleration.
• The framework links rigidity-dependent escape from SNRs with early high-energy particle release and subsequent spectral breaks.
• By incorporating time-varying injection and stochastic source populations, CREDIT models predict narrow, structure-enhanced features in local proton spectra.

Cosmic-Ray Energy-Dependent Injection Time (CREDIT) denotes a class of non-stationary source models in which the time relevant to cosmic-ray release or entry into acceleration depends on particle energy or rigidity. In supernova-remnant (SNR) escape models, CREDIT means that higher-rigidity particles are released earlier and lower-rigidity particles later; in a non-stationary diffusive-shock-acceleration treatment, it means that the highest-energy particles observed at a later epoch must have entered acceleration at earlier times after the supernova explosion. The term should therefore be distinguished from work on the source spectral exponent alone, where “injection” refers to the power-law index of the emitted spectrum rather than to an energy-dependent time history (Stall et al., 2024, Petruk et al., 2017, Lagutin et al., 2017).

1. Terminological scope and principal usages

The cited literature uses the term “injection” in distinct but related ways. In all cases, the central issue is temporal non-stationarity, but the time variable attached to particle production or release is not identical across models.

Reference	Energy/time content	Relation to CREDIT
(Petruk et al., 2017)	Higher-energy particles observed today must have been injected earlier into acceleration; gamma-ray energy maps to post-explosion epoch	Explicit CREDIT concept
(Stall et al., 2024, Stall et al., 9 Jul 2025)	Higher-rigidity cosmic rays escape an SNR earlier; lower-rigidity particles escape later	Explicit CREDIT scenario
(Merten et al., 10 Jun 2026)	Continuous injection with energy-dependent effective escape from an evolving wind bubble	CREDIT-relevant, but not a full explicit CREDIT law
(Lagutin et al., 2017)	Reconstruction of the source spectral exponent $p$ from Galactic propagation	Not about energy-dependent injection time

This terminological distinction is essential. In (Stall et al., 2024) and (Stall et al., 9 Jul 2025), CREDIT is a rigidity-dependent escape history from individual SNRs. In (Petruk et al., 2017), CREDIT links gamma-ray energy to the epoch when parent particles entered the acceleration process. By contrast, (Lagutin et al., 2017) assumes an instantaneous point source with a power-law injection spectrum and studies how propagation through an inhomogeneous interstellar medium modifies the observed spectrum, not how injection time depends on energy.

2. CREDIT as a non-stationary acceleration history in supernova remnants

A formulation of CREDIT appears in the analysis of gamma-ray spectra from SNRs, where the highest-energy hadrons require the longest acceleration time and therefore preserve a memory of the earliest post-explosion injection history. In this picture, low-energy particles and low-energy gamma rays are mainly shaped by later injection, whereas TeV gamma rays originate from particles that began acceleration during the first months after the supernova explosion. The particle transport is modeled with the non-stationary acceleration equation

$$\frac{\partial f}{\partial t}+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,$$

and the injection term

$$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$

where $Q_t(t)$ carries the time dependence of the injection efficiency (Petruk et al., 2017).

In the test-particle regime and $t_1 \gg t_2$, the shock solution is written as

$$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$

with stationary spectrum

$$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$

The steady-state criterion is

$$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$

For a power-law injection history, $Q_t\propto \tau^\beta$, the spectrum becomes approximately

$s=s_f+\alpha\beta,$

with $$\frac{\partial f}{\partial t}+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,$$0. Time-dependent injection therefore changes not only the normalization but also the slope of the accelerated particle spectrum.

The application to IC443 uses SN1987A as a proxy for the parent supernova. From ATCA radio data between 1517 and 8014 days after explosion, the inferred injection history rises smoothly with time and is approximately described during much of the observational interval by

$$\frac{\partial f}{\partial t}+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,$$1

reaching the steady-state level $$\frac{\partial f}{\partial t}+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,$$2 only around $$\frac{\partial f}{\partial t}+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,$$3 days. For the unobserved early phase, the adopted model assumes steady injection up to day 100, then $$\frac{\partial f}{\partial t}+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,$$4 with $$\frac{\partial f}{\partial t}+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,$$5. With the same model, the fit to the IC443 gamma-ray spectrum reproduces the observed data from Fermi, MAGIC, and VERITAS without imposing an ad hoc proton break: the proton spectrum breaks around $$\frac{\partial f}{\partial t}+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,$$6, the gamma-ray spectrum steepens around $$\frac{\partial f}{\partial t}+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,$$7 GeV, the spectrum above $$\frac{\partial f}{\partial t}+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,$$8 TeV is controlled by the earliest injection phase, photons above $$\frac{\partial f}{\partial t}+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,$$9 GeV are influenced by injection before the radio observations of SN1987A began, and the gamma-ray spectrum above $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$0 GeV reflects the first $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$1 years of SNR evolution (Petruk et al., 2017).

3. CREDIT as rigidity-dependent escape from supernova remnants

A second, and now more explicit, formulation of CREDIT treats the source term itself as rigidity dependent in time. In this SNR escape picture, cosmic rays do not all leave the remnant simultaneously; instead, the escape time depends on rigidity $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$2. The source term is

$$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$3

where $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$4 is the rigidity-dependent escape time. The corresponding smooth-source null hypothesis averages over many remnants and yields

$$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$5

In that limit, the escape history is invisible because source discreteness is washed out (Stall et al., 2024, Stall et al., 9 Jul 2025).

The physical picture is tied to the evolution of the maximum confining rigidity. Early in the remnant’s life, magnetic-field amplification is strongest and higher maximum rigidities can be confined. After the Sedov-Taylor transition, the shock slows and $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$6 drops. The escape delay is parameterized as

$$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$7

with

$$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$8

for $$Q(t,x,p)=\frac{\eta n_1 u_1}{4\pi p_i^2}\,\delta(p-p_i)\,\delta(x)\,Q_t(t),$$9, and

$Q_t(t)$0

otherwise. The fiducial values are $Q_t(t)$1, $Q_t(t)$2, $Q_t(t)$3, $Q_t(t)$4 as the physically motivated confinement scale, and source rate $Q_t(t)$5. The injection spectrum is

$Q_t(t)$6

This construction interpolates between CREDIT and burst-like injection. For $Q_t(t)$7, escape is energy dependent; in the limit $Q_t(t)$8, all rigidities are released at the same time. The CREDIT prediction is that young and nearby SNRs produce rigidity-localized enhancements in the local proton spectrum, rather than a broad smooth tilt, because a source contributes sharply when its age and $Q_t(t)$9 align (Stall et al., 2024, Stall et al., 9 Jul 2025).

4. Galactic transport, Green’s functions, and stochastic source populations

The rigidity-dependent SNR CREDIT scenario is embedded in a simplified Galactic diffusion model. The transport equation for protons is

$t_1 \gg t_2$0

with diffusion coefficient

$t_1 \gg t_2$1

The setup assumes that all sources lie in the Galactic disk $t_1 \gg t_2$2, diffusion occurs in a halo with free-escape boundaries at $t_1 \gg t_2$3, no radial boundary is imposed, and only protons above a few GV or above about 10 GeV are considered, so inelastic losses, advection, reacceleration, and convection are neglected in the relevant implementations. Fiducial parameters are $t_1 \gg t_2$4, $t_1 \gg t_2$5 at $t_1 \gg t_2$6, and halo height $t_1 \gg t_2$7 (Stall et al., 2024, Stall et al., 9 Jul 2025).

For a single source at $t_1 \gg t_2$8, the source term is

$t_1 \gg t_2$9

and the Green’s-function solution is written with

$$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$0

The factor $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$1 enforces the free-escape boundary conditions by a mirror-charge construction, and the local flux is

$$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$2

A causality restriction is imposed by excluding source injections outside the observer’s past light cone.

The source population is stochastic. Ages are drawn from a uniform distribution consistent with the Galactic supernova rate, source positions are drawn from an axisymmetric disk distribution depending on Galactocentric radius, the Sun is placed at $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$3, the Galaxy radius is $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$4, and sources are treated as one SNR population injecting protons up to the knee. The break rigidity $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$5 is either fixed within realizations spanning $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$6 or drawn per source from

$$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$7

The spectra are generated by summing Green’s-function contributions from tens of millions of simulated sources, and one implementation uses GPU acceleration via jax (Stall et al., 2024, Stall et al., 9 Jul 2025).

5. Predicted spectral structure and statistical discrimination

The principal observational prediction of the SNR CREDIT scenario is the appearance of narrow rigidity-dependent peaks in the local proton spectrum. Below $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$8, realizations are fairly smooth, with only percent-level deviations from the ensemble mean. Above $$f_o(t,p)=f_o(p)\int_0^\tau Q_t(\tau-\tau')\,\varphi_o(\tau')\,d\tau',$$9, the CREDIT scenario generates pronounced, narrow peaks and rapid structure. For young and nearby sources, the flux can be enhanced by tens of percent in narrow rigidity intervals, and in some realizations even by a few times the average flux in narrow bins. These structures remain visible after rebinning to AMS-02 or DAMPE rigidity bins and exceed expected experimental uncertainties in many cases (Stall et al., 2024, Stall et al., 9 Jul 2025).

The uncertainty model used for detectability estimates approximates statistical errors as

$$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$0

with a $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$1 uncorrelated systematic uncertainty added in quadrature. The key contrast is with the smooth-source null hypothesis, which assumes a continuous spatial source density, steady injection, and no observable imprint of individual source ages or distances. In that limit, the local spectrum is effectively the ensemble mean.

To test distinguishability, a decision tree classifier is trained on three hypotheses: smooth source distribution with uncorrelated experimental errors, burst-like discrete source injection, and CREDIT with $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$2. The classifier takes flux values in rigidity bins above $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$3 as a 58-dimensional input space. It is trained and validated with $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$4 realizations per scenario and 10-fold cross-validation. The confusion matrices are almost diagonal, misclassifications are rare, and varying $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$5 either fixed per realization or randomly per source barely changes the performance. Halving the source rate $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$6, halving the source lifetime $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$7, or smoothing the release time with a Gaussian kernel of width $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$8 does not substantially alter the classification structure for $$f_o(p)=\frac{\eta n_1}{4\pi p_i^3}\frac{3\sigma}{\sigma-1}\left(\frac{p}{p_i}\right)^{-s_f}, \qquad s_f=\frac{3\sigma}{\sigma-1}.$$9 (Stall et al., 2024, Stall et al., 9 Jul 2025).

The interpretive use is explicit. A CREDIT-like classification could be used to infer or constrain the supernova time $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$0 of the source associated with a feature, given assumptions about $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$1, $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$2, $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$3, and $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$4. A burst-like classification would imply a lower bound on $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$5. A smooth classification would challenge the SNR paradigm in the sense stated in the paper: non-detection of source signatures would be unlikely if SNRs are indeed the dominant Galactic cosmic-ray sources (Stall et al., 2024).

6. Related transport problems and non-CREDIT usages

A related but distinct time-dependent escape problem arises in wind bubbles. There, particles are continuously injected at the wind termination shock and propagate through advection and diffusion until escape at the time-dependent position of the forward shock, treated as a free escape boundary. The model is one-dimensional and spherically symmetric, with

$$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$6

so the distance between termination shock and forward shock increases with time. The injected spectrum is $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$7 with a time-dependent high-energy cutoff $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$8. There is no explicit energy-dependent injection-time formula, but the effective escape time becomes energy dependent because high-energy particles diffuse faster and are more likely to reach the moving boundary, while low-energy particles are delayed or trapped. The escaping spectra can be harder than $$\mathcal{I}(\tau)\equiv \int_0^\tau \varphi_o(\tau')\,d\tau' \approx 1.$$9; for Kraichnan diffusion the average hardening is about $Q_t\propto \tau^\beta$0, Bohm diffusion gives the strongest low-energy suppression, and Kolmogorov diffusion the weakest. The trapped low-energy population may contribute to multimessenger radiation and accumulated grammage within the bubble (Merten et al., 10 Jun 2026).

A separate source of confusion is the use of “injection” in studies of the source spectral exponent. In the analysis of the average cosmic-ray injection exponent at Galactic sources, the source is assumed to be instantaneous and point-like,

$Q_t\propto \tau^\beta$1

with diffusion coefficient

$Q_t\propto \tau^\beta$2

In the homogeneous picture,

$Q_t\propto \tau^\beta$3

so if $Q_t\propto \tau^\beta$4 one recovers $Q_t\propto \tau^\beta$5, which for $Q_t\propto \tau^\beta$6 and $Q_t\propto \tau^\beta$7 gives $Q_t\propto \tau^\beta$8. The paper’s main result is that, once interstellar-medium inhomogeneity is included, both anomalous diffusion and normal diffusion in a non-homogeneous medium yield a steeper average source index, $Q_t\propto \tau^\beta$9. This is a result about the average energy spectral exponent at Galactic sources, not about a source model in which higher- or lower-energy cosmic rays are injected at different times (Lagutin et al., 2017).