An analytically iterative method for solving problems of cosmic-ray modulation

TL;DR

This paper introduces an analytically iterative method for solving cosmic-ray modulation problems, effectively approximating solutions with minimal iterations and matching well with existing analytical and numerical results.

Contribution

The paper presents a new iterative approach for solving cosmic-ray modulation equations, incorporating anisotropy and providing accurate solutions with fewer iterations.

Findings

The method accurately reproduces analytical solutions for constant diffusion coefficients.

Two iterations suffice to match numerical solutions closely.

The approach effectively handles different forms of diffusion coefficients.

Abstract

The development of an analytically iterative method for solving steady-state as well as unsteady-state problems of cosmic-ray (CR) modulation is proposed. Iterations for obtaining the solutions are constructed for the spherically symmetric form of the CR propagation equation. The main solution of the considered problem consists of the zero-order solution that is obtained during the initial iteration and amendments that may be obtained by subsequent iterations. The finding of the zero-order solution is based on the CR isotropy during propagation in the space, whereas the anisotropy is taken into account when finding the next amendments. To begin with, the method is applied to solve the problem of CR modulation where the diffusion coefficient $\kappa$ and the solar wind speed $u$ are constants with an Local Interstellar Spectra (LIS) spectrum. The solution obtained with two iterations was compared with an analytical solution and with numerical solutions. Finally, solutions that have only one iteration for two problems of CR modulation with $u = constant$ and the same form of LIS spectrum were obtained and tested against numerical solutions. For the first problem, $\kappa$ is proportional to the momentum of the particle $p$, so it has the form $\kappa=k_0\eta$, where $\eta=\frac{p}{m_0c}$. For the second problem, the diffusion coefficient is given in the form $\kappa=k_0\beta\eta$, where $\beta=\frac{v}{c}$ is the particle speed relative to the speed of light. There was a good matching of the obtained solutions with the numerical solutions as well as with the analytical solution for the problem where $\kappa = constant$.

Figures

Tables

Table 1. Table 1: The estimates of the errors (in per cent) of the force-field solutions and AI solutions in comparison with HPF and with the numerical solutions at 1 au for u = 4 × 10 5 𝑢 4 superscript 10 5 u=4\times 10^{5}  m s -1 , α = 4 𝛼 4 \alpha=4 and for various x 0 subscript 𝑥 0 x_{0} (for more details, see the the text).

	$\kappa =$ constant					$\kappa =k_0\eta$					$\kappa =k_0\beta \eta$
	$\frac{N_0}{\mathrm{HPF}}$	$\frac{N_0+N_1}{\mathrm{HPF}}$	$\frac{N_0+N_1+N_2}{\mathrm{HPF}}$	$\frac{\mathrm{HPF}}{\mathrm{CN}}$		$\frac{N_0}{\mathrm{CN}}$	$\frac{N_0+N_1}{\mathrm{CN}}$	$\frac{N_0}{B-p}$	$\frac{N_0+N_1}{B-p}$		$\frac{N_0}{\mathrm{CN}}$	$\frac{N_0+N_1}{\mathrm{CN}}$	$\frac{N_0}{B-p}$	$\frac{N_0+N_1}{B-p}$
$x_0=1.2$	7.7	-0.001	-0.009	-0.17		10.7	1.13	12.16	2.66		6.58	-0.1	7.77	1.17
$x_0=3$	31	3.9	0.08	-0.04		9.25	0.52	6.54	-1.55		-2.71	-2.77	0.31	0.26
$x_0=5$	55.3	17.1	2.9	0.71		-2.78	-3.03	3.24	3.24		-12.1	-5	-8.05	-1.13
$x_0=6.5$	69	30.6	8.8	0.96		-8.77	-4.63	-7.36	-0.74		-18.2	-5.92	-21.12	-8.47
$x_0=9$	83.9	53.6	26.1	1.94		-21.8	-6.82	-15.01	2.21		-26.7	-6.57	-43.38	-20.51

Table 2. Table 2: The estimates of the errors (in per cent) of the force-field solutions and AI solutions in comparison with the numerical solutions for problems where κ = k 0 ​ η 𝜅 subscript 𝑘 0 𝜂 \kappa=k_{0}\eta and κ = k 0 ​ β ​ η 𝜅 subscript 𝑘 0 𝛽 𝜂 \kappa=k_{0}\beta\eta with x 0 = 1.2 subscript 𝑥 0 1.2 x_{0}=1.2 ( r = 1 𝑟 1 r=1 au, u = 4 × 10 5 𝑢 4 superscript 10 5 u=4\times 10^{5}  m s -1 , k 0 = 5 × 10 18 subscript 𝑘 0 5 superscript 10 18 k_{0}=5\times 10^{18}  m 2  s -1 ), α = 4 𝛼 4 \alpha=4 and for energies higher than η = 0.2 𝜂 0.2 \eta=0.2 ( T = 19 𝑇 19 T=19 MeV) (for more details, see the the text).

	$\kappa =k_0\eta$					$\kappa =k_0\beta \eta$
	$\frac{N_0}{\mathrm{CN}}$	$\frac{N_0+N_1}{\mathrm{CN}}$	$\frac{N_0}{B-p}$	$\frac{N_0+N_1}{B-p}$		$\frac{N_0}{\mathrm{CN}}$	$\frac{N_0+N_1}{\mathrm{CN}}$	$\frac{N_0}{B-p}$	$\frac{N_0+N_1}{B-p}$
$\eta =0.2$ ($T=19$ MeV)	7.36	0.45	8.14	2.08		-12.73	-5.38	-9.6	-2.45
$\eta =0.3$ ($T=41$ MeV)	11.4	1.55	13.23	3.7		-2.98	-3.13	-1.39	-1.54
$\eta =0.4$ ($T=72$ MeV)	11.53	1.44	13.18	3.23		3.22	-1.23	4.27	-0.12
$\eta =0.5$ ($T=111$ MeV)	10.7	1.13	12.16	2.66		6.58	-0.1	7.77	1.17
$\eta =0.6$ ($T=156$ MeV)	9.67	0.84	11.15	2.37		8.06	0.37	9.01	1.39
$\eta =0.7$ ($T=207$ MeV)	8.67	0.61	9.92	1.89		8.45	0.47	9.53	1.65
$\eta =0.8$ ($T=263$ MeV)	7.76	0.44	8.69	1.38		8.25	0.41	9.13	1.37
$\eta =0.9$ ($T=324$ MeV)	6.97	0.32	7.73	1.07		7.77	0.29	8.51	1.09
$\eta =1$ ($T=389$ MeV)	6.28	0.23	6.89	0.84		7.18	0.17	7.82	0.86
$\eta =2$ ($T=1.16$ GeV)	2.68	-0.0008	2.9	0.22		2.94	-0.18	3.31	0.21
$\eta =5$ ($T=3.85$ GeV)	0.64	0.003	0.63	-0.005		0.59	-0.08	0.71	0.04
$\eta =10$ ($T=8.49$ GeV)	0.19	0.006	0.22	0.04		0.15	-0.04	0.18	-0.01
$\eta =20$ ($T=17.85$ GeV)	0.05	0.005	0.08	0.03		0.03	-0.02	0.06	0.01

Equations

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\kappa \frac{\partial f}{\partial r}\right)-u\frac{\partial f}{\partial r}+\frac{1}{3r^2}\frac{\partial }{\partial r}(r^{2}u)\frac{\partial f}{\partial \ln p}=\frac{\partial f}{\partial t},$$

$$f(\text{r},p,t)=\frac{1}{4\pi }\int F(\text{r},\text{p},t)\mathrm{d}\Omega .$$

$$N(\text{r},p,t)=\int F(\text{r},\text{p},t)\mathrm{d}\Omega ,$$

$$\text{J}(\text{r},p,t)=\int F(\text{r},\text{p},t)\text{upupsilon}\mathrm{d}\Omega ,$$

$${\nabla }_r{j}_r+{\nabla }_p{j}_p=-\frac{\partial N(r,p,t)}{\partial t},$$

$$j_r=-\kappa (r,p,t)\frac{\partial N(r,p,t)}{\partial r}-\frac{up}{3}\frac{\partial N(r,p,t)}{\partial p}$$

$$j_p=\frac{up}{3}\frac{\partial N(r,p,t)}{\partial r}$$

$$N=N_0(r,p,t)+N_1(r,p,t)+N_2(r,p,t)+...$$

$$\{\begin{array}{l}{j}_{r,0}(r,p,t)=0,\\ N_0(r_0,p,t)=\mathrm{LIS}.\end{array}$$

$${\nabla }_r{j}_{r,i+1}+{\nabla }_p{j}_{p,i}=-\frac{\partial N_i(r,p,t)}{\partial t},$$

$$\frac{1}{r^2}$$

$$=-\frac{\partial N_i}{\partial t}.$$

$$\{\begin{array}{l}{j}_{r,i+1}(0,p,t)=0,\\ N_{i+1}(r_0,p,t)=0.\end{array}$$

$$\{\begin{array}{l}\frac{\partial N_0(x,\eta )}{\partial x}+\frac{\eta }{3}\frac{\partial N_0(x,\eta )}{\partial \eta }=0,\\ N_0(x_0,\eta )=C{\eta }^{-\alpha }.\end{array}$$

$$\{\begin{array}{l}{N}_0(x,\eta )=f(\eta \exp (-\frac{x}{3})),\\ N_0(x_0,\eta )=C{\eta }^{-\alpha }.\end{array}$$

$$N_{0}=C\eta^{-\alpha}\exp\bigg{[}\frac{\alpha}{3}(x-x_{0})\bigg{]}.$$

$$\frac{1}{x^{2}}\frac{\mathrm{d}}{\mathrm{d}x}x^{2}\Bigg{[}-\frac{\mathrm{d}\Phi_{1}}{\mathrm{d}x}+\frac{\alpha}{3}\Phi_{1}\Bigg{]}+\frac{\alpha(3-\alpha)}{9}\exp\bigg{[}\frac{\alpha}{3}(x-x_{0})\bigg{]}=0.$$

$$\frac{\mathrm{d}\Phi_{1}}{\mathrm{d}x}-\frac{\alpha}{3}\Phi_{1}=\frac{\alpha(3-\alpha)}{9}\exp\bigg{(}-\frac{\alpha}{3}x_{0}\bigg{)}\sum\limits_{n=0}^{\infty}\frac{(\frac{\alpha}{3})x^{n+1}}{n!(n+3)}.$$

$$\Phi_{1}=B\exp\bigg{(}\frac{\alpha x}{3}\bigg{)}+\exp\bigg{(}\frac{\alpha x}{3}\bigg{)}\int\limits_{0}^{x}\exp\bigg{(}-\frac{\alpha x}{3}\bigg{)}\frac{\alpha(3-\alpha)}{9}\exp\bigg{(}-\frac{\alpha}{3}x_{0}\bigg{)}\sum\limits_{n=0}^{\infty}\frac{(\frac{\alpha}{3})x^{n+1}}{n!(n+3)}\mathrm{d}x.$$

$$\Phi_{1}=B\exp\bigg{(}\frac{\alpha x}{3}\bigg{)}+\frac{\alpha(3-\alpha)}{9}\exp\bigg{(}\frac{\alpha}{3}(x-x_{0})\bigg{)}\sum\limits_{n=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{(-1)^{m}(\frac{\alpha}{3})^{n+m}x^{n+m+2}}{n!m!(n+3)(n+m+2)}.$$

$$B=-\frac{\alpha(3-\alpha)}{9}\exp\bigg{(}-\frac{\alpha}{3}x_{0}\bigg{)}\sum\limits_{n=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{(-1)^{m}(\frac{\alpha}{3})^{n+m}x_{0}^{n+m+2}}{n!m!(n+3)(n+m+2)}.$$

$$N_{1}=C\eta^{-\alpha}\exp\bigg{[}\frac{\alpha}{3}(x-x_{0})\bigg{]}\frac{\alpha(3-\alpha)}{9}\sum\limits_{n=0}^{\infty}\sum\limits_{m=0}^{\infty}\frac{(-1)^{m}(\frac{\alpha}{3})^{n+m}(x^{n+m+2}-x_{0}^{n+m+2})}{n!m!(n+m+2)(n+3)}.$$

$$N_2=$$

$$\displaystyle\times\Bigg{[}\frac{n+m+k+2}{(n+m+k+4)(n+m+k+t+3)}(x^{n+m+k+t+3}-x_{0}^{n+m+k+t+3})-\frac{k{x_{0}^{n+m+2}}}{(k+2)(k+t+1)}(x^{k+t+1}-x_{0}^{k+t+1})\Bigg{]}.$$

$$N_{\mathrm{es}}=C\eta^{-\alpha}\frac{F\bigg{(}\frac{2\alpha}{3},2;x\bigg{)}}{F\bigg{(}\frac{2\alpha}{3},2;x_{0}\bigg{)}},$$

$$N_{\mathrm{es},\alpha =3}=C{\eta }^{-3}\exp [x-x_0].$$

$$\{\begin{array}{l}\frac{\partial N_0(x,\eta )}{\partial x}+\frac{1}{3}\frac{\partial N_0(x,\eta )}{\partial \eta }=0,\\ N_0(x_0,\eta )=A{\eta }^{-\alpha }.\end{array}$$

$$\{\begin{array}{l}{N}_0(x,\eta )=f(\eta -\frac{x}{3}),\\ N_0(x_0,\eta )=A{\eta }^{-\alpha }.\end{array}$$

$$N_0=A{\left[\eta -\frac{x-x_0}{3}\right]}^{-\alpha }.$$

$$\frac{1}{x^2}\frac{\partial }{\partial x}{x}^2\left[-\eta \frac{\partial N_1}{\partial x}-\frac{\eta }{3}\frac{\partial N_1}{\partial \eta }\right]=-\frac{\alpha A}{9{\eta }^2}\frac{\partial }{\partial \eta }\left[{\eta }^3{\left(\eta -\frac{x-x_0}{3}\right)}^{-\alpha -1}\right].$$

$$-\frac{\partial N_1}{\partial x}-\frac{1}{3}\frac{\partial N_1}{\partial \eta }=I,$$

Full text

††pubyear: 2017

1 Introduction

The classic form of the equation of transport for cosmic rays (TPE) in a spherically symmetric interplanetary region, without taking into account particle drift and any local sources, is (Parker, 1965; Dolginov $\&$ Toptygin, 1967; Gleeson $\&$ Axford, 1967; Dolginov $\&$ Toptygin, 1968; Jokipii $\&$ Parker, 1970)

[TABLE]

where $f(\textbf{{r}},\>p,\>t)$ is the omnidirectional distribution function, r is the radial distance from the Sun, $p$ is the modulus of the momentum of the particle, $t$ is time, $u$ is the radial solar wind (SW) speed and $\kappa(r,p,t)$ is the effective radial diffusion coefficient, i.e. it is $\kappa_{rr}$, which is the one from the nine elements of the diffusion tensor K.

The omnidirectional distribution function, $f(\textbf{{r}},\>p,\>t)$, relates to the full cosmic-ray (CR) distribution function, $F(\textbf{{r}},\>\textbf{{p}},\>t)$, as

[TABLE]

On the other hand, if particles are scattered by random inhomogeneities so that the angular distribution of the directions of the particles becomes isotropic, then it is possible to apply the diffusion approximation (Morse $\&$ Feshbach, 1953), i.e. to use the moments of $F(\textbf{{r}},\>\textbf{{p}},\>t)$:

[TABLE]

[TABLE]

where $N(\textbf{{r}},\>p,\>t)$ is the phase density of the particles and $\textbf{{J}}(\textbf{{r}},\>p,\>t)$ is the vector of the particle flux in space. The integration in (2), (3) and (4) is over the solid angle element $\mathrm{d}\Omega$ in momentum space. Thus, taking into account (2) and (3), the TPE may be rewritten, as in Dolginov $\&$ Toptygin (1967, 1968), in a form involving $N(\textbf{{r}},p,t)$. Note that in the form was obtained by Dolginov $\&$ Toptygin (1967, 1968) from the kinetic Boltzman equation in which the regular as well as the random components of the interplanetary magnetic fields were taken into account. Such form can be represented as a continuity equation that describes the propagation of particles for a spherically symmetric event in phase space (Dorman et al., 1978; Gleeson $\&$ Webb, 1978):

[TABLE]

where $\nabla_{r}=\frac{1}{r^{2}}\frac{\partial}{\partial r}r^{2}$ is the radial component of the divergence operator and $\nabla_{p}=\frac{1}{p^{2}}\frac{\partial}{\partial p}p^{2}$ is the component of the divergence operator in momentum space. The value

[TABLE]

is the CR flux in space, and

[TABLE]

is the CR flux in momentum space. For stationary case, equation (5) describes CR propagation that is reduced to CR advection by the SW. Particles with different energies are transferred with different rates, which is described by term $-\frac{up}{3}\frac{\partial N(r,\>p,\>t)}{\partial p}$, and, on the other hand, CRs penetrate diffusively into the SW with the rate $\kappa(r,\>p,\>t)\frac{\partial N(r,\>p,\>t)}{\partial r}$. This is accompanied by a change in the particle energy, which is determined by the balance between the oncoming and following encounters of particles with the moving irregularities of the interplanetary magnetic field ($\frac{up}{3}\frac{\partial N(r,\>p,\>t)}{\partial r}$).

The possibility to obtain the exact analytical solutions (1) for various forms of $\kappa$ and $u$ [not included when (1) is greatly simplified for $u\propto\frac{1}{r^{2}}$] has always attracted the attention of researchers. In this way, significant results were obtained, but only for the steady-state case with $u=const$ and involuntary forms of $\kappa$, i.e. a form of $\kappa$ that does not depend simultaneously on $r$, $p$ or $R$ (the particle rigidity) and $\beta$, where $\beta=\frac{\upupsilon}{c}$ is the particle speed relative to the speed of light. In particular, analytic solutions have been derived by Dolginov $\&$ Toptygin (1967, 1968) for $\kappa=$ constant (see Appendix A, equation 23), by Fisk $\&$ Axford (1969) for $\kappa=k_{0}T^{\alpha}r^{b}$ ($b>1$, T is particle kinetic energy), by Webb $\&$ Gleeson (1977) for $k_{0}(p)r^{b}$ (via Green’s theorem and Green’s function) and by Shakhov $\&$ Kolesnik (2008) for $k_{0}r^{b}$ (via Mellin transform). Attempts to obtain the exact analytical solution (1) for arbitrary forms of $\kappa$ and $u$ and, especially, for unsteady cases encounter serious mathematical difficulties. Therefore, the development of the approximate analytical methods but with an acceptable accuracy for solving steady-state as well as time-dependent problems is an urgent problem. Furthermore, often it is necessary to have the qualitative characteristics of the processes and the ability to freely operate by the variables in analytical formulae. This is a big advantage compared to the numerical calculations, which, moreover, are almost always time-consuming.

2 About the method

In this paper, we continue with the development of an analytically iterative method for solving problems of CR modulation. The method was initially proposed by Shakhov $\&$ Kolesnyk (2006, 2009). The point of the method is as follows: Let us represent the particle density as the sum

[TABLE]

where $N_{0}(r,\>p,\>t)$ is the zero-order solution and $N_{1,2,3..}(r,\>p,\>t)$ are amendments to the zero solution. To obtain the zero-order solution, we assume that the anisotropy of the CRs during propagation in space is small, and that there is no radial CR flux at any point of space. This is the force-field approximation, which was initially introduced by Gleeson $\&$ Axford (1968) and furthered in (Gleeson $\&$ Urch, 1973), which presented a new development for its solutions. At the boundary of the CR modulation ($r_{0}$), for example, at the heliosphere boundary, the zero-order solution must satisfy a boundary condition involving its spectrum; for example be equal to some LIS spectrum. So, to find of the zero-order solution, the following conditions must be satisfied:

[TABLE]

The amendments to the zero-order solution can be obtained from (5) by the following recurrence relation:

[TABLE]

where $i$ is an iteration index, which changes from the zero to infinity.

Taking into account (6) and (7), the last relation can be written in the following form:

[TABLE]

Here, in order to obtain $N_{i+1}$, we will be guided by the fact that any radial CR flux (such as from the zero solution and also from the amendments) in the centre of the heliosphere (near the Sun) must be absent, i.e. $j_{r,i}(0,\>p,\>t)=0$. This follows from the spherical symmetry of the problem. It should be noted that the amendments of the fluxes will be small but not equal to zero. This means that as more amendments are obtained, more of the CR anisotropy in the space will be taken into account. And the last condition is for finding $N_{i+1}$. As a common solution is presented in the form (8) and the zero solution $N_{0}$ satisfies the boundary condition at $r_{0}$, this means that all amendments $N_{i+1}$ at $r_{0}$ must be equal to zero.

To sum up the above, for the finding of the amendments necessary to fulfil such conditions,

[TABLE]

Thus, in the first step, using (6) and (9), we may obtain the zero solution [$N_{0}(r,\>p,\>t)$]. Then, inserting it into (2) and applying (12), we obtain the first amendment [$N_{1}(r,p,t)$]. For obtaining the next amendment, it is necessary to substitute $N_{1}(r,\>p,\>t)$ into (2) and again apply (12), and so on. A common solution, as mentioned above, will be expressed in the form (8).

In the following sections of this paper, we test the analytically iterative method (AI) for a constant diffusion coefficient and two shapes of non-constant diffusion coefficient.

We used the stochastic integration B-p method and the Crank-Nicolson method (CN) to verify the solutions obtained by AI. B-p is a stochastic integration method to solve the TPE backwards in time. To find stochastic solutions of the Fokker-Planck equation (FPE), Itos lemma was used to rewrite the TPE into a set of stochastic differential equations. This method historically is treated as a forward-in-time method of solution. The backward-in-time solution was introduced by Kota (1977) and was described in detail in, for example, Zhang (1999). We apply the 1D backward stochastic integration method that is used, for example, by Yamada, Yanagita $\&$ Yoshida (1998). A comparison of the backward and forward methods with a systematic error estimation for realistic LIS was published in Bobik et al. (2016), where the B-p method is described in detail.

CN is a numerical unconditionally stable finite-difference method to solve FPEs (Diaz, 1958). The method was proposed by Fisk (1971) to solve TPEs, and then used by many other authors. Here, the CN method is applied in the form given in Batalha (2012).

To compare the AI method, we use two completely independent numerical methods, to avoid any numerical artefacts or errors in the verification of the method.

3 A constant diffusion coefficient

To verify the AI approach, we start with the simpler case where the diffusion coefficient $\kappa$ is constant. Obtaining a solution by the AI method for $\kappa=$ constant and an LIS at the heliosphere boundary that has the form $Cp^{-\alpha}$ is discussed in Appendix A. A comparison with an analytical solution (see equation 23) that is based on the confluent hypergeometric functions (HPF method) is available for this case, which therefore allows us to very accurately estimate the precision of the AI method. If not mentioned otherwise, all results published in this paper were evaluated for protons in a spherically symmetric heliosphere with radius 100 au, where the SW has a speed of 400 km s*-1*.

In Fig. 1, the phase density function $N$ at $1$ au is shown for $\kappa=5\times 10^{18}$ m2 s*-1*  as are the slopes of the initial spectrum at the heliosphere boundary $\alpha=4$ evaluated by AI (red diamonds) and by HPF (black line), together with $N$ at the modulation boundary (blue line). Both solutions are very similar at $1$ au. The difference between them is of the order of 0.01 per cent at all energies from $\eta=0.01$ ($\eta=\frac{p}{m_{0}c}$, a proton with kinetic energy 47 keV) to $\eta=100$ (a proton with kinetic energy 92.89 GeV). As can be seen from Fig. 1, the modulation is relatively the same for all energies because the diffusion coefficient does not depend on the energy. More precisely, by saying modulation, we mean the ratio between the phase density inside the heliosphere (here at 1 au) and that normalised to the phase density for the same momentum at the modulation boundary (at 100 au). From Fig. 1, we can see that difference between the solution that was obtained by AI and that obtained by HPF at the selected $r$ does not depend on the energy.

We present in Fig. 2 a comparison of the phase density spatial distribution across the heliosphere for selected energies, as evaluated by the AI, HPF, and CN methods.

Presented are the ratios of the phase-density functions evaluated for $\kappa=5\times 10^{18}$ m2 s*-1*  and the slopes of the initial spectra at the heliosphere boundary $\alpha=3$ (left-hand panel) and $\alpha=4$ (right-hand panel). As for the case of constant $\kappa$, the ratio of AI/HPF does not depend on the energy (as can be seen from the AI solution and the exact solution 23); thus, the ratios as functions of the radius for $\eta=1$ (388 MeV) and $\eta=10$ (8.49 GeV) are very similar. Therefore, the green line in Fig. 2 representing the ratio for $\eta=10$ hides the red line of the $\eta=1$ solution in both figure panels.

For the special case when $\alpha=3$ with the mentioned value of $\kappa$, the AI solution is very precise, with a difference between the AI and HPF solutions at a level of $10^{-6}$ per cent. This is due to the fact that for this case, the first amendment $N_{1}$ and the second amendment $N_{2}$ (as well as all higher amendments that will be found) are equal to zero (see Appendix A). As a result, a common solution via AI is described by $N_{0}$ alone (equation 15) and is $N_{0}=C\eta^{-3}\exp[x-x_{0}]$. Moreover, HPF describes equation (24) (see Appendix A) for this case. Therefore, HPF and AI yield identical solutions for $\alpha=3$. The range of numerical error during the evaluation of the solutions is then very small. In Fig. 2, we also show the ratios of the phase-density functions evaluated by CN and HPF. CN, in both cases, attains a precision of $0.01$ per cent (precision as the difference from the analytic HPF solution). You can see that the AI solution for $\alpha=4$ is more precise than the CN solution. We evaluated the solutions for the same set of parameters also by the stochastic B-p method. The results from the B-p method confirmed the results presented by the AI and HPF methods.

The precision of the AI solution depends on $\kappa$ and slightly less on $\alpha$. In Fig. 3, we present the dependence of the ratio between the AI and HPF solutions on $\kappa$ at different distances $r$. For $\alpha=2$ and 4, the AI solutions are precise, of the order of $0.1$ per cent precision level, for diffusion coefficients $\kappa$ greater than $2\times 10^{18}$ m2 s*-1*. As could be expected, for $\alpha=3$, we see a very high precision at every position inside the heliosphere for $\kappa>3\times 10^{17}$ m2 s*-1*. For $\kappa=2\times 10^{18}$ m2 s*-1*, we evaluated the ratio of the AI and HPF solutions in dependence on $\alpha$. The results are presented in Fig. 4. The precision decreases as $\alpha$ becomes more different from $3$ and with decreasing $r$.

It should be noted that the precision of the AI solution may be improved by finding more amendments (as mentioned in the description of the method) and by increasing the order of the summations in the evaluation of the amendments. This is true in the presented case for the first and the second amendments. The latter is understandable, due to the fact that we have applied a Taylor series expansion for the exponential to obtain these amendments. But this has a limit due to limitations of the data types of computing languages. In particular, if a higher value of $\frac{ur}{\kappa}$ is used, then the number of terms in the sums in $N_{1}$ and $N_{2}$ must be increased in order to attain an acceptable accuracy. As a result, this situation significantly increases the time for the evaluation of the amendments.

4 Non-constant diffusion coefficient

4.1 Solution for $\kappa=k_{0}\eta$

The solution for the case of a non-constant diffusion coefficient with shape ${\kappa=k_{0}\eta}$ is derived in Appendix B on the basis of the AI method. Here, it should be noted that although $N_{1}$ (see equation B) is expressed through integrals, but to construct the corresponding curve, it is sufficient to use one of the mathematical packages such as MAPLE, MATHEMATICA, MATLAB, etc. In this case, the result will be obtained faster than the numerical calculation of such a problem. At the same time, numerical calculations (for CN as well as the B-p stochastic method) of the problem should use or apply the program code to obtain the required points. And then based on these points and by means of the software, we can build the appropriate curve.

An analytical solution for this shape of the diffusion coefficient and power-law unmodulated spectrum is not available. In particular, it should be noted that for a similar problem where for $\kappa=k_{0}p^{\alpha}$, an analytical solution was obtained by Cowsik $\&$ Lee (1977) only for the constant SW speed and for a specific form of LIS spectrum, but not for an arbitrary form. Herewith, the form of LIS spectrum that was examined is not a power-law kind.

We compared the AI solutions, i.e. $N_{0}$ and $N_{0}+N_{1}$, with the solutions provided by CN and the B-p stochastic method. By using two independent numerical methods (B-p and CN), we want to ensure the validity of the comparison with the AI method.

In Fig. 5, a comparison of the solutions for energies starting from $\eta=0.15$ ($T=11$ MeV) at 1 au for $\alpha=4.0$ and $k_{0}=5\times 10^{18}$ m2 s*-1*is presented. The zero-order solution $N_{0}$ is indicated in the figure by AI$(N_{0})$, and the sum of the zero-order solution and the first amendment $N_{0}+N_{1}$ is indicated by AI$(N_{0}+N_{1})$. Over a wide range of energies, all four solutions, i.e. $N_{0}$, $N_{0}+N_{1}$, CN and B-p, in Fig. 5 are similar. As can be seen from the right-hand panel of Fig. 5, comparisons of the AI solution with that of the CN method, and of the AI solution with B-p, show that the difference is less than 1 per cent for energies higher than 324 MeV ($\eta>0.9$). The CN and AI solutions are more similar than the AI and B-p solutions. The ratio AI/CN remains within the 1 per cent difference range for energies higher than 133 MeV ($\eta>0.55$). The CN solution and the B-p solution differ by less than 1 per cent over a wide range of energies: from 263 MeV ($\eta=0.8$) to higher energies.

To look more closely at the precision of the AI solution, we compare the solutions provided by the AI and CN methods at 1 au for a range of values of the diffusion coefficient and for the selected slopes of the power-law unmodulated spectrum at 100 au. In Fig. 6, the ratios between the AI and CN solutions at 1 au for different ${\kappa_{0}}$ and ${\alpha=3.0}$, $4.0$ and $5.0$ and selected energies are shown. The dotted lines represent the ratio of the solutions AI$(N_{0})$/CN and the solid lines represent the ratio AI$(N_{0}+N_{1})$/CN. Here, AI$(N_{0})$ is the zero-order solution and AI$(N_{0}+N_{1})$ is the sum of the zero-order solution and the first amendment. The ratios become closer to unity with increasing values of $\kappa_{0}$. We can see from the figure that when $k_{0}$ is larger than $2\times 10^{17}$ m2 s*-1*, then all AI$(N_{0}+N_{1})$ differ from the CN solutions by less than 1 per cent for energies over $\eta=20$ (17.85 GeV) and in the range of $\alpha$ tested. Also, for $k_{0}$ larger than $2\times 10^{18}$ m2 s*-1*, all AI solutions differ from CN solutions by less than 1 per cent for energies over $\eta=2$ (1.16 GeV).

4.2 Solution for $\kappa=k_{0}\beta\eta$

The solution for a non-constant diffusion coefficient with shape ${\kappa=k_{0}\beta\eta}$ and for a power-law unmodulated spectrum is derived in Appendix C. The comparison of the AI solutions with the solutions that were obtained by the CN and B-p methods for energies starting from $\eta=0.1$ ($T=5$ MeV) at 1 au for $k_{0}=5\times 10^{18}$ m2 s*-1*and $\alpha=4.0$ is presented in Fig. 7.

The ratios of the AI solution to CN and B-p solutions for the zero-order solution $N_{0}$ [indicated by AI$(N_{0})$] and for the common solution that consists of the zero solution and the first amendment $N_{0}+N_{1}$ [indicated by AI$(N_{0}+N_{1})$] are shown in the figure’s right-hand panel.

From the inspection of the figure, we can see that the common AI solution differs by less than 1 per cent from the CN solution for energies bigger than 79 MeV ($\eta>0.42$). The B-p solution stays within the 1 per cent difference region from energies $\geq$ 349 MeV ($\eta>0.94$), but for smaller energies, it oscillates around the 2 per cent difference region.

The precision of the AI solution for different values of $k_{0}$ for ${\kappa=k_{0}\beta\eta}$ was similarly reviewed, as in the previous two diffusion cases. As one can see from Fig. 8, the ratio AI/CN goes to unity with increasing values of the two factors, namely with $k_{0}$ and the energy. But with increasing $\alpha$, the situation is slightly the opposite and ambiguous. For example, there is about a 3 per cent difference between the AI and CN solutions when $\alpha=3$ for an energy of 111 MeV ($\eta=0.5$) and $k_{0}=1\times 10^{18}$ m2 s*-1*, and the difference is already about 5.5 per cent when $\alpha=4$ for the same energy and $k_{0}$. But when $\alpha=5$, then the difference again becomes 3 per cent for the same parameters. The situation is similar for the other energies.

For example, when the energy is 1.16 GeV ($\eta=2$) and $k_{0}=4\times 10^{17}$ m2 s*-1*, the deferences are 2 per cent for $\alpha=3$, 1 per cent for $\alpha=4$ and 6.5 per cent for $\alpha=5$. Here, it should be noted that if the modulation parameter $x_{0}$ ($x_{0}=\frac{ur_{0}}{k_{0}}$) is less than 10, i.e. when $k_{0}>6\times 10^{17}$ m2 s*-1*, then the deviation of the AI solution from the CN solution reaches no more than 6 per cent for any checked $\alpha$ and energy. We think this is a good result taking into account that numerical methods also have accuracy limits, as can be seen from the ratio CN/B-p (Fig. 7).

Finally, the dependence of the precision of the solution on ${\alpha}$ at 1 au for two selected ${k_{0}}$ is presented in Fig. 9. The values used for ${k_{0}}$ were chosen as the border values of the modulation range, between a very strong maximum and a weak minimum of the solar cycle. $k_{0}$ was evaluated from the modulation strength $\phi$ (Usoskin et al., 2005; Usoskin, Bazilevskaya $\&$ Kovaltsov, 2011). The weakest modulation during solar cycles 22 and 23 occurred in 1987 February (solar cycle 22) with $\phi=311$ MV. From the average SW speed in 1987 February, $u=5.08\times 10^{5}$ m s*-1*(OMNIweb, 2017), we found that $k_{0}=7.56\times 10^{18}$ m2 s*-1*. The strongest modulation during solar cycles 22 and 23 occurred in 1991 June with $\phi=2016$ MV, $u=5.23\times 10^{5}$ m s*-1*, which leads to $k_{0}=1.2\times 10^{18}$ m2 s*-1*. Using these two values of $k_{0}$, we can show the precision of the evaluated methods for a range of values appearing in the last solar cycles. During very strong modulations in solar maxima, when $k_{0}$ reaches the smallest values, the AI/CN ratios stay within the 1 per cent difference range for energies over 3.85 GeV ($\eta>5$) for all tested values of $\alpha$ (see the left-hand panel in Fig. 9). The precision is one order of magnitude better during solar minima, when $k_{0}$ reaches its highest values. As we can see in the right-hand panel of Fig. 9, the precision is better than 0.1 per cent for energies higher than 1.16 GeV ($\eta>2$) for all tested $\alpha$. It should be noted that the difference is less than the 6 per cent between AI and CN for the solar minima period as well as for the solar maxima period for all tested energies and $\alpha$.

Let us note that we published codes of models used in this paper on the web page www.CRmodels.org.

5 Discussion and Conclusions

In this paper, we have obtained solutions for three steady-state problems of CR modulation by means of an AI method. By summing all three cases, we can reach the following conclusion: The accuracy of these solutions at a given point $r$ is dependent more on the value of the modulation parameter $x_{0}$ ($x_{0}=\frac{ur_{0}}{\kappa}$ for the problem, where $\kappa=$ constant and $x_{0}=\frac{ur_{0}}{k_{0}}$ for problems, where $\kappa=k_{0}\eta$ or $\kappa=k_{0}\beta\eta$) than on the choice of the slope of the initial spectrum $\alpha$.

Thus, for $x_{0}<10$, the obtained solutions have minimal differences in accuracy in comparison with numerical solutions and the exact solution for the problem when $\kappa=$ constant. In particular, for the problem where the diffusion coefficient is a constant, so that $x_{0}=1.2$ (i.e. when $k_{0}=5\times 10^{18}$ m2 s*-1*, $u=4\times 10^{5}$ m s*-1*and $r_{0}=100$ au), the obtained AI solution is even closer to the exact solution than its numerical solution. For a task where the diffusion coefficient is proportional to the momentum, and for the same value of $x_{0}$, we have shown that the AI solution deviates from the numerical CN solution by not more than $1.55$ per cent over the entire considered range of energies and the deviation is less than $1$ per cent starting from energies of $T=133$ MeV ($\eta=0.55$). For the problem where $\kappa=k_{0}\beta\eta$, the deviation of the AI solution from the numerical CN solution is less than $1$ per cent beginning at an energy of $T=79$ MeV ($\eta=0.42$), which is better than for the previous problem. It should be noted that although the comparison of the AI solution with the numerical solution B-p is more advantageous for this case (the B-p solution has a deviation of not more than $2.5$ per cent over all of the investigated range of energies) and less profitable for $\kappa=k_{0}\eta$, here we have provided only one of them because the difference between the two numerical methods can reach $10$ per cent, as seen from Fig. 7. It should also be noted that the situation with the obtained AI solutions at smaller values of $x_{0}$, i.e. $x_{0}<1$, is even better than the above-examined case.

In this paper, we have evaluated the possible $x_{0}$ values under real physical conditions; as an example, we have chosen different cycles of solar activity. The maximum $x_{0}$ value that we have found during maximum solar activity does not exceed $x_{0}=6.5$ (i.e. when $k_{0}=1.2\times 10^{18}$ m2 s*-1*, $u=5.23\times 10^{5}$ m s*-1*and $r_{0}=100$ au). In very rare situations, when there is a very fast SW and a low value of $k_{0}$, then $x_{0}$ may reach a value over 10.

In this case, the difference between the AI and numerical solutions becomes larger but, at the same time, gives us much better accuracy than the force-field solutions ($N_{0}$). This result is shown in Tables 1 and 2, along with estimates of the errors (in per cent) of $N_{0}$ and AI solutions in comparison with HPF (only for the problem where $\kappa=$ constant, there is an exact solution) and numerical solutions at 1 au for $u=4\times 10^{5}$ m s*-1*and $\alpha=4$ for various cases. In this table, the case $\kappa=$ constant corresponds to the exact HPF solution only. Therefore, Table 1 demonstrates estimates for all three problems mentioned above. The parameters were chosen as follows: $x_{0}=1.2,\>3,\>5,\>6.5,\>9$; energy $\eta=0.5$ ($\kappa=k_{0}\eta$ and $\kappa=k_{0}\beta\eta$). Table 2 displays estimates for the cases when $\kappa$ is not constant (for $x_{0}=1.2$) and for various energies that are higher than $\eta=0.2$. These tables show the importance of the inclusion of the amendment $N_{1}$ in comparison with the force-field solution for the selected $x_{0}$ (Table 1) and for the selected $\eta$ (Table 2). The result is even better when the second amendment is also taken into account (for the problem with $\kappa=$ constant). We can see that accuracy deteriorates with the increasing $x_{0}$ (for all the three problems) and with the decreasing $\eta$ (for the problems where $\kappa$ is not constant).

Therefore, we think that using more amendments in this method would produce more accurate solutions for those rare occasions of $x_{0}$, when it is necessary. To show the above and also that the common iterative solution (8) converges with the exact solution, a problem was examined in Appendix E, where LIS has the time dependence. For this purpose, a term $\frac{\partial N_{i}}{\partial t}$ in (2) was taken into account. In addition, it has been shown that the exact solution of the considered problem is unknown because of the difficulty to carry out the inverse Laplace transform for equation (49); therefore, the AI method in this situation, as for another problems, has a big advantage.

So we believe that the obtained solutions for the considered problems by means of the proposed method can be regarded as analytical solutions of such problems with some remarks. The method can be widely applied for solving the various steady-state and non-stationary problems of CR modulation with almost arbitrary dependences of the SW speed, forms of LIS spectrum and the diffusion coefficient ($\kappa$), which is simultaneously dependent on $p$ (or $R$), $\beta$ and $r$. Depending on $u$ and/or $\kappa$, the proposed method may have some restrictions. This is why we have used above the word ‘almost’. To find at least the zero-order solution of (6), we should solve the transcendental equation, which at the same time, depends on $u$, $\kappa$. Usually, it is difficult to solve analytically or even not possible. For example, if $u=$ constant and we choose a form $k$ such as $\kappa=k_{0}(1+\eta+\eta^{2})$, when we insert them into (6), we will obtain an equation that may be solved only by a numerical way. Here, it should be noted that ‘1’ in the form of $\kappa$ ‘prevents’ solving this example analytically. We guess that other the restrictions on the use of the method do not exist.

We think that by means of the proposed approach, it is possible to theoretically examine a stationary model of the spherically symmetric heliosphere to describe the physical processes during CR propagation therein (Kolesnik $\&$ Shakhov, 2012; Kota $\&$ Jokipii, 2014; Fedorov, 2015). The model includes an environment that consists of two adjacent spherically symmetric (relative to the Sun) regions. The internal part of the heliosphere is limited by a termination shock (TS), which is placed at $r_{0}$ from the Sun, the radial SW speed is constant and supersonic, $V_{in}$, but $\kappa$ has a form: $\kappa=k_{1}\beta\eta\frac{r}{r_{0}}$. Beyond TS, the SW speed jumps to a value $V_{in}/3$ (Richardson $\&$ Burlaga, 2013) and spreads with deceleration in the outer part of the heliosphere (the heliosheath), which is limited by the heliopause (HP). Such deceleration can be presented in the form of a power law depending on the heliocentric distance with index ‘n’, i.e. $V_{ihs}=\frac{V_{in}}{3}{(\frac{r}{r_{0}})}^{-n}$; here $\kappa$ has the form $\kappa=k_{2}\beta\eta$, where $k_{1}$ and $k_{2}$ may be determined from data obtained by Voyager 1. Due to the absence of the SW speed inside the interstellar medium (beyond HP), the phase density of the particles will be determined by only the diffusion process under the constant diffusion coefficient. Note, that Kota $\&$ Jokipii (2014) and Fedorov (2015) considered the easier mathematical problem, when $n=2$ ($\mathbf{\nabla}\cdot\textbf{{V}}=0$, i.e. beyond the TS, no adiabatic energy changes occur for CRs). It should be noted here that although $\mathbf{\nabla}\cdot\textbf{{V}}$ for large $n$ in the heliosheth is a very small quantity, it is not equal to zero. In this case, the third term in (1) will be proportional to $\mathbf{\nabla}\cdot\textbf{{V}}\cdot\eta\frac{\partial N}{\partial\eta}$. This means that neglecting the whole term is not correct, as its contribution may be significant for a range of low-energy particles ($\eta<1$). In a numerical way, this effect was shown by Langner et al. (2006), where the symmetrical 2D model for different scenarios, $V_{ihs}$, in the heliosheath was studied. In particular, it was shown there that changing the profiles for $V_{ihs}$ has a noteworthy effect on barrier modulation at low energies and that the barrier at such energies in the heliosheath is wider and deeper for the $V_{ihs}\propto\frac{1}{r^{8}}$ scenario than for $V_{ihs}\propto\frac{1}{r^{2}}$. In this way, using this method, we think that it is possible to derive a solution of the problem for given diffusion coefficients as well as with an arbitrary $n$ to look for possible effects on theory.

The method may be applied for the theoretical description of time-dependent CR modulation too, i.e. when we can observe solar-cycle-related changes in CR intensities due to changes in the modulation environment (Manuel et al., 2011). For this purpose, a time-dependent term needs to be taken into account in (2), and changes in the modulation environment will be described by the diffusion coefficient that varies with time.

We also think that the proposed method is a significant step towards developing an analytical approach to solve 2D problems of CR modulation (dependence on the radial distance and polar angle), taking into account particle drift.

Acknowledgements

We are grateful to the referee for his/her constructive comments.

This work was supported in part by the National Scholarship Programme of the Slovak Republic for the Support of the Mobility of Students, PhD Students, University Teachers, Researchers and Artists (YLK).

This work was partially supported by the Slovak Academy of Sciences, grant no. MVTS JEM-EUSO.

Appendix A $\kappa$ = constant

To obtain a solution of the problem, we introduce the following dimensionless variables: $x=\frac{ur}{\kappa},x_{0}=\frac{ur_{0}}{\kappa}$ and $\eta=\frac{p}{m_{0}c}$. Then, taking into account (6) and (9), we obtain a system of equations for finding the zero-order solution:

[TABLE]

The last system of equations can be rewritten in the following form:

[TABLE]

Then, the solution of the last system has the form

[TABLE]

To obtain the first and next amendments, we use a form such as $N_{i}=C\eta^{-\alpha}\Phi_{i}(x)$, and when we insert $N_{0}$ into (2), we will obtain an equation for $N_{1}$:

[TABLE]

If we expand the exponent in the last equation in a Taylor series, we obtain

[TABLE]

Here, the first condition (12) was applied, i.e. the integration was carried out from 0 to $x$. We will seek a solution of the last equation in the following form:

[TABLE]

After an expansion of the exponent in a series followed by integration, we obtain

[TABLE]

The last condition (12) for finding $B$ was used, i.e. $\Phi_{1}(x_{0})=0$, and then $B$ has the following form:

[TABLE]

And, as a result, the first amendment is

[TABLE]

When we perform similar actions in order to find the second amendment, we will obtain

[TABLE]

An exact solution for this problem was initially obtained by Dolginov $\&$ Toptygin (1967). In our notation, it has the following form:

[TABLE]

where $F(a,b;x)$ is the confluent hypergeometrical function (Bateman $\&$ Erdelyi, 1953). In particular, from the property of the confluent hypergeometrical function, it follows that for the simpler case when $\alpha=3$, the exact solution (23) has the form

[TABLE]

Appendix B $\kappa=k_{0}\eta$

To obtain a solution for the problem, we introduce the following dimensionless variables: $x=\frac{ur}{k_{0}},x_{0}=\frac{ur_{0}}{k_{0}}$ and $\eta=\frac{p}{m_{0}c}$. Then, taking into account (6) and (9), we obtain a system of questions for finding the zero-order solution:

[TABLE]

The last system can be rewritten in the following form:

[TABLE]

Then, the solution of the last system has the form

[TABLE]

After inserting $N_{0}$ into (2), we will obtain an equation for $N_{1}$:

[TABLE]

After differentiating the left-hand side of the last equation, we obtain

[TABLE]

where

[TABLE]

Here, the limits of integration were selected as being from 0 to $x$ in consequence of applying the first condition (12). Now, we consider the integral $I$ separately. If we define $y=\eta-\frac{\xi-x_{0}}{3}$, then the expression for $I$ can be rewritten in the following form:

[TABLE]

Whereas $N_{0}$ has the form (27), then, taking into account the second condition (12), it is necessary to search for an $N_{1}$ of the following form:

[TABLE]

And, as a result, the first amendment takes the following form:

[TABLE]

It should be noted that $N_{1}$ is not defined for $\alpha=0,1,2$.

Appendix C $\kappa=k_{0}\beta\eta$

To obtain a solution for the problem, we introduce the following dimensionless variables: $x=\frac{ur}{k_{0}},x_{0}=\frac{ur_{0}}{k_{0}}$ and $\eta=\frac{p}{m_{0}c}$. Note that for such a problem, the diffusion coefficient $\kappa$ can be rewritten in the following form: $\kappa=k_{0}\frac{\eta^{2}}{\sqrt{\eta^{2}+1}}$. This follows from the fact that $p=\frac{\upupsilon E}{c^{2}}$ and $E=\sqrt{p^{2}c^{2}+m_{0}^{2}c^{4}}$. Taking into account (6) and (9), we obtain a system of equations for finding the zero-order solution:

[TABLE]

The last system can be rewritten in the form

[TABLE]

Then, the solution of the last system has the form

[TABLE]

After inserting $N_{0}$ into (2), we will obtain an equation for $N_{1}$:

[TABLE]

After introducing the variable: $\theta=\sqrt{\eta^{2}+1}$ and simplifying the last equation, we obtain:

[TABLE]

where

$K=\int\limits_{0}^{x}\left(\left[\theta-\frac{x-x_{0}}{3}\right]^{2}-1\right)^{-\frac{\alpha}{2}-1}x^{2}\mathrm{d}x;$

$L=\int\limits_{0}^{x}\left(\left[\theta-\frac{x-x_{0}}{3}\right]^{2}-1\right)^{-\frac{\alpha}{2}-1}x^{3}\mathrm{d}x;$

$M=\int\limits_{0}^{x}\left(\left[\theta-\frac{x-x_{0}}{3}\right]^{2}-1\right)^{-\frac{\alpha}{2}-2}\left[\theta-\frac{x-x_{0}}{3}\right]^{2}x^{2}\mathrm{d}x$.

Note that the limits of integration in $K$, $L$, $M$ were selected as being from 0 to $x$ in consequence of applying of the first condition (12). The detailed calculation of the integrals was explained in Appendix D.

Due to the fact that the solution of the homogeneous equation (38) has the form $f(\theta-\frac{x}{3})$, and taking into account the second condition (12), it is necessary to search for an $N_{1}$ in the following form:

[TABLE]

By taking into account the form for $N_{1}$ (39) and also using (38), we obtain

[TABLE]

Appendix D calculation of the integrals K, L, M

For the calculation of the integrals $K$, $L$, $M$ from Appendix C, let us first consider the integral $I$ that has the following form:

[TABLE]

If we represent $(z^{2}-1)^{\beta}$ as a sum, i.e. $(z^{2}-1)^{\beta}=\sum\limits_{n=0}^{\infty}(-1)^{n}C_{\beta}^{n}\cdot(z^{2})^{\beta-n}$, where $C_{\beta}^{n}$ is the binomial coefficient (which can be represented through the Gamma function $\Gamma$ as $C_{\beta}^{n}=\frac{\Gamma(\beta+1)}{\Gamma(\beta-n+1)n!}$), we can apply here properties of $\Gamma$ such as $\frac{\Gamma(\beta+1)}{\Gamma(\beta-n+1)}=\frac{\Gamma(n-\beta)}{\Gamma(-\beta)}(-1)^{n}$ and see that $I$ assumes the following form:

[TABLE]

If we now use the fact that $\Gamma(z)=z\Gamma(z+1)$ and multiply (42) by $\frac{\Gamma(-\frac{k}{2}+\frac{1}{2}-\beta)}{\Gamma(-\frac{k}{2}-\frac{1}{2}-\beta)}$, then we obtain

[TABLE]

For the last expression, we can apply the definition of hypergeometrical function $F(a,b;c;z)$ (Bateman $\&$ Erdelyi, 1953): $F(a,b;c;z)=\sum\limits_{n=0}^{\infty}\frac{\Gamma(a+n)}{\Gamma(a)}\frac{\Gamma(b+n)}{\Gamma(b)}\frac{\Gamma(c)}{\Gamma(c+n)}\frac{z^{n}}{n!}$. Finally, we obtain

[TABLE]

And now, to use (44) to find the integrals $K$, $L$, $M$, we substitute $\theta-\frac{x-x_{0}}{3}=z$ in these integrals. As a result, these integrals take the following form:

$K=27\int\limits_{\theta-\frac{x-x_{0}}{3}}^{\theta+\frac{x_{0}}{3}}\left(z^{2}-1\right)^{-\frac{\alpha}{2}-1}(z^{2}-2\omega z+\omega^{2})\mathrm{d}z;$

$L=-81\int\limits_{\theta-\frac{x-x_{0}}{3}}^{\theta+\frac{x_{0}}{3}}\left(z^{2}-1\right)^{-\frac{\alpha}{2}-1}(z^{3}-3\omega z^{2}+3\omega^{2}z-\omega^{3})\mathrm{d}z;$

$M=27\int\limits_{\theta-\frac{x-x_{0}}{3}}^{\theta+\frac{x_{0}}{3}}\left(z^{2}-1\right)^{-\frac{\alpha}{2}-2}(z^{4}-2\omega z^{3}+\omega^{2}z^{2})\mathrm{d}z$,

where $\omega=\theta+\frac{x_{0}}{3}$. Eventually, if we apply (44) to each of these integrals, we obtain

[TABLE]

[TABLE]

[TABLE]

Appendix E $\kappa$ = constant and LIS has the time dependence

Let us consider a problem when the diffusion coefficient ($\kappa$) is constant inside the heliosphere, but LIS has a form $N_{0}(r_{0},\>\eta,\>t)=f(t)\eta^{-\alpha}$, i.e. LIS has dependence on the time and momentum of the particle on the boundary of the heliosphere $r_{0}$. For example, such dependence may be considered during a supernova explosion, when the density of the particles on the border of the heliosphere varies with time. To obtain solutions, i.e. as an exact solution as well as by means of the method, we introduce the following dimensionless variables: $x=\frac{ur}{\kappa},x_{0}=\frac{ur_{0}}{\kappa}$, $\eta=\frac{p}{m_{0}c}$ and $\tau=\frac{u^{2}t}{\kappa}$.

The exact solution of the problem may be obtained if to solve a system of equations:

[TABLE]

Here, the first equation is TPE in a form (1), taking into account the dimensionless variables, and the second equation is the boundary condition. Applying the Laplace transform, i.e. $N(x,\>\eta,\>s)=\int\limits_{0}^{\infty}\exp(-s\tau)N(x,\>\eta,\>\tau)\mathrm{d}\tau$ to the last system, we obtain

[TABLE]

Let us note that during integrating of the right-hand side of the last system, we used the fact that there was no phase density at the initial time, i.e. $N(x,\>\eta,\>0)=0$. The solution of the last system has the following form:

[TABLE]

If we expand the exponent and the confluent hypergeometrical function in the last equation in a Taylor series and restrict ourselves by the second-order $x$ and $x_{0}$ in such expansion, we obtain

[TABLE]

Now, divide the numerator of (E) ($1+\frac{\alpha}{3}x+(\frac{1}{6}s+\frac{1}{18}\alpha+\frac{1}{27}\alpha^{2})x^{2}+...$) by the denominator of (E) ($1+\frac{\alpha}{3}x_{0}+(\frac{1}{6}s+\frac{1}{18}\alpha+\frac{1}{27}\alpha^{2})x_{0}^{2}+...$) in a column

[TABLE]

And as a result, the analytical image of the problem retaining until the second order of $x$ and $x_{0}$ has the following form:

[TABLE]

To obtain the analytical solution for the problem, it is necessary to apply the inverse Laplace transform, i.e. $f(\tau)=\frac{1}{2\pi i}\int\limits_{\sigma_{1}-i\infty}^{\sigma_{1}+i\infty}\exp(s\tau)F(s)\mathrm{d}s$. Note that in contrast to (49), this approach can be readily applied for the last equation, subject to the restriction that has been imposed above. As a result, we obtain an analytic solution for the problem:

[TABLE]

Now we try to find a solution to the problem by the AI method. Taking into account (6) and (9), we obtain a system of equations for finding the zero-order solution:

[TABLE]

Then, the solution of the last has the form:

[TABLE]

Note that the image of the zero-order solution (54) has a similar form to (15) for a problem, which was detailed in Appendix A. To obtain the first amendment, we use a form such as $N_{1}=\eta^{-\alpha}\Phi_{1}(x)$, and when we insert $N_{0}$ into (2), we obtain an equation for the image of $N_{1}$:

[TABLE]

When we perform the same actions as were done in Appendix A while the obtaining $N_{1}$, we finally obtain

[TABLE]

If we expand the exponents in a series after retaining until the second order of $x$ and $x_{0}$ as it has been done in finding analytical solution, we obtain

[TABLE]

It should be noted that here we restricted ourselves to only the first amendment because the following amendments will contain $x$ and $x_{0}$ higher than the second order, as it can be seen in Appendix A.

It can be seen that every term of the image of the AI solution (E) coincides with the corresponding term of the analytical image (51). Of course, own solutions of these images are identical to each other, i.e. (52) is a solution for both (E) and (51). Furthermore, if we examine our problem till the first order $x$ and $x_{0}$, we see that the first three terms of the image (E) (that included in the $N_{0}$) coincide with the first three terms of image (51). This indicates that the common iterative solution (8) corresponds and converges with the exact solution. In other words, using more amendments by means of the AI method would bring it closer to the exact solution of the considered problem.