Labels

Thursday, October 11, 2018

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


No comments:

Post a Comment