The solution of cosmic ray propagation equation

transport equation

$$\nabla  \cdot (D \nabla N_e) + \frac{\partial}{\partial E}(\frac{dE}{dt}N_e) + Q_e = \frac{\partial N_e}{\partial t}$$

energy change of particle
$$E(t+dt) = E(t) + \frac{dE}{dt}dt$$
$$dt = \frac{dE}{dE/dt}$$
$$\int_{t_1}^{t_2} dt = \int_{E_1}^{E_2}  \frac{dE}{dE/dt}$$

if $dE/dt \equiv -b(E) = -aE^2$  and $t_2>t_1$
$$t_2 - t_1 = \frac{1}{a}[ \frac{1}{E_2} - \frac{1}{E_1}]$$
$$E(t_1) = \frac{E(t_2)}{1-a(t_2-t_1)E(t_2)}$$
and
$$E(t_2) = \frac{E(t_1)}{1+a (t_2-t_1)E(t_1)}$$

green function
$$\begin{equation}\nonumber  G_e(\vec{r},\vec{r}',E,E',t,t') = \frac{1}{[\pi 4 f(E,E')]^{3/2}b(E)} \exp[\frac{-(\vec{r}'-\vec{r})^2}{4f(E,E')}]\cdot \delta(t'-t+h(E,E')) \end{equation}$$

$f(E,E') = \int_{E'}^E  D(u)\frac{du}{du/dt}$ and $h(E,E') = \int_{E'}^E \frac{du}{du/dt}$

$\delta$ function
$$\delta(g(x)) = \sum_i  \frac{\delta(x-x_i)}{|g^\prime(x_i)|}$$
$t>t'$    the time axis  ------t'---------t----->
$\delta(t' - t + h)$
$$h = \int_{E'}^{E}  \frac{dE}{dE/dt} = - \frac{1}{a} [ \frac{1}{E^\prime} - \frac{1}{E} ]$$

$g(E') = t'-t+h = t'-t + \frac{1}{a}[\frac{1}{E}-\frac{1}{E'}]$

$g(E') = t'-t + \frac{1}{a}[\frac{1}{E}-\frac{1}{E'}] = 0$

$g'(E') = \frac{1}{a{E^\prime}^2}$

$$\begin{equation} \delta(t' - t + h)  =  \delta(t' - t  + [- \frac{1}{a}  (\frac{1}{E^\prime} - \frac{1}{E} )  ])                                  =  \delta(E' -\frac{E}{1 -a (t-t')E}) \cdot [a (\frac{E}{1-a(t-t')E})^2] \end{equation}$$

The solution of transport equation is given by
$$\begin{equation}\nonumber \begin{split}  N_e(\vec{r},E,t)  &= \int \int \int d\vec{r}' dE' dt' G_e(\vec{r},\vec{r}',E,E',t,t')Q_e(\vec{r}',E',t') \\                               &= \int dt' \int d\vec{r}'  \frac{b(E')}{b(E)} \frac{Q_e(\vec{r}',E',t')}{(4\pi f(E,E'))^{3/2}} \exp{[\frac{-(\vec{r}' - \vec r)^2}{4f(E,E')}]} \end{split} \end{equation}$$
In this equation, the $E'$ is given by $E' =  \frac{E}{1-a(t-t')E}$.


In fact, the solution $$N(\vec{r},E,t)  = \int dt' \int d\vec{r}'  \frac{b(E')}{b(E)} \frac{Q(\vec{r}',E',t')}{(4\pi f(E,E'))^{3/2}} \exp{[\frac{-(\vec{r}' - \vec r)^2}{4f(E,E')}]}$$
is general. It doesn't  depend on the assumption of the form of energy loss.
if we expand $\delta(t'-t + \int_{E'}^{E} \frac{dE}{dE/dt})$ directly, we can get the general solution. In this case $E'$ is a function of $E$ and $t-t'$.




The integration of spatial part
To simplify our mark, we let the $\vec r =0$. It means we take the position of observer as original point.
if the spatial part of source term $Q$ is $q(r')$, the integration for spatial part is
$$SI = \int d\vec{r}' q(\vec r') \exp( \frac{-{r^\prime}^2}{4f(E,E')}) =  \int q(\vec{r}')  \exp( \frac{-{r^\prime}^2}{4f(E,E')}) {r^\prime}^2dr'd\Omega$$

Assume the source locates at $\vec r_s$ and the position with respect to source is $\vec r_0$ .
For the geometry, one can refer the figure 1 in the paper https://arxiv.org/pdf/1009.0894.pdf.

Then $\vec r' = \vec r_s + \vec r_0$ and ${r^\prime}^2 = r_s^2 + r_0^2 - 2r_s r_0 \sin(\theta)\cos(\phi)$.

Note that the $\vec r_s$ is fixed,  $f(\vec r') = f(\vec r_0)$ and $d\vec r' = d\vec r_0$.

If the particles are homogeneous injected by the sphere surface with radius $R$, the source term becomes $q(\vec r') = \frac{1}{4\pi R^2}\delta( r_0 -R)$.
I mark $4f(E,E')$ as L.
$$\begin{equation}\nonumber \begin{split}  SI \cdot 4\pi R^2 &= \int  \exp(-[r_s^2 + r_0^2 - 2r_s r_0 \sin(\theta)\cos(\phi)]/L) r_0^2 \delta(r_0-R)dr_0 \sin(\theta)d\theta d\phi \\                               &= \int R^2 \exp( -[r_s^2 + R^2]/L)\cdot  \exp(2r_s R \sin(\theta)\cos(\phi)/L) \sin(\theta)d\theta d\phi \end{split} \end{equation}$$
let $2r_s R/L = c$.

how to compute the integration $\int_0^{2\pi}\int_0^{\pi}  \exp(c \sin(\theta)\cos(\phi)) \sin(\theta)d\theta d\phi$  ?

The author of paper On the point-source approximation of nearby cosmic ray sources says the integration is $4\pi/c \sinh(c)$. I checked this result by numerical integration with python.

#for example

import scipy.integrate as integrate
from math import *
def f(x,y):
    return exp(-10*sin(x)*cos(y))*sin(x)
integrate.nquad(f,[[0,pi],[0,2*pi]]) 
#output (13839.636596576718, 4.0309916130354395e-05)
4*pi/10*sinh(10)
#output 13839.63659657671

$$\int_0^{2\pi}\int_0^{\pi} \exp(c \sin(\theta)\cos(\phi) = \frac{4\pi}{c}\sinh(c)$$
$$SI = \frac{2}{r_s R} \exp(-\frac{r_s^2+R^2}{4f(E,E')}) \sinh(\frac{r_s R}{2f(E,E')}$$
$$\begin{equation}\nonumber \begin{split} N(\vec{r},E,t) &= \int dt' \frac{b(E')}{b(E)} Q(t',E') \frac{1}{4\pi^{3/2}f^{1/2}r_s R} \exp(-\frac{r_s^2+R^2}{4f(E,E')})\sinh(\frac{r_sR}{2f(E,E')}) \end{split} \end{equation}$$
$\int_{t'}^{t} dt = \int_{E'}^{E}  \frac{dE}{dE/dt}$


REFERENCE
solution
On the point-source approximation of nearby cosmic ray sources 
numerical integration
https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html
Dirac function
https://en.wikipedia.org/wiki/Dirac_delta_function

Others
On the contribution of nearby sources to the observed cosmic-raynuclei