Predicting the turbulent transport of cosmic rays via neural networks
A fast artificial neural network is developed for the prediction of cosmic
ray transport in turbulent astrophysical magnetic fields. The setup is trained
and tested on bespoke datasets that are constructed with the aid of
test-particle numerical simulations of relativistic cosmic ray dynamics in
synthetic stochastic fields. The neural network uses, as input, particle and
field properties and estimates transport coefficients 10^7 faster than standard
numerical simulations with an overall error of ~5% .
D. I. Palade
# Predicting the turbulent transport of cosmic rays via neural networks

###### Abstract

A fast artificial neural network is developed for the prediction of cosmic ray transport in turbulent astrophysical magnetic fields. The setup is trained and tested on bespoke datasets that are constructed with the aid of test-particle numerical simulations of relativistic cosmic ray dynamics in synthetic stochastic fields. The neural network uses, as input, particle and field properties and estimates transport coefficients \(10^{7}\) faster than standard numerical simulations with an overall error of \(\sim 5\%\).

cosmic ray, magnetic turbulence, turbulent transport, neural networks

## I Introduction

Cosmic rays (CRs) are charged particles originating mostly in galactic supernovae remnants and being accelerated by the Fermi mechanism in shock waves [1; 2] although, more exotic sources are possible [3]. Due to their (ultra)relativistic energies, CRs can permeate systems of all sizes, from the heliosphere, to galaxies and even the extra-galactic medium [3; 4]. Consequently, CRs carry important information to our understanding of fundamental high-energy physics, astrophysical magnetic fields, the structure of astrophysical media, space weather, etc. [5; 6; 7].

The dynamics of CRs is a long-standing fundamental problem of astrophysics with many open questions that are difficult to solve for several reasons. During their journey from sources to detection, CRs interact mainly with magnetic fields. The latter are highly turbulent [8] and, via non-linear interactions with particles, lead to consistent (anomalous) transport phenomena beyond ballistic motion [9]. The CR's energies cover a wide range of magnitudes (from \(MeV\) to \(10^{20}eV\)[2]) while the turbulent features of the magnetic fields are quite diverse, spanning different space-scales, anisotropies, and fluctuation amplitudes [10]. This richness in physical regimes opens a variety of possible transport types ranging from subdiffusive to superdiffusive in the perpendicular and parallel directions, with complicated dependencies [9; 10; 11; 12].

Despite this discouraging picture, a lot of progress has been made in the past decades in the topic of turbulent CRs transport. Quasilinear approaches [13] and nonlinear extensions [14] gained a lot of attention due to their technical simplicity and ability to provide insight in the physical mechanisms involved. Unfortunately, such models are known to fail in a lot of relevant cases such as the limit of strong turbulence or the 90 _degree problem_. A much more accurate approach, which gained momentum in the scientific community, is the method of test-particle simulations either in synthetic [15; 16; 17; 18] or MHD generated [19; 20] turbulent magnetic fields. Within this approach, the dynamics of CRs is mimicked at the statistical level with an ensemble of fictitious particles that allows for a direct evaluation of the transport coefficients. Its main drawback is the relatively high numerical effort required, which is hardly compatible with the diversity of possible physical regimes. It is clear that the astrophysics community would benefit from a fast and accurate method of prediction in their quest to understand the CR turbulent transport.

A helping hand might come from another scientific front. In the recent years, we have witnessed the rise of artificial intelligence (AI) methods, in particular machine learning (ML) techniques [21; 22; 23], that are able to provide inferences in various mathematical problems. From the multitude of ML methods, the reader is directed towards artificial neural networks (ANNs) [22; 23]. The promise of ML is that, regardless of the chosen technique, if sufficient data is available, it can learn from that data and make fast and reliable prognoses for unknown cases. Such abilities would be equivalent to having analytical expressions at hand and would bypass the technical difficulties that arise in doing many simulations of CR dynamics.

Given this scientific context, the purpose of the present work is to illustrate the methodology for developing an artificial neural network designed for predictions of cosmic ray turbulent transport in astrophysical magnetic fields. For the training and testing phases, a database is constructed with the aid of test-particle numerical simulations in synthetically generated fields. Since the main purpose is methodological, we restrict ourselves here to the evaluation of perpendicular diffusion [9] in a relatively limited range of physical parameters. Nonetheless, it is the hope of the author that this work will open a path towards more extensive databases and, consequently, more potent ANNs.

The rest of this manuscript is organized as follows. The Theory section II describes, briefly, the CR transport model (II.1), the test-particle numerical method (II.2), the general architecture of ANNs (II.3) and the construction of a training and testing database (II.4). The Results Section (III) presents the convergence properties and the predictive power of the ANN. In the Con
clusions section (IV), the main findings are resumed and perspectives are discussed.

## II Theory

The ingredients of an ANN are the programming structure and the training dataset. For our problem, the latter is constituted of input-output pairs representing field-particle properties and diffusion coefficients. Such data is obtained with the numerical method of test particles that move according to a transport model. In this section, these elements, shown schematically in Fig. (1), are discussed in reverse order.

### The transport model

A cosmic ray is a charged particle characterized by its rest mass \(m\), charge \(q\), position \(\mathbf{r}(t)\) and velocity \(\mathbf{v}(t)=d\mathbf{r}(t)/dt\) in a global reference frame obeying the relativistic Newton-Lorentz equation in the presence of an astrophysical magnetic field \(\mathbf{B}\):

\[\frac{d(\gamma\mathbf{v})}{dt}=\frac{q}{m}\left(\mathbf{v}\times\mathbf{B} \right). \tag{1}\]

The kinetic energy of the particle is \(T=mc^{2}(\gamma-1)\) where the Lorentz factor \(\gamma=(1-v^{2}/c^{2})^{-1/2}\), \(c\) is the speed of light and any electric field is neglected \(\mathbf{E}=0\). A cartesian system of coordinates \(\mathbf{r}\equiv(x,y,z)\) is defined. The magnetic field \(\mathbf{B}\) is decomposable into an average constant component \(B_{0}\) along \(Oz\) and a fluctuating part \(\mathbf{b}\) which is turbulent. Relative to this decomposition we coin the \(Oz\) direction as being "parallel" \(\hat{e}_{z}\equiv\hat{e}_{\parallel}\) and the plane \((x,y)\) as "perpendicular". Naturally, for any wavevector \(\mathbf{k}=(k_{x},k_{y},k_{z})\), \(k_{z}\equiv k_{\parallel}\) and \((k_{x},k_{y})\equiv\mathbf{k}_{\perp}\). The velocity \(\mathbf{v}=v_{\parallel}\hat{e}_{\parallel}+\mathbf{v}_{\perp}\) can be used to define the pitch angle \(\mu=v_{\parallel}/v\), where \(v^{2}=|\mathbf{v}|^{2}\).

The turbulent component \(\mathbf{b}\) of the total magnetostatic field \(\mathbf{B}=B_{0}\hat{e}_{\parallel}+\mathbf{b}\) is represented in the paradigm of a 2D model [16; 18; 24]:

\[\mathbf{b}(\mathbf{r})=\hat{e}_{\parallel}\times\nabla_{\perp}a_{\parallel} (\mathbf{r}). \tag{2}\]

The magnetic vector potential \(a_{\parallel}\) is Gaussian, zero-averaged \(\langle a_{\parallel}(\mathbf{r})\rangle=0\) and homogeneous [15; 17; 24]. The last property is evident in the spectrum of its fluctuations:

\[\langle\tilde{a}_{\parallel}(\mathbf{k})\tilde{a}_{\parallel}(\mathbf{k}^{ \prime})\rangle=S(\mathbf{k})\delta(\mathbf{k}+\mathbf{k}^{\prime}). \tag{3}\]

By \(\tilde{a}_{\parallel}(\mathbf{k})\) was denoted the Fourier transform of \(a_{\parallel}(\mathbf{r})\), while \(\langle\cdot\rangle\) stands for statistical average. The overall turbulence amplitude is defined as \(b=\sqrt{(\mathbf{b}^{2}(\mathbf{r}))}\) and the spectrum is assumed to be of Kolmogorov type with parallel-perpendicular anisotropic dependencies [15; 18; 24]:

\[S(\mathbf{k}_{\perp},k_{\parallel})\sim\frac{(k_{\perp}\lambda_{\perp})^{q} \left(1+(k_{\parallel}\lambda_{\parallel})^{2}\right)^{-s/2}}{k_{\perp}(1+(k_ {\perp}\lambda_{\perp})^{2})^{(s+q)/2}}, \tag{4}\]

where \(q=3,s=5/3\) and \(\lambda_{\parallel},\lambda_{\perp}\) are "bend-over" scales, intimately related to coherence/integral scales [15].

Within this transport model, realizations of the turbulent field \(\mathbf{b}\) drive associated CR trajectories \(\mathbf{r}(t)\) via the eq. of motion (1). The ensemble of trajectories can be used to derive, as statistical averages, the diffusion transport coefficients in any direction \(Ox\):

\[D_{x,x}(t)=\frac{\langle\left(x(t)-x(0)\right)^{2}\rangle}{2t}. \tag{5}\]

In this paper, we are discussing regimes that are strictly diffusive and for which the quantity of interest is the asymptotic perpendicular diffusion coefficient \(D_{\perp}=\lim_{t\rightarrow\infty}D_{x,x}(t)\). In this context, the problem of CR turbulent transport is equivalent to asking: _what are the diffusion coefficients \(Y=D_{\perp}\) for any given set of particle-field input parameters \(X=(T,\mu,b,\lambda_{\perp},\lambda_{\parallel})\)_?

For practical reasons related to numerical implementation, all quantities involved in the transport model (eqns. (1)-(4)) are scaled as follows. The magnetic fields \(\mathbf{b}\rightarrow\mathbf{b}B_{0}\), the time \(t\to\omega_{c}^{-1}\), the velocities \(\mathbf{v}\rightarrow\mathbf{v}\mathbf{c}\), the kinetic energy \(T\to Tmc^{2}\), the space scales \((\mathbf{r},\lambda_{i})\rightarrow(\mathbf{r},\lambda_{i})\rho_{L}^{0}\) and the wave-vectors \(\mathbf{k}\rightarrow\mathbf{k}/\rho_{L}^{0}\). The cyclotron frequency is defined as \(\omega_{c}=|q|B_{0}/m\) while the "bare" Larmor radius \(\rho_{L}^{0}=mc/(|q|B_{0})\). Consequently, the diffusion coefficients are scaled as \(D\to Dm_{0}c^{2}/(qB_{0})\).

### Test-particle numerical simulations

The present work employs the method of test-particle simulations [15; 17] (or direct numerical simulations [18; 25; 26]) for the calculus of diffusion coefficients. The main idea is that we can mimic the reality of turbulent dynamics through a numerical representation of

Figure 1: Schematic view of the elements required to construct the ANN.
the transport model described in the previous subsection (II.1). A finite ensemble of \(N_{p}\) random fields \(\{\mathbf{b}\}\) with appropriate spectral properties (distribution and spectra) is constructed. For each realization of this ensemble, a CR trajectory \(\mathbf{r}(t)\) can be computed by solving eq. (1). In the limit of many particles \(N_{p}\rightarrow\infty\), the transport coefficients are given by statistical averages over the numerical ensemble of trajectories \(\{\mathbf{r}(t)\}\) (5). More details on the method can be found in [26].

The magnetic vector potential, \(a_{\parallel}\), is assumed Gaussian, zero-averaged, and homogeneous; thus, it can be constructed numerically in each realization with a Fourier-like (harmonic) decomposition [27; 26; 17]:

\[a_{\parallel}(\mathbf{r})=b\lambda_{\perp}\sqrt{\frac{2}{N_{c}}}\sum_{j=1}^{N _{c}}\sin(\mathbf{k}^{j}\mathbf{r}+\alpha^{j}) \tag{6}\]

where \(\alpha^{j}\) are independent random initial phases distributed uniformly in the interval \([0,2\pi)\). The partial waves \(\mathbf{k}^{j}\equiv\left(\mathbf{k}^{j}_{\perp},k^{j}_{\parallel}\right)\) are independent random variables generated with the use of the acceptance-rejection algorithm [28] for the PDF \(S(\mathbf{k})\). It can be shown that the representation (6) is, indeed, zero-averaged with the correct spectrum (4). In the limit of many partial waves \(N_{c}\rightarrow\infty\), the Gaussianity is also achieved (provided by the _central limit theorem_).

Solving eq. (1) is a standard numerical problem which is tackled here with a relativistic Boris pusher [29]. Magnetic field values are evaluated in each realization directly using eqn. (6). Within the present scaling, the gyro-period is \(2\pi\) and the time step is chosen \(\delta t=2\pi/10\). This level of time discretization was found to be sufficient for an accurate prediction of the Larmor rotation and CR dynamics.

Practical experience has shown that using \(N_{p}=10^{4}\) and \(N_{c}=500\) in a simulation is enough to achieve acceptable Gaussianity, Eulerian, and Lagrangian convergence, thus reliable results on transport.

### Artificial neural networks

In this section, the generic features of an ANN architecture that will be used for predicting the transport of CRs in turbulent magnetic media are described. Its structure is represented graphically in Fig. (2). For more details, see [30].

A general ANN is a computing setup consisting of \(m+2\) "layers", each layer \(i\) contains \(n_{i}\) "neurons" and each neuron \(j\) is described by a numerical value \(Z_{i,j}\), a "bias" \(\beta_{i,j}\), an "activation function" \(f_{i,j}\) and some "weights" \(W_{i,j,k}\). The latter connects recursively neurons from neighbouring layers, thus, \(W\) is an irregular tensor of \((m+1)\times n_{i}\times n_{i-1}\) dimension. The first layer, \(i=0\), contains the values of the input variable \(X\), \(Z_{0,j}=X(j)\). The last layer, \(i=m+1\) contains the values of the output variable \(Y\), \(Z_{m+1,j}=Y(j)\). The rest are coined "hidden layers". Within an ANN there is a feed-forward recurrent relation between layers:

\[Z_{i,j}=f_{i,j}\left(\sum_{k=1}^{n_{i-1}}W_{i,j,k}.Z_{i-1,k}+\beta_{i,j}\right) \tag{7}\]

The purpose of an ANN is to model the true relation (which is unknown) between the input variable \(X\) and the output \(Y\). This is achieved by finding appropriate weights and biases \(\{W_{i,j,k},\beta_{i,j}\}\) that minimize an error function \(E\) between the output of the network, \(Z_{m+1}\), and the real output, \(Y\). This must be done over a dataset of \(N_{d}\) pairs \(\{X^{p},Y^{p}\}_{p=1,N_{d}}\) which should be _reasonably large_. The error function \(E\) is defined as Minkowski error:

\[E=\frac{1}{n}\sum_{p=1}^{N_{d}}\sum_{j=1}^{n_{m+1}}\left(Z_{m+1,j}^{p}-Y^{p}(j )\right)^{2}. \tag{8}\]

The error minimization is reached iteratively, over multiple "epochs", using ADAM [31], a modified version of gradient descent algorithm [32]. In this approach, during each iteration (epoch), the weights and biases \(\{W_{i,j,k},\beta_{i,j}\}\) (denoted generically \(\theta\)), are updated using a "learning"/"training" procedure:

\[v_{\theta} =\alpha_{1}v_{\theta}+(1-\alpha_{1})\left(\frac{\partial E}{ \partial\theta}\right)^{2} \tag{9}\] \[m_{\theta} =\frac{\alpha_{2}m_{\theta}}{1-\alpha_{1}^{epoch}}+\frac{1- \alpha_{2}}{1-\alpha_{1}^{epoch}}\left(\frac{\partial E}{\partial\theta}\right)\] (10) \[\theta =\theta-\gamma\frac{m_{\theta}}{\varepsilon+\sqrt{v_{\theta}}} \tag{11}\]

where \(\gamma,\alpha_{1},\alpha_{2},\varepsilon\) are parameters controlling the learning rate and the convergence of the ANN. The "moments" \(v_{\theta},m_{\theta}\) are initialized to zero, while the weights and biases \(\theta\equiv\{W_{i,j,k},\beta_{i,j}\}\) are normal random variables with zero average and unit variance. The ADAM algorithm (9)-(11) should find, asymptotically, the global minimum of

Figure 2: The architecture of a general ANN with 5 input and 2 output variables.
the error function \(E\). This means that the ANN should predict output values \(Z_{m+1}^{p}\) as suitable approximations for the real \(Y^{p},\forall p=\overline{1,N_{d}}\) and beyond.

Based on numerical experience, the optimal values used for the present ANN are \(m=2\) hidden layers and \(n_{0}=3,n_{1}=10,n_{2}=5,n_{3}=1\) neurons in each layer. The activation functions are \(f_{i,j}(z)=\tanh(z),\forall i=\overline{1,m},j=\overline{1,n_{i}}\) and \(f_{m+1,j}(z)=z^{2},\forall j=\overline{1,n_{m+1}}\) for the hidden, respectively, output layers. The parameters for ADAM are \(\gamma=0.001,\alpha_{1}=0.9,\alpha_{2}=0.95,\varepsilon=10^{-7}\).

### Defining the database

At this point, we must recall how the question of CR turbulent transport fits into the formal description of an ANN. The eqns. of the model (1)-(5) have a set of free parameters, namely \(\left(T,\mu,b,\lambda_{\perp},\lambda_{\parallel}\right)\). In our simulations, we chose \(\mu=0.5\) and \(\lambda_{\parallel}=10\lambda_{\perp}\), values relevant to observations from solar winds. The rest, \(X=(T,b,\lambda_{\perp})\) represents the input variable from ANN's perspective. For each set of parametric values \(X\), associated perpendicular diffusion coefficients are obtained as output \(Y=D_{\perp}\). In order to train and test the ANN, we need to construct a database of \(N_{d}\) pairs \((X^{p},Y^{p})_{p=1,N_{d}}\).

The latter is obtained through the test-particle method described in Section (II.2). An appropriate number of \(N_{d}=5\times 10^{3}\) simulations has been performed. Note that this task is far more demanding than constructing and training the ANN, since it requires a much larger volume of CPU computing time.

The \(n_{0}-\)dimensional space of \(X\) parameters must be truncated to a finite domain. This is done by choosing a cube of limiting values: \(T\in(0,10)\) which for protons corresponds to a kinetic energy \((0eV,10GeV)\), \(b\in(0,1)\) corresponding to \((0,B_{0})\) amplitude of fluctuations, \(\lambda_{\perp}\in(0.2,20)\). Note that for a proton, \(T=1\) corresponds to \(\approx 1GeV\) energy and a scaled Larmor radius \(\rho_{L}\approx 1.7\), thus, comparable with most values of correlation lengths, \(\lambda_{\perp}\).

In order to avoid spurious biases, the values \(\{X^{p}\}\) are generated randomly uniform inside the parametric hypercube. It is important to observe that the database can be completed with analytical results. For example, at \(b=0\to D_{\perp}=D_{\parallel}=0\). Approximately \(5-10\%\) of the database can be filled with such exact results. If this fraction is too large, the hypercube becomes non-uniformly filled, and the ANN learns predominately the analytical values, which are of no interest.

From the total database, \(90\%\) of the \((X,Y)\) pairs are used for training, while the remaining \(10\%\) will serve as testing grounds for the assessment of ANN's predictive abilities.

## III Results

### Turbulent transport

The transport model has been used in the framework of the test-particle method to evaluate CR trajectories and, consequently, perpendicular diffusion coefficients \(Y=D_{\perp}\) for different random combinations of free input parameters \(X\equiv(T,b,\lambda_{\perp})\) values. The equations of motion (1) are fully relativistic and do not rely on the guiding-center approximation [18]. Consequently, particle trajectories describe the full Larmor rotation induced by the average magnetic field \(B_{0}\) as well as the "scattering" in the fluctuating field \(\mathbf{b}\) which is 2D. This level of description is needed whenever the characteristic length scales of fluctuations, in our case \(\lambda_{\perp}\), are comparable with the Larmor radius \(\rho_{L}\). This is precisely the case for the database since for most energies \(\rho_{L}\sim 1\) (in the scaling described in Section (II.1)) while \(\lambda_{\perp}\in(0.2,20)\).

In Figs. (3a)-(3b) are shown typical CR trajectories in the perpendicular plane (Fig. (3a)) and in full 3D space (Fig. (3a)). Due to the fact that the pitch angle has been set to \(\mu=0.5\) in all simulations, the particles propagate along the parallel direction with larger velocities than the perpendicular guinding-center drifts. The Larmor rotation and the scattering due to fluctuations are obvious.

The running diffusion coefficient \(D_{\perp}(t)\) related to the mean-square-displacement of particles, \(\langle x^{2}(t)\rangle\), is computed accordingly with eq. (5). Another, more frequently used, formula for diffusion is \(d(t)=\partial_{t}\langle x^{2}(t)\rangle/2\). The reason for which \(D_{\perp}(t)\) is used here instead of \(d_{\perp}(t)\) is due to the fact that we capture Larmor rotations in our model. This component of motion exhibits large fluctuations at small times in \(d_{\perp}(t)\). \(D_{\perp}(t)\), on the other hand, due to its functional form, suppresses the effect of rotations very fast as well as the statistical fluctuations at long times. Nonetheless, both forms converge asymptotically to the same value. All these features can be seen in Fig. (4) where two typical time dependent diffusion profiles are shown for the same set of free parameters.

It must be emphasized that our database describes only pure diffusive regimes (in the perpendicular plane). In general, anomalous behaviour is possible as it was described, many times, in literature [7; 18; 16; 20; 9]. In fact, the finite parallel correlation length \(\lambda_{\parallel}=10\lambda_{\perp}\) is responsible for the decay of the Lagrangian correlation and, consequently, the saturation of the running diffusion to a constant value, thus, diffusive transport. The mechanism is simple: the parallel motion of CRs in fields with parallel dependence induces an effective decorrelation time \(\tau_{c}\approx\lambda_{\parallel}/v_{\parallel}\). Describing subdiffusive or superdiffusive transport remains an important task for future developments of ANNs.
### Convergence properties of the ANN

Once the database has been build and the ANN constructed (programmed), we are all set to start the training phase. A random configuration of initial weights and biases \(\{W_{i,j,k},\beta_{i}\}\) is chosen and the minimization procedure (via ADAM algorithm) begins. The random nature of initial \(\{W_{i,j,k},\beta_{i}\}\) is needed to avoid setting the ANN in a configuration point from which won't be able to converge to the global minimum. Due to the same reason, different initializations lead to distinct paths of convergence. There is also a question of weather the true global minimum is achieved [21] or the algorithm gets stuck a neighbouring local minimum. Regardless, the practical experience has shown that the asymptotic states are appropriate approximations of the real minimum.

In Fig. (5) it is shown a typical evolution of the error (8) relative to the norm of output values. The learning stage is fast, approximately \(5s\) on a typical personal CPU. One must note that the error function is almost always decreasing, but the evolution can manifest periods of quasi-plateau, as the one between \(400-900\) epochs in Fig. (5). For this reason, it is important to allow for long periods of minimization, in order to avoid confusing a local minima (the plateau) with the global one.

The number of particles propagated within a single simulation in this work is set to \(N_{p}=10^{4}\). Looking at Fig. (4), it might seem that this value is unnecesary since the diffusion coefficient \(D_{\perp}(t)\) saturates to a nice, constant, plateau without temporal fluctuations. While this is true, this does not mean that the ensemble is sufficiently well represented and, consequently, that the asymptotic value is statistically robust (without fluctuations). In other words, if \(N_{p}\) is too small, for two identical simulations we might obtain relatively different diffusions. This would have a detrimental impact on the

Figure 4: Typical perpendicular diffusion coefficients computed in two ways.

Figure 5: Error \(E\) evolution during the iterative training of the ANN with ADAM.

Figure 3: Typical CR trajectories in turbulent 2D magnetic fields shown in the perpendicular plane (a) and in full 3D geometry (b).
database, introducing a spurious numerical fluctuation and making the job of convergence much more difficult.

When the ADAM algorithm reaches the final plateau and an acceptable small error, one is ready to enter the testing phase. Now, the ANN, that is to say the values of \(\{W_{i,j,k},\beta_{i}\}\) obtained after minimization, can be used to evaluate diffusion coefficients for the remaining 10% input, \(X^{test}\), of the database that was not used in the training stage. The predicted values, \(Z^{test}_{m+1}\), are compared with exact output diffusions \(Y^{test}\) from the database. Fig. (6) shows the histogram of relative errors between ANN's predictions and real output, \(Err=100\left(Z^{test}_{m+1}/Y^{test}-1\right)\).

While there are points with relatively large errors (usually those of very small diffusion \(Y^{test}\ll 1\)), the overall accuracy of the ANN can be estimated with the second moment of the histogram and it is close to 5%.

When comparing computing times between test-particle simulations, \(t_{sim}\), and ANN, \(t_{ANN}\), the latter shows its true power: for a single diffusion coefficient on a two-processor CPU \(t_{sim}\approx 2s\) while \(t_{ANN}\approx 10^{-7}s\). Thus, ANNs are \(10^{7}\) times faster than numerical simulations.

### Making predictions with the ANN

The histogram of errors obtained in the testing phase (6) suggests a global error of \(\sim 5\%\). But this is not enough to ensure an acceptable accuracy of the ANN. In fact, there are points which have departed more than 10% from the exact values. Moreover, these errors might be a serious liability if they are clustered in some special way inside the parametric space of input variables.

In order to test whether this is the case, we must look for reduced spaces inside the database and evaluate if the ANN's predictions are a good fit for the real data. We do that on two distinct levels: first we look at two-dimensional dependencies of \(D_{\perp}\) varying two parameters at a time, than, we look at simple dependencies between \(D_{\perp}\) and only one parameter.

In Figs. (7a)-(7c) we compare ANN output (red dots) with exact values (blue surface, obtained with test-particle simulation). Fig. (7a) uses as variables the magnetic turbulence amplitude \(b\) and the perpendicular length \(\lambda_{\perp}\) while the energy is set to \(T=1\). Fig. (7b) sets \(b=0.2\) and varies \(\lambda_{\perp}\) and \(T\), while in Fig. (7c) \(\lambda_{\perp}=5\). As one can see, all predicted values (red) lie close to the surface of exact diffusions (blue) in all cases. Thus, there is no special parametric domain where errors tend to accumulate and the ANN exhibits similar accuracy across the entire database.

In Figs. (8a)-(8c) we set two parameters to constant values and vary only the third, as it follows: Fig. (8a) uses \(b=0.2,\lambda_{\perp}=5\), Fig. (8b) \(T=1,\lambda_{\perp}=5\) and Fig. (8c) \(T=1,b=0.2\). The overlap between predictions (red) and exact data (blue) are, again, consistently accurate across all simulations. Furthermore, at the level of such single parameter dependencies, we can understand how using an ANN tool could enhance our ability to in

Figure 6: Relative error distribution on the testing database.

Figure 7: ANN predicted diffusion (red dots) in comparison with exact, test-particle derived, data (blue surface) keeping one parameter constant at a time (\(T=1\)-(a), \(b=0.2\)-(b), \(\lambda_{\perp}=5\)-(c)).
vestigate and understand the physical processes at play. For example, having easy access to the red curve in Fig. (8a) allows one to observe the maxima in the diffusion profile and infer about the existence of two competing mechanisms in the influence of particle energy. One is related to the monotonic increase of parallel velocity \(v_{\parallel}\) with \(T\) which increases the transport, while the other is connected to the finite Larmor radius effects that tend to decrease diffusion [18].

## IV Conclusions

The present work described a methodology for building artificial neural networks designed for predictions of cosmic ray transport in turbulent magnetic fields. The ANN developed here has a standard architecture with two hidden layers and \(tanh\)/quadratic activation functions. It uses the ADAM algorithm for optimizations in the learning phase. The input data in the training/testing database consists of the values of free parameters for CRs and turbulence, chosen inside a (hyper)cube of convenience. The output is represented by associated diffusion coefficients. The values of the latter are evaluated using a transport model which is numerically tackled with the aid of test-particle simulations.

The learning stage showed fast and good convergence properties. In the testing phase, a good fit between exact and ANN predicted data was found with overall errors of \(\approx 5\%\). The predictive power of the ANN is proven by the dependencies of the transport coefficients on individual parameters. The most stringent feature of the network is its speed, since it is able to compute diffusion on bulk cases approximately \(10^{7}\) times faster than test-particle simulations.

One might argue that the overall database constructed here is too limited: the model of turbulent fields (2D with a specific spectrum and connected correlation lengths \(\lambda_{\parallel}=10\lambda_{\perp}\)) is too simple to represent most astrophysical regimes and the numerical range of parameters (energy, correlation lengths, field strengths) is quite narrow. This is all true, but the purpose of this work was not to exhaust wide physical regimes. It was the author's intention to present a proof of concept for the use of machine learning techniques in the field of astrophysics, applied to a particular problem.

Hopefully, the community will explore the use of this tool in future works using more sophisticated transport models and wider databases, thus, developing more potent ANNs. It is likely that, given the numerical effort required in building large databases, collaborations will be suitable.

_A word of caution._ Developing ANNs to evaluate quantities of interest is not real knowledge. It does not constitute insight into any physical processes or quantitative mechanisms. Nonetheless, for practical applications where predictions are important, ANNs might be the choice of preference. Finally, even for more academically oriented questions, having a fast tool at hand to explore different regimes could be invaluable in the quest for disentangling the physical mechanisms at work.

## Statement

The database and the ANN developed in this study are available upon reasonable request from the authors.

Figure 8: ANN predicted diffusion (red dots) in comparison with exact, test-particle derived, data (blue dots) keeping two parameters constant at a time (\(b=0.2,\lambda_{\perp}=5\)-(a), \(T=1,\lambda_{\perp}=5\)-(b), \(T=1,b=0.2\)-(c)).
## Acknowledgements

This research was partially supported by Romanian Ministry of Research, Innovation and Digitalization under Romanian National Core Program LAPLAS VII - contract no. 30N/2023.
