11institutetext: Université de Paris 13
† LIPN, ∗ LAGA,
99, avenue Jean-Baptiste Clément, 93430 Villetaneuse, France
<EMAIL_ADDRESS>

Authors’ Instructions

# Towards Parallel CFD computation
for the ADAPT framework

Imad Kissami${}^{\dagger}{}^{\ast}$ Christophe Cerin† Fayssal Benkhaldoun∗
Gilles Scarella∗

###### Abstract

In order to run Computational Fluid Dynamics (CFD) codes on large scale
infrastructures, parallel computing has to be used because of the
computational intensive nature of the problems. In this paper we investigate
the ADAPT platform where we couple flow Partial Differential Equations and a
Poisson equation. This leads to a linear system which we solve using direct
methods. The implementation deals with the MUMPS parallel multi-frontal direct
solver and mesh partitioning methods using METIS to improve the performance of
the framework. We also investigate, in this paper, how the mesh partitioning
methods are able to optimize the mesh cell distribution for the ADAPT solver.
The experience gained in this paper facilitates the move to the 3D version of
ADAPT and the move to a Service Oriented view of ADAPT as future work.

###### Keywords:

Unstructured mesh, Mesh partitioning, Parallel direct solver, Multi-frontal
method, MUMPS, METIS, Multi-physics, Multi-scale and Multilevel Algorithms.

##  1 Introduction and context of the work

In the last recent decades CFD (Computational Fluid Dynamics) used to play an
important role in industrial designs, environmental impact assessments and
academic studies. The aim of our work is to fully parallelize an unstructured
adaptive code for the simulation of 3D streamer propagation in cold plasmas.
As a matter of fact, the initial sequential version of the Streamer code needs
up to one month for running typical benchmarks.

The generation of streamer discharges is described by coupling electrostatic
to the motion of charged particles (electrons, positive and negative ions).
The electrostatic is represented by a Poisson equation for the electric
potential and the motion of particles is described by a set of convection-
diffusion-reaction equations.

We aim to parallelize both the linear solver issued from the Poisson equation
and the evolution equation using domain decomposition and mesh adaptation at
’the same time’. To our knowledge this is the first time that such a challenge
is considered. Intuitively speaking, separating the two steps may add delays
because we need to synchronize them and to align the execution time on the
slowest processor. Considering the two steps simultaneously offers the
potential to better overlap different computational steps.

Other authors have tackled the parallelization of either the linear solvers or
the evolution equation usually without mesh adaptation. In [1] the authors
consider the parallelization of a linear system of electromagnetic equation on
non adaptive unstructured mesh. Their time integration method leads to the
resolution of a linear system which requires a large memory capacity for a
single processor.

In [2] the authors introduce a parallel code which was written in C++
augmented with MPI primitives and the LIS (Linear Iterative Solver) library.
Several numerical experiments have been done on a cluster of 2.40 GHz Intel
Xeon, with 12 GB RAM, connected with a Gigabit Ethernet switch. The authors
note that a classical run of the sequential version of their code might easily
exceed one month of calculation. The improvement of the parallel version is
due to the parallelization of three parts of the code: the diagonalization of
the Schrödinger matrix, advancing one step in the Newton–-Raphson iteration,
and the Runge-–Kutta integrator.

The authors in [3] focused only on the parallelization of a linear solver
related to the discretization of self-adjoint elliptic partial differential
equations and using multigrid methods.

It should also be mentioned that authors in [4] have successfully studied
similar problems to those presented in this paper. The associated linear
systems have been solved using iterative Gmres and CG solvers [5]. One
difference is that in our work we consider direct methods based on LU
decomposition using the MUMPS solver.

Moreover, the direct solution methods generally involve the use of frontal
algorithms in finite element applications. The advent of multi—-frontal
solvers has greatly increased the efficiency of direct solvers for sparse
systems. They make full use of high performance software layers such as
invoking level 3 Basic Linear Algebra Subprograms (BLAS) [6, 7] library. Thus
the memory requirement is greatly reduced and the computing speed greatly
enhanced. Multi—-frontal solvers have been successfully used both in the
context of finite volume, finite element methods and in power system
simulations.

The disadvantage of using direct solvers is that the memory size increases
much more quickly than the problem size itself. To circumvent this problem,
out-of-core multi—-frontal solvers have been developed which have the
capability of storing objects of the resolution on the disk during
factorization. Another viable alternative is to use direct solvers in a
distributed computing environment. MUMPS [8, 9] is among the fastest parallel
general sparse direct solvers available under public domain.

In order to deal with complex geometries and fluid flows, a large number of
mesh cells should be used. Therefore the parallel computing paradigm has to be
introduced for coping with this large number of mesh cells. In the parallel
computing field, several factors, such as the load balancing, the number of
neighboring sub-domains and the halo cells (cells which are at the boundary of
sub-domains), affect the performance. A well balanced load has the potential
to reduce the amount of waiting processors and the partitioning method can
optimize the distribution of the mesh cells across processors, thus it will
improve the performance of the parallel application too.

In this paper, we also use METIS [10], an open source mesh partitioning
software, and we incorporate it into the streamer code of ADAPT [11], to solve
the evolution PDE, in order to study the impact of mesh partitioning on the
parallel version of ADAPT that we are currently developing.

Moreover we carry out experiments using the parallel multi—-frontal direct
solver (MUMPS) with matrices extracted from the streamer code of the ADAPT
platform. The general strategy is as follows. The linear system of equations
is evaluated on different processors corresponding to the local grid assigned
to the processor. The right hand side vector is assembled on the host
processor and it is injected into the MUMPS solver. At the last step of the
MUMPS solver, the solution is assembled centrally on the host processor. This
solution is then broadcast to all the processors. We discuss later the pro and
cons of such a strategy.

The organization of the paper is the following. In section 2 we introduce the
aim of the ADAPT platform and its positioning. In section 3 we discuss about
parallel approaches to solve evolution equation coupled with Poisson equation,
show the strategy to parallelize our code and do some experiments with METIS.
In section 4 we analyze different parts of our code in terms of speedup and
efficiency and we provide numerical results that show the efficiency our work.
Section 5 concludes the paper.

##  2 The ADAPT framework

###  2.1 Overview

ADAPT [11] is an object oriented platform for running numerical simulations
with a dynamic mesh adaptation strategy and coupling between finite element
and finite volume methods.

ADAPT, as a CFD (Computational Fluid Dynamics) software package, has been
developed for realizing numerical simulation on an unstructured and adaptative
mesh for large scale CFD applications.

In this paper, in order to tackle the long running time necessary for
numerical simulation, we study the MUMPS and METIS toolkits and we integrate
them into the ADAPT framework. In fact we perform code coupling between the
two main steps of the studied problem.

The existing ADAPT implementation is a sequential C++ code for each
phenomenon, and the code requires a huge CPU time for executing a 3D
simulation. For example, the 3D streamer code may run up to 30 days before
returning results, for this reason we decided to parallelize the code. This
paper is an important step into this direction.

###  2.2 Working environment

To realize all the experiments (sequential and parallel) we worked on the MAGI
cluster÷footnotehttp://www.univ-paris13.fr/calcul/wiki which is located at the
University of Paris 13, and also on the Ada cluster of
Idris111http://www.idris.fr, one of the national computing facility in France.

##  3 Parallel approach

In this paper, we parallelize the evolution equation coupled with the Poisson
equation, written as:

$\left\\{\begin{aligned} \frac{\partial u}{\partial t}+F(V,u)=S,\\\ \\\ \Delta
P=b.\end{aligned}\right.$ (1)

given that $F(V,u)=div(u.\overrightarrow{V})-\Delta u$, $S=0$, and
$V=\overrightarrow{\nabla}P$, the previous system gives :

$\left\\{\begin{aligned} \frac{\partial u}{\partial
t}+div(u.\overrightarrow{V})=\Delta u\\\ \Delta P=b.\end{aligned}\right.$ (2)

The first equation is discretized using the finite volume method on an
unstructured triangular mesh. The time-integration of the transport equation
is performed using an explicit scheme. The discretized form of Poisson
equation leads to a linear algebraic system. We obtain the following set of
equations:

$\left\\{\begin{aligned} u_{i}^{n+1}=u_{i}^{n}-\frac{\Delta
t}{\mu_{i}}\underbrace{\sum_{j=1}^{m}u_{ij}\overrightarrow{V}_{ij}\overrightarrow{n}_{ij}|\sigma_{ij}|}_{Rez\\_conv}+\frac{\Delta
t}{\mu_{i}}\underbrace{\sum_{j=1}^{m}\overrightarrow{\nabla}u_{ij}\overrightarrow{n}_{ij}|\sigma_{ij}|}_{Rez\\_dissip}\\\
A.\overrightarrow{P}^{n}=\overrightarrow{b}^{n}\end{aligned}\right.$ (3)

where m is the number of faces of volume $\mu_{i}$, $\overrightarrow{n}_{ij}$
is the unit normal vector of the face $\sigma_{ij}$ (face between volumes
$\mu_{i}$ and $\mu_{j}$) and $|\sigma_{ij}|$ is its length. Other variables
denoted by subscript ij represent variables on the face $\sigma_{ij}$.

A is a large sparse matrix, the coefficients of the matrix depend only on the
grid topology.

#### Computation of $u_{ij}$:

In figure 1 we took the example when the mesh is splitted into two subdomains,
and we compute $u_{ij}$ on the face $\sigma_{ij}$. For a given face
$\sigma_{k}$, suppose $T_{i}$ and $T_{j}$ are respectively the cells at the
left and right of $\sigma_{k}$. Let’s note $\sigma_{k}=\sigma_{ij}$ and
$u_{k}=u_{ij}$. The algorithm 1 shows how the $u_{ij}$ is computed on face
$\sigma_{ij}$.

Figure 1: Computation of $u_{ij}$ in mesh splitted into two subdomains

1 F:number of faces;

2 for  _k:=1 _to_ F_ do

3 if  _$dot(V_{k}.n_{k}) >=0$_ then

4 $u_{k}=u_{i}$;

5

6 end if

7 else

8 if  _$\sigma k$ is inner faces_ then

9 $u_{k}=u_{j}$;

10

11 end if

12 else if  _$\sigma k$ is halo faces_ then

13 $u_{k}=u_{h}$;/* $u_{h}:$ halo value sent by neighbor subdomain */

14 end if

15

16 end if

17

18 end for

Algorithm 1 Compute $u_{ij}$ on face $ij$

#### Computation of

$\overrightarrow{\nabla}u_{ij}$: In figure 3 we take the same mesh, and we
compute $\overrightarrow{\nabla}u_{ij}$ of the face $\sigma_{ij}$. The
algorithm 2 shows how the $\overrightarrow{\nabla}u_{ij}$ is computed on
volume $\mu_{ij}$. The diamond cell in 2 is constructed by connection of
centers of gravity $(i,j)$ of cells $T_{i}$, $T_{j}$ which shares the face
$\sigma_{ij}$ and its endpoints $A,B$.

Figure 2: Diamond cell in 2D

$\overrightarrow{\nabla}u_{ij}=\frac{1}{2\mu(D_{\sigma_{ij}})}[(u_{A}-u_{B})\overrightarrow{n}_{LR}|\sigma_{LR}|+(u_{j}-u_{i})\overrightarrow{n}_{ij}|\sigma_{ij}|]\\\
$ (4)

1$u_{node}$:double;

2 $N$ :number of nodes;

3 $Alpha$ :weight coming from the least square method;

4 for  _n:=1 _to_ N _ do

5 for  _c:=1 _to_ inner cells around node n_ do

6 $u_{node}(n)+=Alpha*u_{j}(c)$;

7

8 end for

9 for  _m:=1 _to_ halo cells around node_ do

10 $u_{node}(n)+=Alpha*u_{h}(m)$;

11

12 end for

13

14 end for

15return $u_{node}$;

Algorithm 2 Compute value at node ($u_{node}$)  Figure 3: Computation of
$\nabla{u_{ij}}$ in mesh splitted into two subdomains

1F:number of faces;

2 $\overrightarrow{\nabla}u$:Vector2d;

3 for  _k:=1 _to_ F _ do

4 A:first node of face;

5 B:second face node;

6 i:center of gravity of cell $T_{i}$;

7 j:center of gravity of cell $T_{j}$;

8 $mes=\frac{1}{2\mu(D_{\sigma_{ij}})}$;

9 if  _Inner faces_ then

10
$\overrightarrow{\nabla}u(k)=mes*(u_{node}(A)-u_{node}(B))\overrightarrow{n}_{LR}|\sigma_{LR}|+(u_{j}-u_{i})\overrightarrow{n}_{ij}|\sigma_{ij}|$;

11

12 end if

13 else if  _Halo faces_ then

14
$\overrightarrow{\nabla}u(k)=mes*(u_{node}(A)-u_{node}(B))\overrightarrow{n}_{LR}|\sigma_{LR}|+(u_{h}-u_{i})\overrightarrow{n}_{ih}|\sigma_{ih}|$;

15

16 end if

17

18 end for

19return $\overrightarrow{\nabla}u$;

Algorithm 3 Compute ${\nabla}u_{ij}$ on face $ij$

#### Computation of the linear system:

The system is solved directly by LU decomposition with an implementation for
sparse matrices. We use the Intel MKL library222https://software.intel.com/en-
us/intel-mkl#pid-3374-836 solvers and UMFPACK [12].

Note that in traditional approches, one splits the mesh to solve both the
evolution equation and the linear system in each sub-domain using iterative
methods with added communications. In our work we use the same methodology in
partitioning the mesh but the linear system is solved with a direct method.

#### Parallelization:

For the parallelization of CFD simulations, ADAPT employs the domain
decomposition method. The 2D unstructured mesh in Figure 4 was decomposed into
eight subdomains using the METIS algorithm, these partitions have
approximatively the same size that we may consider as a good property because
the workload will be balanced on the homogenous processors of our platform.

Figure 4: Decomposition of the computational domain into 8 subdomains using
the METIS method

Algorithm 4 on page 10 provides with a piece of pseudocode that shows how the
parallelization is done in the streamer code according to the coupling of
evolution equation with Poisson equation. We can see that most parts of the
code are parallel ones (line 19, line 9 to 11 and line 22 to 26) except
reading and splitting mesh in the beginning and between lines 14 and 15.
Indeed, at this step, we construct matrix A (eq. 2) that will be used to solve
the linear system in line 19. Matrix A is computed one time because it depends
only on the mesh.

Our contribution will serve in the future to tackle the 3D version of streamer
which includes the same type of equations, same strategy to solve the linear
system. The difference is the structure of the grid because we will work with
Tetrahedra instead of triangles, which will make the problem much more
complicated due to big differences in handling the mesh in the 3D case.

1 W:double;

2 if  _rank==0_ then /* for the master processor */

3 Read mesh data;

4 Split mesh with METIS;

5 Distribute mesh to all processors;

6

7 end if

8W=0;

9 for  _each rank_ do /* for each processor */

10 Initialize conditions and create constants;

11 Send the informations of halos cells to neighbor subdomains;

12 Apply boundary conditions;

13

14 end for

15if  _rank==0_ then

16 Construct matrix of linear system;

17 Split the matrix and send part of each processor;

18

19 end if

20for  _each iteration_ do

21 for  _all rank_ do /* for all processors */

22 Solve linear system using MUMPS;

23 end for

24 for  _each rank_ do

25 Send the informations of halos cells to neighbor subdomains;

26 Apply boundary conditions;

27 Compute fluxes of convection, diffusion and source term;

28 Update solution : $W^{n+1}=W^{n}+\Delta
t*(rez\\_conv+rez\\_dissip+rez\\_source)$;

29 Save results in parallel way using Paraview;

30

31 end for

32

33 end for

Algorithm 4 Algorithm of parallel ADAPT

##  4 Application to streamer equations

When non-ionized or low ionized matter is exposed to high intensity electric
field, non-equilibrium ionization processors (so-called discharges or
streamers) occur. Because of the reactive radicals they emit, streamers are
used for the treatment of contaminated media like exhaust gasses, polluted
water or bio-gas.

###  4.1 The governing equation

The streamer consists of a convection-diffusion-reaction for the electron
density, an ordinary differential equation for the positive ion density
coupled by the Poisson’s equation for the electric potential. The model is
given by the following equations:

$\displaystyle\frac{\partial n_{e}}{\partial
t}+div(n_{e}\overrightarrow{v_{e}}-De\overrightarrow{\nabla}n_{e})=S_{e},$ (5)
$\displaystyle\ \ \frac{\partial n_{i}}{\partial t}=S_{e},$ (6)
$\displaystyle\Delta V=-\frac{e}{\epsilon}(n_{i}-n_{e})\ ,$ (7)
$\displaystyle\ \ \vec{E}=-\overrightarrow{\nabla}V\ ,$ (8)

Where $V$ is a potential of electric field E, $\epsilon$ is the dielectric
constant, e the electron charge, $n_{e}$ and $n_{i}$ are the number densities
of electrons and positive ions, the drift velocity of electrons is
$v_{e}=v_{e}(E)$ and $D_{e}=D_{e}(E,v_{e})$ is the diffusion coefficient. The
source terms depend on the electron drift velocity and the electron density
$S_{e}=S_{e}(v_{e},n_{e})$.

###  4.2 Numerical method

The equations of the model are discretized using the finite volume method on
an unstructured triangular mesh. The time-integration of the transport
equations is performed using an explicit scheme. The discretized form of
Poisson’s equation consists of a linear algebraic system that is solved with a
direct method at each time step during the time-integration. We approximate
the equation for the electron density (eq. 3) by the following finite volume
method:

$\displaystyle\frac{\partial n_{e}}{\partial t}+\frac{1}{\mu(T)}\
\oint_{\partial
T}(n_{e}\overrightarrow{v_{e}}-D_{e}\overrightarrow{\nabla}n_{e})\overrightarrow{n}\mathrm{d}s=S_{e}$
(9)

where $\mu(T)$ is the volume of the cell T, $\overrightarrow{n}$ the outward
unit normal vector to the faces of the cell T.

The Poisson’s equation (eq. 5) is discretized by a central type approximation
which leads to a system of linear equation, as follows:

$\displaystyle A.\overrightarrow{V}^{n+1}=\overrightarrow{b}^{n+1}$ (10)

$A$ is a matrix of coefficients, $\overrightarrow{V}$ is a vector of unknowns
(its dimension is equaled to the total number of cells) and
$\overrightarrow{b}$ is a vector of right hand side.

###  4.3 Parallel results

Table 1 summarize our results. The speedup in this example shows a good
scalability of the present method for solving the convection-diffusion
equation coupled with Poisson equation problems using mesh with 529240 cells.
For practical applications, the computation time could be reduced from
49h54min (one computing core) down to 5min (1024 computing cores). This test
is made on MAGI cluster at Paris13.

Table 1: Execution time (in s) of different parts of parallel 2D streamer code using mesh with 529240 cells Compute | Total | Convection | Diffusion | Linear solver
---|---|---|---|---
cores |  |  |  |
1 | 49 h 54 min 48 s | 02 h 51 min 04 s | 13 h 06 min 00 s | 33 h 57 min 44 s
2 | 25 h 06 min 27 s | 01 h 22 min 57 s | 06 h 41 min 02 s | 17 h 02 min 27 s
4 | 12 h 34 min 35 s | 00 h 42 min 07 s | 03 h 22 min 04 s | 08 h 30 min 24 s
8 | 06 h 27 min 18 s | 00 h 22 min 13 s | 01 h 46 min 26 s | 04 h 18 min 38 s
16 | 03 h 39 min 45 s | 00 h 12 min 37 s | 01 h 01 min 40 s | 02 h 25 min 26 s
32 | 01 h 50 min 50 s | 00 h 08 min 03 s | 00 h 29 min 17 s | 01 h 13 min 29 s
64 | 01 h 01 min 41 s | 00 h 03 min 59 s | 00 h 17 min 05 s | 00 h 40 min 36 s
128 | 00 h 32 min 44 s | 00 h 01 min 52 s | 00 h 08 min 22 s | 00 h 22 min 29 s
256 | 00 h 18 min 16 s | 00 h 01 min 02 s | 00 h 04 min 34 s | 00 h 12 min 39 s
512 | 00 h 10 min 07 s | 00 h 00 min 33 s | 00 h 02 min 52 s | 00 h 06 min 42 s
1024 | 00 h 05 min 14 s | 00 h 00 min 18 s | 00 h 01 min 33 s | 00 h 03 min 23 s

#### Scalability measurement:

There are two basic ways to measure the scalability of parallel algorithms:

  1. (a)

Strong scaling with a fixed problem size but the number of processing elements
are increased.

  2. (b)

Weak scaling with problem size and compute elements increasing concurrently.

The speedup is defined as

$sp=\frac{t_{b}}{t_{N}},$

and the ideal speedup, $sp_{ideal}$, is naturally equivalent to the number of
compute cores. Therefore, the strong scaling efficiency can also be calculated
by

$\frac{s_{p}}{sp_{ideal}}$

which shows how far is the measured speedup from the ideal one.

In the this work, the strong scaling measurement is employed and the scaling
efficiency is calculated by:

$\frac{t_{b}N}{t_{N}N_{b}}100,\,N\geq N_{b}$

where $t_{b}$ is the wall clock time of a base computation, $N_{b}$ is the
number of compute cores used for the base computation, $t_{N}$ is the wall
clock time of a computation with $N$ compute cores.

Figures 5 and 6 illustrate the speedup and the strong scaling efficiency in
the computations of our illustrate example with $t_{b}$ the wall clock times
by one compute core. A more detailed injection into the data of this example
provides interesting information about the computational efficiency of the
individual parts of the parallel simulation. The results of the strong scaling
tests, plotted in Fig 6 depict that the parallel speedup of the different
parts increases up to 1024 computational cores, so our parallel strategy shows
its efficiency.

Figure 5: Speedup of different parts of streamer code using mesh with 529240
cells Figure 6: Efficiency of different parts of streamer code using mesh with
529240 cells Figure 7: Plasma discharge in sequential code (left) and parallel
one (right)

###  4.4 Comparison between sequential and parallel code

| Compute | Execution
---|---|---
| cores | time
Mesh using 43674 elements | 1 | 04h30min
Mesh using 529240 elements | 1 | 150h15min
64 | 03h59min
Table 2: Comparison between sequential and parallel code

In Figure 7 a comparison is made between sequential 2D streamer code and
parallel counterpart. We can also observe the propagation of streamer in two
cases: more the mesh is fine more the results are accurate (without
oscillations). In table 2, the parallel code takes the same running time
(using 64 MPI processors) as the sequential code but the input size is 12
times greater.

To conclude, our strategy although not classical shows to be very efficient,
and provides with relevant properties regarding performance.

##  5 Conclusions and future work

In this work, we have introduced an effective parallelization of the ADAPT
platform for CFD applications. We present the parallelization of the
convection-diffusion equation and the linear system. The originality of this
work is to solve linear system using a direct method whereas existing studies
mainly use iterative methods to solve linear system in this kind of problems.
The workflow is realized: (1) using external tools such as METIS for mesh
partitioning in order to enable the computational load balance for the global
assembly, (2) using MPI to assure communication, (3) using MUMPS solver to
solve the system of linear equations. The current approach with MUMPS shows
significant advantages in terms of avoiding the problems of pre-conditioners,
and doesn’t need a lot of iterations to converge to the solution. In summary,
the most important advantages of the presented scheme are:

  * •

Efficient memory usage distribution over computer nodes,

  * •

Good scalability to solve convection-diffusion equation,

  * •

Stable, fast and good scalability of linear solvers provided by MUMPS.

We observe that the campaign of numerical tests provides expected results,
meaning that our parallel implementation goes into the right direction at the
moment.

For the future, we plan to investigate in depth the mesh partitioning problems
and to conduct experiments with specific tools. We may target PAMPA [13] and
the question is how such tool can impact the overall performance of ADAPT
software. The challenge is to mix mesh partitioning and dynamic mesh
adaptations. How to manage load balancing in this context? When dealing with
three dimensional problems, the question is even more critical. At least, as
explained in the discussion section our current work will serve for solving 3D
streamer problems and a service oriented view of the ADAPT framework.

## Acknowledgment

The funding supports of this work is the EnCoMix AAP SPC project
(ANR-11-IDEX-05-02 Ref: SPC/JFG/2013‐031). The experiments conducted in this
work were done on the nodes B500 of University Paris 13 MAGI Cluster and
available at http://www.univ-paris13.fr/calcul/wiki/

## References

  * [1] Franck Assous, J. Segré, and Eric Sonnendrücker.  A domain decomposition method for the parallelization of a three-dimensional maxwell solver based on a constrained formulation.  Mathematics and Computers in Simulation, 81(11):2371–2388, 2011\.
  * [2] Francesco Vecil, José M. Mantas, María J. Cáceres, Carlos Sampedro, Andrés Godoy, and Francisco Gámiz.  A parallel deterministic solver for the schrodinger-poisson-boltzmann system in ultra-short dg-mosfets: Comparison with monte-carlo.  Computers & Mathematics with Applications, 67(9):1703 – 1721, 2014\.
  * [3] Yvan Notay and Artem Napov.  A massively parallel solver for discrete poisson-like problems.  Journal of Computational Physics, 281(0):237 – 250, 2015.
  * [4] Pawan Kumar, Stefano Markidis, Giovanni Lapenta, Karl Meerbergen, and Dirk Roose.  High performance solvers for implicit particle in cell simulation.  Procedia Computer Science, 18:2251–2258, 2013.
  * [5] Youcef Saad and Martin H. Schultz.  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [6] L. S. Blackford, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry, M. Heroux, L. Kaufman, A. Lumsdaine, A. Petitet, R. Pozo, K. Remington, and R. C. Whaley.  An updated set of basic linear algebra subprograms (blas).  ACM Transactions on Mathematical Software, 28:135–151, 2001.
  * [7] Ronald F. Boisvert and Jack Dongarra.  Preface to the special issue on the basic linear algebra subprograms (BLAS).  ACM Trans. Math. Softw., 28(2):133–134, 2002.
  * [8] Patrick Amestoy, Iain S. Duff, Jean-Yves L’Excellent, and Jacko Koster.  MUMPS : A general purpose distributed memory sparse solver.  In Tor Sørevik, Fredrik Manne, Randi Moe, and Assefaw Hadish Gebremedhin, editors, Applied Parallel Computing, New Paradigms for HPC in Industry and Academia, 5th International Workshop, PARA 2000 Bergen, Norway, June 18-20, 2000 Proceedings, volume 1947 of Lecture Notes in Computer Science, pages 121–130. Springer, 2000.
  * [9] Patrick Amestoy, Abdou Guermouche, Jean-Yves L’Excellent, and Stéphane Pralet.  Hybrid scheduling for the parallel solution of linear systems.  Parallel Computing, 32(2):136–156, 2006.
  * [10] George Karypis and Vipin Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM J. Sci. Comput., 20(1):359–392, December 1998.
  * [11] F. Benkhaldoun, J. Fort, K. Hassouni, and J. Karel.  Simulation of planar ionization wave front propagation on an unstructured adaptive grid.  Journal of Computational and Applied Mathematics, 236:4623–4634, 2012.
  * [12] Timothy A. Davis.  Algorithm 832: UMFPACK V4.3 - an unsymmetric-pattern multifrontal method.  ACM Trans. Math. Softw., 30(2):196–199, 2004.
  * [13] Cédric Lachat, François Pellegrini, and Cécile Dobrzynski.  PaMPA: Parallel Mesh Partitioning and Adaptation.  In 21st International Conference on Domain Decomposition Methods (DD21), Rennes, France, June 2012. INRIA Rennes-Bretagne-Atlantique.


# Semblance study of force-field and convection-diffusion solutions with
observations of solar magnetic phenomena

M. Enríquez-Vargas  Instituto de Ciencias Físicas, Universidad Nacional
Autónoma de México,
62210, Cuernavaca, Morelos, México  Jorge U. Pérez Romero  Facultad de
Ingeniería, Universidad Nacional Autónoma de México,
04510, Ciudad de México, México

###### Abstract

A quantitative analytical model based on the semblance method between the
modulation factor with solar phenomena is proposed. Different Local
Interstellar Spectra (LIS) have been computed to introduce into a transport
equation solution which in turn have been introduced into the atmospheric
yield function (Caballero-Lopez & Moraal, 2012), that allows to compute a
Cosmic Rays (CR) Modulation Factor. The results were as expected; there are
correlation between modulation factor and sunspots, and anticorrelation
between modulation factor and mean magnetic field.

A transport equation’s solution is necessary to compute atmospheric yield
function, in this case the used transport equation’s solutions were
convection-diffusion and force field. Both solutions offer similar models, yet
the force field solution shows a higher correlation value in the semblance
than the convection-diffusion solution.

Several LIS were also computed because they are introduced into the transport
equation solutions. The used LIS were Lagner, Potgieter and Webber LIS in
2003, Burguer and Potgieter LIS in 2000, Garcia-Munoz, Mason and Simpson LIS
in 1975 and Ghelfi, Barao, Derome and Maurin LIS in 2017. Those LIS were used
because they have a model for different nuclear composition: Helium and
Hydrogen. The LIS with more changes when is introduced into the semblance is
Garcia-Munoz, Mason and Simpson in 1975.

The atmospheric yield function also need data from some neutron detectors that
have changes according their cut-off rigidity value, hence two analyses were
made. In the first analysis were selected three stations from different
locations: close from the Arctic Circle as OULU station, in the Antarctic as
SNAE and close of the equator as AATB station. In the second analysis, five
different stations were averaged by arithmetic mean to observe and analyze
what is the general behaviour of arrival cosmic rays on the Earth. The
computed semblance models are very similar with significant differences that
are shown in an analysis of zeros.

Many images were produced when the analyses were computed, then the main
images, that are used as examples, are shown along the paper, while the rest
can be found in the appendix.

Cosmic Rays — Convetion-Diffusion solution — Force field solution — Semblance
— Geophysics — Local interstellar spectal — Atmospheric yiel funtion

111Computational resources can be found in https://github.com/JorUrie/TRC.git

##  1 Introduction

To generate the proposed model, data from neutron monitors were utilized and
diverse local interstellar spectra were computed. These spectra were
introduced into the atmospheric yield function to determine the cosmic ray
modulation factor.

Cosmic rays are electrically charged particles capable of reaching velocities
as high as the speed of light, in other words, they are atoms, that lack
electrons and contain high amounts of energy. Primary cosmic radiation is
composed of protons, $\alpha$-particles (two protons and two neutrons), $\
{}^{-4}$He and heavy nuclei.

Cosmic rays enter to the Earth’s magnetosphere and other layers from all
directions (primary cosmic rays). The primary CRs collide with particles in
the Earth’s atmosphere, these violent collisions cause particles to break down
into smaller particles (secondary cosmic rays). This process is repeated
multiple times until the particles lose sufficient energy to become elemental
particles or are unable to break down due to energy loss. This phenomenon is
known as an air shower (Auger et al., 1939).

High-energy cosmic rays are produced by stochastic collisions within a
magnetically turbulent plasma such as a supernova, and there are also low-
energy cosmic rays from similar stars to our Sun (Sekido & Elliot, 1985).

Cosmic rays originating from sources beyond our solar system are detected on
Earth’s surface, and their detection is affected by solar activity, which is
reliant on the solar cycles. Furthermore, the Sun also generates cosmic rays
within its corona.

The semblance has been used to compute the desired model. Semblance is a
geophysical tool mainly used in exploration seismology and it can be defined
as a correlation in Fourier terms or other transform between two traces, that
is, semblance is a measure to compute the similarity between two traces.

##  2 Transport equation

The used solutions of transport equation in this document is based in
Caballero-Lopez and Moraal (Caballero-Lopez & Moraal, 2004) article.

The Parker’s equation is given by:

$\frac{\partial f}{\partial t}+\frac{1}{4\pi
p^{2}}\nabla\cdot\vec{S}+\frac{1}{p^{2}}\frac{\partial}{\partial
p}(p^{2}\dot{p})f=Q$ (1)

where $S$ is the current differential density:

$\vec{S}=4\pi p^{2}(C\vec{V}f-K\cdot\nabla f)$ (2)

The equation 1 is in terms of the omnidirectional dis- tribution function
$f(r;P;t)$, $p$ is the particle’s moment, $t$ is the time, $r$ is the radial
distance, $Q$ is the source of cosmic rays and $V$ the solar wind speed. $K$
is the diffusion tensor, that contains a parallel component $k_{\parallel}$,
perpendicular $k_{\perp}$ and transverse $k_{T}$ and they describe the drifts.
In the equation 2, $C$ (from the equation 3) is known as the Compton-Getting
factor (Gleeson & Axford, 1968). This factor was introduced into the transport
equation, because the flow, which is isotropic in the reference system at rest
and is not the same in the solar wind system.

$C=-\frac{1}{3}\frac{\partial lnf}{\partial lnp}$ (3)

The adiabatic energy loss is denoted as:

$\dot{p}=-\frac{p}{3}\nabla\cdot\vec{V}$ (4)

Substituting 3 and 4 in 1 two possible solutions could be gotten:

$\begin{split}\frac{\partial f}{\partial
t}+\nabla\cdot(\underbrace{C\vec{V}f}_{a}\underbrace{-K\cdot\nabla
f}_{b})+&\underbrace{\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}\left(\frac{p}{3}\vec{V}\cdot\frac{\nabla
f}{f}\right)\right]}_{c}f=Q\end{split}$ (5)

##  3 Force Field

In the equation 5 the term $a$ corresponds to the convection, $b$ to the
diffusion and drifts of the cosmic rays and $c$ are the changes of energy. To
solve the previous equation, different approximations have been considered. If
steady state, without sources and neglect the energy losses is assumed, from
the equation 5 can be obtained:

$CVf-K\cdot\nabla f=cte=0$ (6)

For the equation 6 is valid, the term of adiabatic losses of the equation 7 is
required to be smaller than the terms convective and diffusive:

  * •

$Vr=k<<1$ valid for high energy cosmic rays, since k is proportional to
energy.

  * •

$(1/f)\partial f/\partial r<<C/r$ valid for small $r$, this means, in the
internal heliosphere.

Considering spherical symmetry, 8 is:

$CVf-k\frac{\partial f}{\partial r}=0$ (7)

The force field equation can be rewrite in terms of rigidity knowing that
$P=cp/q$:

$\ \frac{\partial f}{\partial r}+\frac{VP}{3k}\frac{\partial f}{\partial P}=0$
(8)

From:

$df(r,P)=\frac{\partial f}{\partial r}dr+\frac{\partial f}{\partial P}dP$ (9)

dividing by $dr$:

$\ \frac{df}{dr}=\frac{\partial f}{\partial r}dr+\frac{\partial f}{\partial
P}\frac{dp}{dr}$ (10)

Considering the equation 8:

$\ df=0$ (11)

The equation 11 implies that $f(r,p)$ is constant, equal to its value at the
boundary, along an outline with characteristic equation $dP/dr=VP/3k$, in the
space $(r,P)$. In the second term of the equation 8, the quantity $VP/k$ has
units of potential per unit length, that is, field units, that is why it is
called a force field solution. Considering the diffusion coefficient with the
form:

$\ k(r,P)=\beta k_{1}(r)k_{2}(P)$ (12)

the solution of the characteristic equation remains as:

$\begin{split}\int_{P}^{P_{b}(r,P)}{\frac{\beta\left(P^{{}^{\prime}}\right)k_{2}\left(P^{{}^{\prime}}\right)}{P^{{}^{\prime}}}dP^{{}^{\prime}}}&=\int_{r}^{r_{b}}{\frac{V\left(r^{{}^{\prime}}\right)}{3k_{1}(r^{{}^{\prime}})}dr^{\prime}}\equiv\phi(r)\end{split}$
(13)

where $\phi$ is the force field parameter. When $k2\propto P$ and
$\beta\approx 1$, the solution is reduced to the most commonly form utilized:

$\ \phi=P_{b}-P$ (14)

##  4 Convection-Diffusion

From the equation 5, $Q$ is not the adiabatic loss in a rest system:

$Vf-k\frac{\partial f}{\partial r}=0$ (15)

The outcome is:

$f=f_{b}e^{-M}$ (16)

Where $M=\int_{r}^{r_{b}}\frac{Vdr}{k}$ is the modulation function, $\phi$ is
the force field parameter. Its outcome is:

$M=\frac{3\phi}{\beta k_{2}}$ (17)

Where $M$ does not have unites.

Both solutions are introduced into the atmospheric yield solution, which is
necessary to compute the modulation factor. However LIS is also necessary in
order to compute the atmospheric yield solution.

##  5 Local Interstellar Spectral

The Local Interstellar Spectrum is a mathematical model created from
observations or simulations that describe the arriving of the Cosmic Rays at
outside the Heliosphere in energy terms.

The used LIS in these analyses were Lagner, Potgieter & Webber LIS in 2003,
Burguer & Potgieter LIS in 2000, Garcia-Munoz, Mason & Simpson LIS in 1975 and
Ghelfi, Barao, Derome & Maurin LIS in 2017.

###  5.1 Lagner, Potgieter & Webber LIS in 2003

This LIS is based on analytic solution of the Parker equation (Parker, 1965),
in Voyager 1, Voyager 2 and Pioneer data at 70 Astronomical Units (AU) close
to the Terminal Shock (TS) during minimum solar between 1987 and 1997.

The modulation developed by Parker included drift effects, that are related to
the polarity of solar magnetic field:

$\begin{split}\frac{\partial f}{\partial t}=&-\left(V+\left\langle
v_{D}\right\rangle\right)\cdot\nabla f+\nabla\cdot\left(K_{S}\cdot\nabla
f\right)\\\ &+\frac{1}{3}(\nabla\cdot V)\frac{\partial f}{\partial
InP}+j_{source}\end{split}$ (18)

Where $j_{source}$ is the local source. Into the equation was introduced
parallel, perpendicular, and asymmetric drift coefficients to describe the
gradient and derivates of curvature in large-scale Heliosphere Magnetic Field
(HMF).

The local sources are not considered and dependent time into spherical
coordinates system like a combination model of diffusive shock acceleration
and drift modulation in two spatial dimensions, neglecting any azimuthal
dependence and contemplating that is symmetric around the equatorial field
(Langner et al., 2003). The model offers several CR intensities that are
interpreted as particles made of different composition. When the modulation is
close or outside from TS is computed in rigidity terms. The protons model is
the following:

$J_{H}=\frac{2.1E_{\frac{k}{n}}^{-2.8}}{1+5.85E_{\frac{k}{n}}^{-1.2}+1.18E_{\frac{k}{n}}^{-2.54}}$
(19)

And the Helium model is:

$J_{H}=\frac{1.075E_{\frac{k}{n}}^{-2.8}}{1+3.9E_{\frac{k}{n}}^{-1.09}+0.90E_{\frac{k}{n}}^{-2.54}}$
(20)

###  5.2 Burguer & Potgieter LIS in 2000

The LIS was compared to ULYSSES data from September 1994 to July 1995 with a
Fast Scanning of Latitude (FSL). The Parker equation was analytically solved
considering an omnidirectional distribution (two dimensions), the source is
located 100 UA, the Solar Wind’s (SW) velocity is $400\ \frac{km}{s}$ in
equatorial plane and is increasing to $800\frac{km}{s}$ in polar zones. The
HMF angle sheet is $15{{}^{\circ}}$, which is a good value to fit for solar
modulation. The model was computed by Bibber (Bieber, 1999):

$f_{IS}(R)=\left\\{\begin{array}[]{c}1.9\times 10^{4}R^{-2.78}\ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ If\ R\geq 7GV\\\
\exp(9.472-1.999R-0.6938R^{2}+0.2988R^{3}-\\\ 0.04714R^{4})\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ If\ R<7\ GV\end{array}(\right)$ (21)

The model can also be expressed in terms of proton composition:

$J_{H}=\frac{1.9\times 10^{4}P^{-2.78}}{1+0.4866P^{-2.51}}$ (22)

$J_{He}=\frac{\left(3.8\times
10^{4}P^{-2.78}\right)}{\left(1+0.9732P^{-2.51}\right)}$ (23)

###  5.3 Garcia-Munoz, Mason and Simpson LIS in 1975

To compute this LIS was consulted IMP-5, IMP-7, and IMP-8 data, because the
satellites can capture several particles and their energy rage is wider, then
the select data were electron > $100MeV$ with quiet-time measurements.

The Parker’s equation was analytically solved considering symmetric spherical
coordinates system in stationary state and drawing on Fisk method (Fisk et
al., 1974). The LIS is the following:

$j=A(T+Bexp(-CT))^{\gamma}$ (24)

Where $j$ is particle flux, $T$ is kinetic energy and parameters are showing
in the table 1:

Species | A | B | C | $\gamma$
---|---|---|---|---
H | $9.9\times 10^{8}$ | 780 | $2.5\times 10^{-}4$ | 2.65
${}^{4}He$ | $1.4\times 10^{8}$ | 660 | $1.4\times 10^{-}4$ | 2.77
C | $1.8\times 10^{6}$ | 620 | $5.2\times 10^{-}4$ | 2.68
Table 1: Equation 7 parameters Garcia-Munoz et al. (1975)

###  5.4 Ghelfi, Barao, Derome and Maurin LIS in 2017

Used data were AMS-01 and AMS-02, BESS-Polar and PAMELA, because they were
recent (in that epoch), but Voyager 1 data was also considered.

First, the author purposes the simplest modulation model to link unmodulated
(IS) to modulated (TOA) qualities, that is a force field approximation Ghelfi
et al. (2017):

$\frac{E^{TOA}}{A}=\frac{E^{IS}}{A}-\frac{|Z|}{A}\phi$ (25)
$J^{TOA}\left(E^{TOA}\right)=\left(\frac{p^{TOA}}{p^{IS}}\right)^{2}\times
J^{IS}\left(E^{IS}\right)$ (26)

Where $E$ is total energy, $p$ is the momentum, and
$J\equiv\frac{dJ}{dE_{\frac{k}{n}}}$ is the differential flux per kinetic
energy per nucleon $E_{\frac{k}{n}}$.

An analysis $\chi^{2}$ was necessary to fix the TOA flux with all the species
$N_{j}(i)$ in this $t_{j}$, over all possible energy $E_{k}(i,\ j)$ and is
giving by:

$\chi^{2}=\sum_{t_{i}}{\sum_{N_{j(i)}}{\sum_{E_{k(i,j)}}\frac{\left(J^{TOA}\left(J_{j}^{IS},\
\phi_{j},E_{k}\right)-data_{ijk}\right)^{2}}{\sigma_{ijk}}}}$ (27)

Where IS parameters are free. For the previous mentioned was used power laws
in total energy (O’Neill, 2006) or rigidity (Shikaze et al., 2007). TOA data
for H and He are old and have more inconsistencies than new data, thus has
applied a $\chi^{2}$ analysis to minimize errors, therefore was used package
MINUIT minimisation package (James, 1998) from Root CERN libraries Sartini et
al. (2010). Markov Chain Monte Carlo analysis (MCMC) was implemented from the
GreAT package (Putze & Derome, 2014) to determine the correlation between
credible intervals for $\phi$ and IS flux. MCMC analysis is based on Bayes
theorem. Finally, all the obtained data is represented as an interpolated
logarithmic-polynomial function, which is shown in the equation 28.

$\begin{split}\log_{10}\left(J_{IS}\right)=\left\\{\begin{array}[]{ll}\sum_{i=0}^{12}c_{i}\left(\frac{\log_{10}\left(E_{\frac{k}{n}}\right)}{\log_{10}\left(800_{\frac{GeV}{n}}\right)}\right)^{i}&\text{if
}E_{\frac{k}{n}}<800\frac{GeV}{n}\\\\[10.0pt]
C_{0}-C_{1}\left(\frac{\log_{10}\left(E_{\frac{k}{n}}\right)}{\log_{10}\left(800_{\frac{GeV}{n}}\right)}\right)&\text{if
}E_{\frac{k}{n}}\geq 800\frac{GeV}{n}\\\ \end{array}\right.\end{split}$ (28)

In this paper were just applied four LIS, however more models were also
considered, and it is worthwhile to mention them.

###  5.5 Boschini, Della Torre, Gervasi and more LIS in 2018

For the calculation of this LIS was used GALPROP and HeldMod to simulate CR
data, but real data was required, so Voyager 1, BESS, PAMELA, AMS-01 and
AMS-02 data were consulted.

GALPROP code use astronomical, particle physics and nuclear information to
predict CR flux, X-ray, synchrotron emissions and their polarisation (Strong
et al., 2007). Heliosphere propagation was computed by GALPROP and thanks to
MCMC analysis, data was more realistic, after the information was introduced
in HelMod program to obtain modulated data.

Previous data was compared with real information from the satellites with
different energy content. The used information was a combination of Voyager 1,
AMS-02, CREAM-I and ATIC-02, but CREAM-I and AMS-02 had better adjust because
they offer minimal error. To compute the analytical LIS was used MCMC analysis
from Eureqa222http://www.nutonian.com/products/eureqa/ and the analytical
model take the next form:

$F(R)\times R^{2.7}=\left\\{\begin{array}[]{c}\sum_{i=0}^{5}{a_{i}R^{i}}\ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ R\leq 1GV\\\
b+\frac{c}{R}+\frac{\left(d_{1}\right)}{d_{2}+R}+\frac{\left(e_{1}\right)}{e_{2}+R}+\frac{\left(f_{1}\right)}{f_{2}+R}+\\\
gR\ \ \ \ \ \ \ R\geq 1GV\\\ \end{array}\right.\ $ (29)

Where the parameters are showing in the table 2.

Figure 1: Table 2. Equation 12 parameters (Boschini et al., 2017)

###  5.6 Vos & Potgieter LIS in 2015

The LIS was computed using observations from several PAMELA data during a
solar minimum of 2006 to 2009. This LIS was computed from the Parker equation
(Parker, 1965), which was solved as a ray distribution function in rigidity
terms, time, position in 3D, in a heliocentric spherical coordinates where
polar angles is equal to $90^{\circ}$. All short-term modulation effects are
neglected (solar minimum conditions), average particle drift velocity is
caused by gradients and curvature in Heliospheric Magnetic Field (HMF), there
is a symmetric diffusion tensor and adiabatic energy changes except in
heliosheath (Potgieter et al., 2012).

The LIS is the following (Vos & Potgieter, 2015):

$j_{LIS}=2.70\left(\frac{E^{1.12}}{\beta^{2}}\right)\left(\frac{E+0.67}{1.67}\right)^{-3.93}$
(30)

Where $E$ is kinetic energy, $\beta=\frac{v}{c}$ particle velocity (like light
velocity), $j_{LIS}=P^{2}f$.

###  5.7 Moskalenko, Strong, Ormes and Potgieter LIS in 2002

Cosmic Ray data was simulated by DTUNUC Monte Carlo code (Ferrari et al.,
1996) and data was recollected from two previous works (Tan & Ng, 1983) and
(Tan & Ng, 1983). The Parker equation (Parker, 1965) was solved through Crank-
Nicholson numerical solution (Crank & Nicolson, 1947) and accorded with the
author a model more realistic was sought, then leaky box code was applied
(weighted-slab). GALPROP was also used to compute 3D cosmic rays.

The Helium LIS was approximated to a Force Field, because of low energy
approximation of the $\phi$ potential modulation has been selected using the
CLIMAX neutron monitor.

To ensure that the spectrum is similar to the simulated values, a power law
dependent on kinetic energy based on the convection-diffusion solution has
been approximated, which cannot only be described by a single function.

Finally, the LIS was subjected to a $\chi_{n}^{2}$ function to observe the
quality of the adjusted data for each individual measurement (Moskalenko et
al., 2002) is just shown in the equation 31:

$\begin{split}J_{H}=\left\\{\begin{array}[]{l}\exp\left(4.64-0.08\left(\log\left(E_{\frac{k}{n}}\right)\right)^{2}-2.91\sqrt{E_{\frac{k}{n}}}\right),\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ E_{\frac{k}{n}}\leq
1GeV\\\
\exp\left(3.22-2.86\log\left(E_{\frac{k}{n}}\right)-\frac{1.5}{E_{\frac{k}{n}}}\
\right),\\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
E_{\frac{k}{n}}>1GeV\\\ \end{array}\right.\ \end{split}$ (31)

The LIS previously mentioned are shown in the figure 2 into the rigidity
purpose parameters by each author:

Figure 2: Local Interstellar Spectral purpose by (Vos & Potgieter, 2015),
(Ghelfi et al., 2017), (Burger et al., 2000), (Boschini et al., 2017),
(Garcia-Munoz et al., 1975), (Moskalenko et al., 2002) and (Langner et al.,
2003)

##  6 Semblance

The semblance is the cross-relation between two traces, the computed value is
interpreted as the similarity of two points contained in the signals that lie
in the Fourier domain, however in this analysis was used Continuous Wavelet
Transform (CWT) that represents a wave in two dimensions and measure the
similarity between two signals.

CWT is given:

$CWT(u,s)=\int_{-\infty}^{\infty}{h(t)\frac{1}{|s|^{0.5}}\Psi^{*}\left(\frac{t-u}{s}\right)dt}$
(32)

Where $s$ is scale, $u$ is displacement, $\Psi$ is the mother wavelet , and ∗
means complex conjugate and $t$ is the coordinate in time (Cooper & Cowan,
2008) the equation 32 can also be expressed in spatial terms (Teolis, 2017):

$\Psi(x)=\frac{1}{\pi f_{b}}e^{2\pi if_{c}x}e^{-\frac{x^{2}}{f_{b}}}$ (33)

The conventional semblance is a normalized coherence measure, first computed
in 1969 (Taner & Koehler, 1969) A time later, Neidell and Taner (Neidell &
Taner, 1971) found, that the coherence and semblance were two different
operators, so semblance was expressed as shown in the equation 34:

$S_{NT}\
[i]=\frac{\sum_{j=i-M}^{i+M}\left(\sum_{k=0}^{N-1}{q[j,k]}\right)^{2}}{N\sum_{j=i-M}^{i+M}{\sum_{k=0}^{N-1}{q[j.k]^{2}}}}$
(34)

where $i,\ j$ are time sample indices, $k$ is a trace number, $q[j,\ k]$ is
the trace amplitude at time index $j$ a trace number $k$ of the NMO-corrected
gather.

To reduce and make more smoothly the decays often is applied a boxcar filter,
that can be rewritten as (Luo et al., 2015):

$S_{C}\
[i]=\frac{\sum_{j}{h[i-j]\left(\sum_{k}{q[j,k]}\right)^{2}}}{N\sum_{j}{h[i-j]\sum_{k}{q[j.k]^{2}}}}$
(35)

There is an alternative expression for the conventional semblance, that can be
expressed as a normalized correlation coefficient, but first a trace $r[j]$
should be defined as a summation over trace number:

$r[j]\equiv\sum_{k}{q[j,k]}$ (36)

Therefore:

$C_{rq}[i]\equiv\sum_{j}{h[i-j]\sum_{k}{r[j]q[j,k]}}$ (37)
$C_{rr}[i]\equiv\sum_{j}{h[i-j]\sum_{k}{r[j]^{2}}}$ (38)
$C_{qq}[i]\equiv\sum_{j}{h[i-j]\sum_{k}{q[j,k]^{2}}}$ (39)

Now, the three equations 37, 38 and 39 can be rewritten as a conventional
semblance (Luo et al., 2015):

$S_{C}[i]=\frac{C_{rq}[i]^{2}}{C_{rr}[i]C_{qq}[i]}$ (40)

##  7 Count rates from ground-based detectors

The neutrons detected at ground level are the product of atmospheric showers
caused by the collision of a primary particle, that collides with a molecule
in the atmosphere. The neutron count is given by the following equation:

$\begin{split}N\left(P_{c},\ \ x,\ \ t\right)=&\
\int_{P_{c}}^{\infty}\left(\frac{-dN}{dP}\right)dP=\\\ &\sum_{i\
}{\int_{P_{c}}^{\infty}{S_{i}(P,\ x)j_{i}(P,t)dP}}\end{split}$ (41)

Where $j_{i}(P,t)$ is the spectrum of the primary species above the
atmosphere, and $S_{i}(P,x)$ is the atmospheric yield function due to this
species. $P_{c}$ is the cut-off rigidity, which is the necessary minimum
energy for a particle to enter the Earth’s magnetosphere. The quantity $dN=dP$
is the differential counting rate of the instrument inside the atmosphere. The
yield function, that was used in this analysis is the following (Caballero-
Lopez & Moraal, 2012):

$S_{H}=\frac{-(dN/dP)}{J_{H}(P)+1.584F(P)j_{He}}$ (42)

where $F(P)$ is the ratio between the yield functions of $He$ and $H$ reported
by (Clem & Dorman, 2000):

$\ F(P)=F_{0}(P^{a}_{0}+P^{a})^{(\gamma_{1}-\gamma_{2}/a)}P\gamma^{2}$ (43)

the values in the equation are: $F_{0}=2.0$, $\gamma_{1}=0$, $\gamma_{2}=10$,
$a=1.4$ and $P_{0}=0.45$.

The equation 42 allows to compute a modulation factor and will be used to
obtain the semblance. The $J_{H}$ and $J_{He}$ terms in the equation 42 are
the transport equation solution (convection-diffusion or force field) that are
in terms of a LIS previously mentioned. A code on Python has been programmed
and can be found here: https://github.com/JorUrie/TRC.git

##  8 First Analysis

The first analysis consists in three stations located in different latitudes (Table 2) and compute their modulation factor, then semblance between modulation factor with sunspots or mean solar magnetic field will be computed. Data were downloaded from the Neutron Monitor DataBase (NMDB)333NMDB, sunspots data from Sunspot Number444Sunspot Number | SILSO (sidc.be) (Figure 3) and mean solar magnetic field from The Wilcox Solar Observatory555WSO - The Wilcox Solar Observatory (stanford.edu) (Figure 4).

Neutron Monitor | Cut-off Rigidity [GeV]
---|---
SNAE | 0.73
OULU | 0.81
AATB | 5.90
Table 2: Cut-off Rigidity from the stations used in the first analysis Figure
3: Sunspots data obtained from Sunspot index and Long-term Solar Observations.
Figure 4: Solar Magnetic Field date obtained from The Wilcox Solar
Observatory.

The Figure 3 is about the number of sunspots recorded, that are related to
solar maximum and minimum. Data is from May, 1957 to July, 2022. In the same
figure, solar cycles can be observed, in other words, solar maximums belong to
the years 1958, 1970, 1980, 1990, 2014 and the next solar maximum could be in
2024. Solar minimum belongs to the years 1965, 1976, 1986, 1997 and 2010. With
the previous information, 11-year cycle can be perceived. In 1958 was the
larggest recorded maximum, while the smallest maximum was in 1970, as soon as
minimums have similar values and there are not much differences.

Mean magnetic field information belongs to the dates from May 1975 to December
2015. The information has cusps and looks to not have changes, but in
September 1981 was the biggest registered maximum, while in March 1991 was the
smaller registered value. If a line is traced along the graphic following a
”regression”, maximums and minimums could be seen, they would be mostly
inversely proportionate to sunspots.

The selected stations were SNAE (south), OULU (north) and AATB (equator). They
were deemed appropriate based on their respective latitudes in Earth (Figure
24). SNAE, OULU and AATB data are showed in the figures 5, 6 and 7. Data were
downloaded and normalized by contained higher value.

Figure 5: SNAE data obtained from NMDB. Figure 6: OULU data obtained from
NMDB. Figure 7: AATB data obtained from NMDB.

The previous mentioned figures 5, 6 and 7 are normalized data, because the
data is obtained from the neutral monitors on the Earth, however the particles
have already modified by modulation phenomena and they have also lost energy
(secondary cosmic rays), for that reason a modulation factor is necessary to
compute.

Non-normalised data have been introduced into the atmospheric yield function,
which requires a LIS and a transport equation solution (convection-diffusion
or force field) to compute the modulation factor.

The figures from 8 to 15 are the modulation factors for the same station with
every previously mentioned LIS and a transport equation solution. Their
unities are $Counts/s$ and $Date$, because the modulation factor is the amount
of particles before to enter the atmosphere and be influenced by particles or
the magnetosphere of the Earth (primary species) in a given date. Modulation
factor would not be necessary if the satellites would take desired
measurements, but the technology and the conditions do not allow it.

Figure 8: OULU station. Modulation factor with Convection-Diffusion solution
using Lagner, Potgieter and Webber LIS in 2003. Figure 9: OULU station.
Modulation factor with Force Field solution using Lagner, Potgieter and Webber
LIS in 2003. Figure 10: OULU station. Modulation factor with Convection-
Diffusion solution using Burguer and Potgieter LIS in 2000. Figure 11: OULU
station. Modulation factor with Force Field solution using Burguer and
Potgieter LIS in 2000. Figure 12: OULU station. Modulation factor with
Convection-Diffusion solution using Garcia-Munoz, Mason and Simpson LIS in
1975. Figure 13: OULU station. Modulation factor with Force Field solution
using Garcia-Munoz, Mason and Simpson LIS in 1975. Figure 14: OULU station.
Modulation factor with Convection-Diffusion solution using Ghelfi, Barao,
Derome and Maurin LIS in 2017. Figure 15: OULU station. Modulation factor with
Force Field solution using Ghelfi, Barao, Derome and Maurin LIS in 2017.

The previous compute was applied to all selected stations. The spectral
modulation factors have a similar pattern for each LIS and are also similar to
the downloaded data, for example, the figure 6 and figure 8, can be observed
the high values in 1987, near to 2004 and 2011, those occurred for the figures
8 to 15, in fact for the figures 8, 10, 12 and 14 between values $1.2$ and
$1.4$ $Counts/s$, two points are always present and are also the highest
values. For the figures 9, 11, 13 and 15, they have the same pattern and the
two highest points are also present, but between $1.5$ and $1.75$ $Counts/s$,
so the difference of the spectral modulation factors using convection-
diffusion solution and force field solution is the amplitude.

From the figures 8 to 15 also looks like to the figure 3, because there are
more extragalatic cosmic rays introduced into the Earth as a result of the low
solar modulation.

For the three stations, that have been selected according to their latitude,
the semblance between modulation factor and sunspots using several LIS are
showed, however only for OULU station, the rest models can be found in the
appendix.

Figure 16: OULU Station a) Semblance between Modulation Factor using
Convection-Diffusion with Lagner, Potgieter and Webber LIS in 2003 vs Sunspot
data. b) Semblance between Modulation Factor using Force Field with Lagner,
Potgieter and Webber LIS in 2003 vs Sunspot data. c) Difference between the
two previous semblances. Figure 17: OULU Station a) Semblance between
Modulation Factor using Convection-Diffusion with Lagner, Potgieter and Webber
LIS in 2003 vs Mean Magnetic Field data. b) Semblance between Modulation
Factor using Mean Magnetic Field with Lagner, Potgieter and Webber LIS in 2003
vs Sunspot data. c) Difference between the two previous semblances. Figure 18:
OULU Station a) Semblance between Modulation Factor using Convection-Diffusion
with Burguer and Potgieter LIS in 2000 vs Sunspot data. b) Semblance between
Modulation Factor using Force Field with Burguer and Potgieter LIS in 2000 vs
Sunspot data. c) Difference between the two previous semblances. Figure 19:
OULU Station a) Semblance between Modulation Factor using Convection-Diffusion
with Burguer and Potgieter LIS in 2000 vs Mean Magnetic Field data. b)
Semblance between Modulation Factor using Force Field with Burguer and
Potgieter LIS in 2000 vs Mean Magnetic Field data. c) Difference between the
two previous semblances. Figure 20: OULU Station a) Semblance between
Modulation Factor using Convection-Diffusion with Garcia-Munoz, Mason and
Simpson LIS in 1975 vs Sunspot data. b) Semblance between Modulation Factor
using Force Field with Garcia-Munoz, Mason and Simpson LIS in 1975 vs Sunspot
data. c) Difference between the two previous semblances. Figure 21: OULU
Station a) Semblance between Modulation Factor using Convection-Diffusion with
Garcia-Munoz, Mason and Simpson LIS in 1975 vs Mean Magnetic Field data. b)
Semblance between Modulation Factor using Force Field with Garcia-Munoz, Mason
and Simpson LIS in 1975 vs Mean Magnetic Field data. c) Difference between the
two previous semblances. Figure 22: OULU Station a) Semblance between
Modulation Factor using Convection-Diffusion with Ghelfi, Barao, Derome and
Maurin LIS in 2017 vs Sunspots data. b) Semblance between Modulation Factor
using Force Field with Ghelfi, Barao, Derome and Maurin LIS in 2017 vs
Sunspots data. c) Difference between the two previous semblances. Figure 23:
OULU Station a) Semblance between Modulation Factor using Convection-Diffusion
with Ghelfi, Barao, Derome and Maurin LIS in 2017 vs Mean Magnetic Field data.
b) Semblance between Modulation Factor using Force Field with Ghelfi, Barao,
Derome and Maurin LIS in 2017 vs Mean Magnetic Field data. c) Difference
between the two previous semblances.

The computed models offer similar results among them (figures $a)$ and $b)$
for each semblance figure), they even look like if they were the same model,
but thanks to the difference between semblances (figures $c)$) is possible to
observe, that there are little differences mainly in the figures 21 and 20. If
both figures previously mentioned are comparative with the rest of respective
semblances in this section, the main difference is in the interval 50 to 100
in horizontal axis and 40 to 80 in the wavelength, where blue spot change his
size and form.

There is correlation between sunspots and modulation phenomena, which can be
observed in the figures 16, 18, 20 and 22 highlighted in red. There is anti-
correlation between the mean solar magnetic field and the modulation phenomena
observed in blue colour. Both behaviours are maintained independently of the
Earth’s latitude.

##  9 Second Analysis

The same calculus was performed, but the spectral factor modulations were
average while taking into account the LIS. The data from the following
stations were used: KERG, MOSC, OULU, THUL, and HRMS (Table 3).

Neutron Monitor | Cut-off Rigidity [GeV]
---|---
THUL | 0.3
OULU | 0.81
KERG | 1.14
MOSC | 2.43
HRMS | 4.58
Table 3: Cut-off Rigidity from the used stations in the second analysis

The Figure 24 are the position of the used stations in both analyses.

Figure 24: Position of the neutron monitors around the world used in both
analyses

The Figures 25 to 32 are the spectral averaged using convection-diffusion and
force field solutions with Lagner, Potgieter & Webber LIS in 2003, Burguer &
Potgieter LIS in 2000, Garcia-Munoz, Mason & Simpson LIS in 1975 and Ghelfi,
Barao, Derome & Maurin LIS in 2017.

Figure 25: Averaged Modulation Factor among KERG, HRMS, MOSC, OULU and THUL
stations. Modulation factor with Convection-Diffusion solution using Lagner,
Potgieter & Webber LIS in 2003. Figure 26: Averaged Modulation Factor among
KERG, HRMS, MOSC, OULU and THUL stations. Modulation factor with Force Field
solution using Lagner, Potgieter & Webber LIS in 2003.. Figure 27: Averaged
Modulation Factor among KERG, HRMS, MOSC, OULU and THUL stations. Modulation
factor with Convection-Diffusion solution using Burguer and Potgieter LIS in
2000. Figure 28: Averaged Modulation Factor among KERG, HRMS, MOSC, OULU and
THUL stations. Modulation factor with Force Field solution using Burguer and
Potgieter LIS in 2000. Figure 29: Averaged Modulation Factor among KERG, HRMS,
MOSC, OULU and THUL stations. Modulation factor with Convection-Diffusion
solution using Garcia-Munoz, Mason and Simpson LIS in 1975. Figure 30:
Averaged Modulation Factor among KERG, HRMS, MOSC, OULU and THUL stations.
Modulation factor with Force Field solution using Garcia-Munoz, Mason and
Simpson LIS in 1975. Figure 31: Averaged Modulation Factor among KERG, HRMS,
MOSC, OULU and THUL stations. Modulation factor with Convection-Diffusion
solution using Ghelfi, Barao, Derome and Maurin LIS in 2017. Figure 32:
Averaged Modulation Factor among KERG, HRMS, MOSC, OULU and THUL stations.
Modulation factor with Force Field solution using Ghelfi, Barao, Derome and
Maurin LIS in 2017.

The Figures 33 to 48 are the semblance between the aforementioned modulation
factors with sunspots and mean magnetic field data. The figures 34, 36, 38, 40
and 42 are the percentage of contained zeros into the difference between
semblance (Figures $c)$), that were calculated to observed with more detail
which LIS has the best adjusted in order to find the most reliable model.

Figure 33: Average among HRMS, KERG, MOSC, OULU and THUL a) Semblance between
Average Modulation Factor using Convection-Diffusion with Lagner, Potgieter
and Webber LIS in 2003 vs Sunspots data. b) Semblance between Average
Modulation Factor using Force Field with Lagner, Potgieter and Webber LIS in
2003 vs Sunstpots data. c) Difference between the two previous semblances.
Figure 34: Contained zeros into Figure 33 c). Figure 35: Average among HRMS,
KERG, MOSC, OULU and THUL a) Semblance between Average Modulation Factor using
Convection-Diffusion with Lagner, Potgieter and Webber LIS in 2003 vs Solar
Magnetic Field data. b) Semblance between Average Modulation Factor using
Force Field with Lagner, Potgieter and Webber LIS in 2003 vs Solar Magnetic
Field data. c) Difference between the two previous semblances. Figure 36:
Contained zeros into Figure 35 c). Figure 37: Average among HRMS, KERG, MOSC,
OULU and THUL a) Semblance between Average Modulation Factor using Convection-
Diffusion with Burguer and Potgieter LIS in 2000 vs Sunspots data. b)
Semblance between Average Modulation Factor using Force Field with Burguer and
Potgieter LIS in 2000 vs Sunstpots data. c) Difference between the two
previous semblances. Figure 38: Contained zeros into Figure 37 c). Figure 39:
Average among HRMS, KERG, MOSC, OULU and THUL a) Semblance between Average
Modulation Factor using Convection-Diffusion with Burguer and Potgieter LIS in
2000 vs Solar Magnetic Field data. b) Semblance between Average Modulation
Factor using Force Field with Burguer and Potgieter LIS in 2000 vs Solar
Magnetic Field data. c) Difference between the two previous semblances. Figure
40: Contained zeros into Figure 39 c). Figure 41: Average among HRMS, KERG,
MOSC, OULU and THUL a) Semblance between Average Modulation Factor using
Convection-Diffusion with Garcia-Munoz, Mason and Simpson LIS in 1975 vs
Sunspots data. b) Semblance between Average Modulation Factor using Force
Field with Garcia-Munoz, Mason and Simpson LIS in 1975 vs Sunstpots data. c)
Difference between the two previous semblances. Figure 42: Contained zeros
into Figure 41 c). Figure 43: Average among HRMS, KERG, MOSC, OULU and THUL a)
Semblance between Average Modulation Factor using Convection-Diffusion with
Garcia-Munoz, Mason and Simpson LIS in 1975 vs Solar Magnetic Field data. b)
Semblance between Average Modulation Factor using Force Field with Garcia-
Munoz, Mason and Simpson LIS in 1975 vs Solar Magnetic Field data. c)
Difference between the two previous semblances. Figure 44: Contained zeros
into Figure 43 c). Figure 45: Average among HRMS, KERG, MOSC, OULU and THUL a)
Semblance between Average Modulation Factor using Convection-Diffusion with
Ghelfi, Barao, Derome and Maurin LIS in 2017 vs Sunspots data. b) Semblance
between Average Modulation Factor using Force Field with Ghelfi, Barao, Derome
and Maurin LIS in 2017 vs Sunstpots data. c) Difference between the two
previous semblances. Figure 46: Contained zeros into Figure 45 c). Figure 47:
Average among HRMS, KERG, MOSC, OULU and THUL a) Semblance between Average
Modulation Factor using Convection-Diffusion with Ghelfi, Barao, Derome and
Maurin LIS in 2017 vs Solar Magnetic Field data. b) Semblance between Average
Modulation Factor using Force Field with Ghelfi, Barao, Derome and Maurin LIS
in 2017 vs Solar Magnetic Field data. c) Difference between the two previous
semblances. Figure 48: Contained zeros into Figure 47 c).

In this case, an average modulation factor of several neutron detectors around
the word was computed to analyse the general behaviour between sunspots and
mean solar magnetic field with modulation phenomena. The results were utterly
similar to the first analysis, in fact the analysis of the contained zeros
confirms the hypothesis, that there is correlation between sunspots and
modulation phenomena and anti-correlation between the mean magnetic field and
modulation phenomena. The analysis of the zeros can be summarised in the Table
4.

Analysis of contained zeros | Sunspots | Mean Magnetic Field
---|---|---
|  <0 |  >0 |  <0 |  > 0
Webber & Lockwood | 47 | 53 | 46 | 54
Burguer | 47 | 53 | 46 | 54
Garcia-Munoz | 49 | 51 | 47 | 53
Maurin | 49 | 51 | 46 | 54
Table 4: Analysis of contained zeros from difference between semblance using
transport equation solution convection-diffusion and force field and different
LIS

##  10 Conclusions

Both analyses allow quantifying the correlation between modulation factor with
sunspots and anti-correlation between modulation factor with mean magnetic
field. The models are mostly similar with both transport equation solutions
and the local interstellar spectrum, because just shown light changes.

For the first analysis, the modulation factors yield similar outcomes and are
solely dependent on the rigidity factor. The correlation between sunspots and
modulation factor is more pronounced for stations closer to the poles, such as
OULU and SNAE, than the equatorial AATB station. However, this pattern is not
replicated for the semblance between the modulation factor and the mean
magnetic field. The station near the equator shown more anti-correlation than
near the poles. The pattern is similar for every semblances, which were used
different LIS and transport equation solution.

Thanks to the second analysis, the first analysis could be affirmed and
generalized, the force field solution is also more compatible with neutron
detectors data, Additionally any modulation factor solution will be similar to
the Neutron detector data and just will have minimal amplitude variations.

In the second analysis, it is possible to observe similarity among the graphs
of contained zeros into the difference of semblances, particularly when the
modulation factors have been computed by different LIS. However if semblance
is computed using the Garcia-Munoz, Mason and Simpson LIS in 1975 shown more
differences, that is visible in the analysis of contained zeros for semblance
between averaged modulation factor and mean magnetic field.

According to the table 4 semblance between sunspots and averaged modulation
factors are really similar and the percentage of contained zeros are akin if
Webber & Lockwood and Burguer and Potgieter LIS in 2000 are used, same
happened for Garcia-Munoz and Ghelfi, Barao, Derome and Maurin LIS in 2017.

Therefore, the LIS that have the best adjust to the semblance model are
Lagner, Potgieter and Webber LIS in 2003, Burguer and Ghelfi, Barao, Derome
and Maurin LIS in 2017 using force field solution because they offer more
consistent values as shown in the table 4.

With both proposed models can be asserted with greater certainty that exist a
correlation between modulation phenomena and sunspots, because the solar
activity is less allowed to detect more rays from the Sun. The anti-
correlation between modulation phenomena and mean magnetic field is due to the
high solar activity, therefore there are less solar rays that affect the
arisen of cosmic rays from others source.

Thanks to this proposed quantitative model, a more precise predictive model of
cosmic rays with solar activity can be compute.

## References

  * Auger et al. (1939) Auger, P., Ehrenfest, P., Maze, R., Daudin, J., & Fréon, R. A. 1939, Reviews of Modern Physics, 11, 288, doi: 10.1103/RevModPhys.11.288
  * Bieber (1999) Bieber, J. 1999, in International Cosmic Ray Conference, Vol. 7, 26th International Cosmic Ray Conference (ICRC26), Volume 7, 61
  * Boschini et al. (2017) Boschini, M. J., Della Torre, S., Gervasi, M., et al. 2017, ApJ, 840, 115, doi: 10.3847/1538-4357/aa6e4f
  * Burger et al. (2000) Burger, R. A., Potgieter, M. S., & Heber, B. 2000, Journal of Geophysical Research: Space Physics, 105, 27447, doi: https://doi.org/10.1029/2000JA000153
  * Caballero-Lopez & Moraal (2004) Caballero-Lopez, R. A., & Moraal, H. 2004, Journal of Geophysical Research: Space Physics, 109, doi: https://doi.org/10.1029/2003JA010098
  * Caballero-Lopez & Moraal (2012) —. 2012, Journal of Geophysical Research: Space Physics, 117, doi: https://doi.org/10.1029/2012JA017794
  * Clem & Dorman (2000) Clem, J., & Dorman, L. 2000, Space Science Reviews, 93, doi: 10.1023/A:1026508915269
  * Cooper & Cowan (2008) Cooper, G. R. J., & Cowan, D. R. 2008, Computers and Geosciences, 34, 95, doi: 10.1016/j.cageo.2007.03.009
  * Crank & Nicolson (1947) Crank, J., & Nicolson, P. 1947, Mathematical Proceedings of the Cambridge Philosophical Society, 43, 50–67, doi: 10.1017/S0305004100023197
  * Ferrari et al. (1996) Ferrari, A., Ranft, J., Roesler, S., & Sala, P. R. 1996, Zeitschrift für Physik C: Particles and Fields, 71, 75, doi: 10.1007/s002880050149
  * Fisk et al. (1974) Fisk, L. A., Goldstein, M. L., Klimas, A. J., & Sandri, G. 1974, ApJ, 190, 417, doi: 10.1086/152893
  * Garcia-Munoz et al. (1975) Garcia-Munoz, M., Mason, G. M., & Simpson, J. A. 1975, ApJ, 202, 265, doi: 10.1086/153973
  * Ghelfi et al. (2017) Ghelfi, A., Barao, F., Derome, L., & Maurin, D. 2017, Astronomy & Astrophysics, 605, C2, doi: 10.1051/0004-6361/201527852e
  * Gleeson & Axford (1968) Gleeson, L. J., & Axford, W. I. 1968, Canadian Journal of Physics Supplement, 46, 937, doi: 10.1139/p68-388
  * James (1998) James, F. a. 1998.  https://cds.cern.ch/record/2296388
  * Langner et al. (2003) Langner, U. W., Potgieter, M. S., & Webber, W. R. 2003, Journal of Geophysical Research: Space Physics, 108, doi: https://doi.org/10.1029/2003JA009934
  * Luo et al. (2015) Luo, X., Zhang, M., Potgieter, M., Feng, X., & Pogorelov, N. V. 2015, ApJ, 808, 82, doi: 10.1088/0004-637X/808/1/82
  * Moskalenko et al. (2002) Moskalenko, I. V., Strong, A. W., Ormes, J. F., & Potgieter, M. S. 2002, The Astrophysical Journal, 565, 280, doi: 10.1086/324402
  * Neidell & Taner (1971) Neidell, N. S., & Taner, M. T. 1971, GEOPHYSICS, 36, 482, doi: 10.1190/1.1440186
  * O’Neill (2006) O’Neill, P. M. 2006, Advances in Space Research, 37, 1727, doi: 10.1016/j.asr.2005.02.001
  * Parker (1965) Parker, E. N. 1965, Planet. Space Sci., 13, 9, doi: 10.1016/0032-0633(65)90131-5
  * Potgieter et al. (2012) Potgieter, M., De Simone, N., Vos, E., Boezio, M., & Di Felice, V. 2012, in 39th COSPAR Scientific Assembly, Vol. 39, 1526
  * Putze & Derome (2014) Putze, A., & Derome, L. 2014, Physics of the Dark Universe, 5, 29, doi: 10.1016/j.dark.2014.07.002
  * Sartini et al. (2010) Sartini, L., Simeone, F., Pani, P., et al. 2010, Nuclear Instruments and Methods in Physics Research A
  * Sekido & Elliot (1985) Sekido, Y., & Elliot, H. 1985, Early history of cosmic ray studies: Personal reminiscences with old photographs.  https://www.osti.gov/biblio/7169147
  * Shikaze et al. (2007) Shikaze, Y., Haino, S., Abe, K., et al. 2007, Astroparticle Physics, 28, 154, doi: 10.1016/j.astropartphys.2007.05.001
  * Strong et al. (2007) Strong, A. W., Moskalenko, I. V., & Ptuskin, V. S. 2007, Annual Review of Nuclear and Particle Science, 57, 285, doi: 10.1146/annurev.nucl.57.090506.123011
  * Tan & Ng (1983) Tan, L. C., & Ng, L. K. 1983, ApJ, 269, 751, doi: 10.1086/161084
  * Tan & Ng (1983) Tan, L. C., & Ng, L. K. 1983, Journal of Physics G: Nuclear Physics, 9, 227, doi: 10.1088/0305-4616/9/2/015
  * Taner & Koehler (1969) Taner, M. T., & Koehler, F. 1969, GEOPHYSICS, 34, 859, doi: 10.1190/1.1440058
  * Teolis (2017) Teolis, A. 2017, Computational Signal Processing with Wavelets, doi: 10.1007/978-3-319-65747-9
  * Vos & Potgieter (2015) Vos, E. E., & Potgieter, M. S. 2015, ApJ, 815, 119, doi: 10.1088/0004-637X/815/2/119

##  11 Appendix

Figure 49: HRMS station. Modulation factor with Convection-Diffusion using
Lagner, Potgieter & Webber LIS in 2003

Figure 50: HRMS station. Modulation factor with Force Field using Lagner,
Potgieter & Webber LIS in 2003

Figure 51: HRMS station. Modulation factor with Convection-Diffusion using
Burguer and Potgieter LIS in 2000

Figure 52: HRMS station. Modulation factor with Force Field using Burguer and
Potgieter LIS in 2000

Figure 53: HRMS station. Modulation factor with Convection-Diffusion using
Garcia-Munoz, Mason and Simpson LIS in 1975

Figure 54: HRMS station. Modulation factor with Force Field using Garcia-
Munoz, Mason and Simpson LIS in 1975

Figure 55: HRMS station. Modulation factor with Convection-Difussion using
Ghelfi, Barao, Derome and Maurin LIS in 2017 in 2017

Figure 56: HRMS station. Modulation factor with Force Field using Ghelfi,
Barao, Derome and Maurin LIS in 2017 in 2017

Figure 57: SNAE Station a) Semblance between Modulation Factor using
Convection-Diffusion with Lagner, Potgieter and Webber LIS in 2003 vs Sunspot
data. b) Semblance between Modulation Factor using Force Field with Lagner,
Potgieter and Webber LIS in 2003 vs Sunspot data. c) Difference between the
two previous semblances. Figure 58: SNAE Station a) Semblance between
Modulation Factor using Convection-Diffusion with Lagner, Potgieter and Webber
LIS in 2003 vs Mean Magnetic Field data. b) Semblance between Modulation
Factor using Force Field with Lagner, Potgieter and Webber LIS in 2003 vs Mean
Magnetic Field data. c) Difference between the two previous semblances. Figure
59: SNAE Station a) Semblance between Modulation Factor using Convection-
Diffusion with Burguer and Potgieter LIS in 2000 vs Sunspot data. b) Semblance
between Modulation Factor using Force Field with Burguer and Potgieter LIS in
2000 vs Sunspot data. c) Difference between the two previous semblances.
Figure 60: SNAE Station a) Semblance between Modulation Factor using
Convection-Diffusion with Burguer and Potgieter LIS in 2000 vs Mean Magnetic
Field data. b) Semblance between Modulation Factor using Force Field with
Burguer and Potgieter LIS in 2000 vs Mean Magnetic Field data. c) Difference
between the two previous semblances. Figure 61: SNAE Station a) Semblance
between Modulation Factor using Convection-Diffusion with Garcia-Munoz, Mason
and Simpson LIS in 1975 vs Sunspot data. b) Semblance between Modulation
Factor using Force Field with Garcia-Munoz, Mason and Simpson LIS in 1975 vs
Sunspot data. c) Difference between the two previous semblances. Figure 62:
SNAE Station a) Semblance between Modulation Factor using Convection-Diffusion
with Garcia-Munoz, Mason and Simpson LIS in 1975 vs Mean Magnetic Field data.
b) Semblance between Modulation Factor using Force Field with Garcia-Munoz,
Mason and Simpson LIS in 1975 vs Mean Magnetic Field data. c) Difference
between the two previous semblances. Figure 63: SNAE Station a) Semblance
between Modulation Factor using Convection-Diffusion with Ghelfi, Barao,
Derome and Maurin LIS in 2017 in 2017 vs Sunspots data. b) Semblance between
Modulation Factor using Force Field with Ghelfi, Barao, Derome and Maurin LIS
in 2017 in 2017 vs Sunspots data. c) Difference between the two previous
semblances. Figure 64: SNAE Station a) Semblance between Modulation Factor
using Convection-Diffusion with Ghelfi, Barao, Derome and Maurin LIS in 2017
in 2017 vs Mean Magnetic Field data. b) Semblance between Modulation Factor
using Force Field with Ghelfi, Barao, Derome and Maurin LIS in 2017 vs Mean
Magnetic Field data. c) Difference between the two previous semblances.