{'fraction_non_alphanumeric': 0.021812749003984062, 'fraction_numerical': 0.00702191235059761, 'mean_word_length': 5.021289355322339, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
CRT: A numerical tool for propagating ultra-high energy cosmic rays through Galactic magnetic field models

Deflection of ultra high energy cosmic rays (UHECRs) by the Galactic magnetic field (GMF) may be sufficiently strong to hinder identification of the UHECR source distribution. A common method for determining the effect of GMF models on source identification efforts is backtracking cosmic rays. We present the public numerical tool CRT for propagating charged particles through Galactic magnetic field models by numerically integrating the relativistic equation of motion. It is capable of both forward- and back-tracking particles with varying compositions through pre-defined and custom user-created magnetic fields. These particles are injected from various types of sources specified and distributed according to the user. Here, we present a description of some source and magnetic field model implementations, as well as validation of the integration routines.


Introduction
Ultra-high energy cosmic rays are thought to be charged protons or nuclei, due to constraints on the photon flux at the highest energies [1]. These particles will therefore be deflected by galactic and extragalactic magnetic fields during propagation from their sources. Any type of UHECR source identification must account for this deflection. Studies including Stanev [2], Medina Tanco et al. [3], Harari et al., [4], Yoshiguchi et al. [5], Prouza andŠmída [6], Takami et al. [7], and Kachelrieß et al. [8] have used the method of backtracking cosmic rays through various GMF models as a source discrimination tool. However, these studies investigated particular configurations of the GMF with fixed quantities for parameters such as the overall normalization and scale heights. Additionally, the propagation codes used in each study were individually developed making a direct comparison of the propagation accuracy between the codes difficult.
Two commonly cited public propagation codes for cosmic rays are GALPROP 1 [9] and CRPropa 2 [10]. However, the focus and design of these codes differ from CRT, which focuses on minimizing computation time for tracking UHECRs of arbitrary composition. These codes are designed to investigate the effects of interaction and energy loss during cosmic ray propagation. While GALPROP contains a detailed exponential model of the GMF, particles are propagated by solving the diffusion equation which is more suited for low energy (GeV-TeV) cosmic ray studies in the GMF. CRPropa numerically integrates the equations of motion but implements only extragalactic turbulent magnetic fields. In particular, the CRPropa code is not optimized for convoluted GMF models and only propagates cosmic rays with Z = 1.
Presented here is a public tool for propagating UHECR events through various models of the GMF. This program implements components of field models used in previous studies and many possible source distributions. The user is capable is specifying many simulation parameters such as tracking method (forward or back), magnetic field model and source configurations, and integration routine settings. This is accomplished through a simply-designed interface. Output includes useful information about the injection and observation site coordinates and deflection magnitude. Header information in the output file provides a record of the simulation settings. The code is actively maintained with revision control using subversion 3 . A website 4 containing public releases of the code is accessible under credentials 5 . Complete simulation information beyond that described here can be found in README text files accompanying each version.

Overview
The numerical tool CRT is written in C++ and has only a single dependency on the GSL random number generation libraries [11]. CRT is capable of forward-and back-tracking cosmic rays from their injection sources to a detector or magnetic field boundary. Command line options set the global parameters for a specified program run, such as backtracking status and detector size. The desired injection sources and magnetic field models, as well as their parameter values, are specified by the user in an external configuration file processed upon program execution.
UHECR events are injected at user-specified sources and propagated through a user-specified magnetic field by integrating the relativistic Lorentz equation. If the backtracking mode is selected, the source is located at the Galactic position of the Sun within the GMF model. CRT implements 5 th order adaptive Runge-Kutta integration routines for determining the cosmic ray trajectories. The term "adaptive" arises from checks performed prior to completing an integration step. The step is recalculated with a smaller stepsize if the induced errors are larger than the user-specified tolerance. These routines are more fully described in Numerical Recipes [12]. Trajectories can be recorded step-by-step or at specific locations 6 , such as the boundary of the detector. Cosmic rays are propagated until they reach the boundary of the field model, integration time limit, or detector.
Energy losses are not considered during propagation. The GZK process [13,14] involves cosmic ray interactions with the cosmic microwave background (CMB) resulting in photo-pion production for protons or photodisintegration of heavy nuclei and is the dominant energy loss on UHECRs due to the CMB constituting the dominant particle species in extragalactic space. The same process can occur within the Galaxy due to interactions with the ambient starlight and infrared radiation from stars and warm dust. However in the Galactic center, the infrared photon density for energies above threshold energy is estimated to be of order 1 cm −3 or less (see, e.g., [15]) and decreases with galactocentric distance. The CMB remains the dominant photon field n cmb = 411 cm −3 even in galactic interstellar space, but UHECR interactions are infrequent due to the reduced pathlength. UHECRs may travel of order tens of kpc through the Galaxy as opposed to Mpc in intergalactic space. Above threshold energy, the interaction cross-section σ pγ is roughly constant at 10 −1 mb. The probability for a single interaction is P = Ln cmb σ pγ , so for a typical pathlength through the Galactic center L ≈ 10 kpc, P ≈ 10 −3 . Figure 1: Detection (blue) and injection (orange if forwardtracking; blue if backtracking) disc implementation in CRT. This image depicts a scenario where a cosmic ray is forwardtracked and detected from an extragalactic source. The black line traces the particle trajectory from the injection disc to the detection disc. Effects from energy loss by synchrotron emission due to the strong magnetic fields observed in the Galactic center region (B ≈ tens of mG) (see, e.g., [16,17]) are also negligible at EeV energies and higher. Energy loss due to inverse Compton scattering is unimportant at ultra-high energies.
CRT implements a single detector scheme. It is constructed as a disc located at the Sun's position in the Galactic plane and oriented such that the surface normal always lies parallel to the particle trajectory. The detector is shown as a blue disc at the center of the coordinate system in Figure 1.
CRT also includes additional executables specifically designed for testing the integration routines and parameterizing the magnetic field models. Use of the integration accuracy test program is demonstrated in Section 5. For a user-specified GMF configuration, the executable getfield outputs information about the magnetic field useful for producing three dimensional grid plots of the field strength and direction. An example plot is shown in Figure 2.

Magnetic Field Models
The public releases of CRT incorporate descriptions of many magnetic field models, but is modularly designed so that the user may easily write and include additional field models, a central feature in its design. It is possible to use any number and combination of default and user-written models simultaneously during the execution of the program; the user specifies the desired configuration in the external configuration file.
Field models included in the public releases are: Logarithmic Spiral, Dipole, Toroidal, Ring, and a simple random field. The user is capable of specifying every parameter value of these field models. These individual field models can be combined to create more complex field models, such as that studied by Prouza andŠmída [6] for example, and also combined with any usercreated field models. These default models are fully described in documentation accompanying each release version. The use of these magnetic field models is entirely optional by the user. Only models specified by the user in the external configuration will be used during program execution.
Poloidal and Toroidal fields are thought to exist as components of the Galactic halo field away from the Galactic plane and so are included. Logarithmic Spiral field models have seen extensive use in previous cosmic ray and radio starlight studies (see, e.g., [2][3][4]6,18,19]). The Logarithmic Spiral field model incorporates different symmetry parameterizations about the Galactic rotation axis and plane, respectively, and radial dependencies from two commonly studied field model : Stanev [2] and Harari-Mollerach-Roulet (HMR) [4]. In total, CRT includes the descriptions of 8 unique types of Logarithmic Spiral models available for use. The Galactic magnetic field also appears to have a significant random component (see, e.g., [16,18]), which can be approximated with a simple random field model. The simple random field in CRT is constructed by populating individual spherical magnetic cells each with a single direction and magnitude chosen from gaussian distributions. The gaussian distributions are described by parameters specified by the user. Cell boundaries are respected during individual cosmic ray propagation, ensuring a proper magnetic field calculation at each position.
It is important to reiterate that the user has complete freedom in choosing which magnetic field models to include during program execution. The user could select every field model listed above, or two, or even none and instead choosing to recompile CRT with their own source code. The creation of any user-written field model simply requires constructing and specifying the field model in a fashion similar to the included models. The construction can be as simple as amending the radial dependence of the field strength equation of a particular default spiral model or a user implementing a highly complex field model in new source code.

Sources and Spectra
CRT includes many possible source descriptions allowing for a variety of source distributions. As with the magnetic field models, CRT is capable of incorporating any combination of included and user-written sources and energy spectra during the program execution. The user is capable of specifying all parameter values of every source type. All sources inject cosmic rays with energies sampled from a power law distribution, Some source types are capable of sampling energies from a broken power law distribution described by two spectral indices and a break energy. Cosmic rays injected from any source type are capable of both forward-and back-tracking. One of the included source types, Isotropic sources generate events with initial velocity directions sampled isotropically across the sky. Once a direction has been selected, the injection site is determined within a disc of user-specified size lying within the tangent plane to that direction. In the case of forward-tracked cosmic rays, events are injected just beyond the magnetic field model boundary from a random location on the disc in order to mimic an extragalactic flux.
The user is capable of specifying the sizes of the injection disc for each source in the external configuration file and the detection disc as a runtime command line parameter. In a forward tracking scenario, injection disc radii larger than roughly kpc effectively mimic an impinging extragalactic cosmic ray flux on the Galaxy. Detector discs with radius of order 0.1 kpc efficiently detect forwardtracked events. The Sun is estimated to be roughly 20 pc above the Galactic plane, as inferred from stellar surveys (see, e.g., [20][21][22]). This detector radius value accounts for the error on the Sun's galactic position and gives a good compromise between accumulating detected events and properly simulating the local GMF. Figure 1 illustrates the orientation of the discs during runtime. For a simulation incorporating forwardtracking, the orange disc represents the injection disc of an extragalactic source. A particle trajectory is shown originating from a randomly selected position on the injection disc, propagating through a magnetic field, and impacting the detector disc. At each integration step, the detector disc is oriented such that its normal is parallel to the particle trajectory.
The Isotropic source type is one of many included source types in the public release. Additional source types are capable of injecting particles from any location or adjusting the injection direction to account for estimated deflection by a simple extragalactic turbulent magnetic field model 7 . There are also source types that incorporate sky coverage of earth-based detectors 8 . This coverage depends on the declination and limiting zenith angle of the detector, which can be easily changed within the source code. These sources are identical to other source types but their source code incorporates the additional coverage factor in selecting event directions.
All source types are more fully described on the website and in documentation accompanying the source code.

Validation
The accuracy of the integration routines has been tested by tracking charged particles in a uniform field. Particles traveling initially perpendicular to the field will execute simple circular motion whereas particles traveling initially parallel will propagate undeflected. The integration accuracy can be estimated by determining the difference between the exact and calculated positions and velocities as functions of time. The residuals are dR = | R exact − R calc | and dV = | V exact − V calc |. Figures 3 and 4 show the residuals as functions of time for 1 EeV protons injected into a uniform 1 µG magnetic field. The initial velocity vectors are parallel and perpendicular to the magnetic field, respectively.
In the parallel injection case (Figure 3), the particle propagates undeflected and the velocity residual always zero. The position residual at the beginning of propagation is dominated by 7 CRT does not propagate the particles through extragalactic space. For small scale fields, the deflection manifests as a gaussian smearing of the cosmic ray directions. 8 Available only for backtracking.  numerical accuracy. In this simulation, the particle travels parallel to a coordinate axis, so that at large times, the exact and calculated coordinate values are each large values. The calculation of the residual is the difference of two large values, leading to rounding uncertainty at large times.
The perpendicular case ( Figure 4) depicts a different scenario. The velocity residual shows the accumulated effect of numerical accuracy over time. The residual decreases at particular times, but generally shows a trend toward increasing values. The position residual is directly influenced by the velocity residual since the integration routine uses the calculated velocity values to determine the position. The behavior of the position residual plot closely mimics that of the velocity residual with negligible time delay. The residual value is initially near numerical accuracy, then slowly begins increasing around the same time as the velocity residual.
The position residuals in both injection cases are well below 10 −10 for an integration time of 10 6 yr, during which the total pathlength of the particles is roughly 300 kpc. The diameter of the Milky Way as resolved by its distribution of stars and gas is about 20-30 kpc and the magnetic field is thought to be roughly the same size. The Larmor radius for a UHECR in a uniform field is, where Z and E are the cosmic ray charge and energy in EeV, and B the field strength in µG. At E ≈ 10 EeV the Larmor radius becomes comparable to the Galactic diameter and only high charge values or magnetic fields with considerable field strengths can contain the UHECRs.
In realistic GMF models, the field vector varies significantly in strength and direction over length scales of 10 2−3 pc such that deflection induced from one region of magnetic field may be destructively interfered by the next region. These field models will influence the cosmic ray trajectory in a manner somewhat similar to random fields, so that the trajectory still maintains quasi-rectilinear behavior even at relatively low energy (E ≈ few EeV).

Example Simulations
Here we illustrate some useful applications of the output of CRT. Upon completion of tracking a particle, information is recorded including the initial and final conditions of the particle, the deflection magnitude, source properties, and detection status. The header information contains complete information about the simulation settings as well as the program call at the command line.
In the first example, 10 4 protons are backtracked from an isotropic observed distribution with a differential energy spectrum dN/dE ≈ E −3 with minimum and maximum energies E min = 50 EeV and E max = 200 EeV. The Galactic magnetic field is represented by the HMR BSS A configuration shown in Figure 2. The field experiences exponential attenuation away from the Galactic plane with two scale heights: 0.5 kpc and 3.5 kpc. Figure 5 shows a Hammer-Aitoff projected skymap in Galactic coordinates of the density of injection directions. These directions are selected from an isotropic sky distribution. Figure  6 indicates the inferred source direction of the particles based on the direction of the velocity vector upon crossing the magnetic field model boundary. Structure is apparent at specific Galactic longitudes, particularly ℓ ≈ 45 • , ±10 • , −50 • , due to the semi-regular spiral structure towards the central region of the field model. The nature of the field model induces the low energy particles to spiral along the magnetic field lines before exiting. The magnetic field model also deflects nearly all particles out of the Galactic plane due to increased pathlength across the plane through regions of large field strength. Figure 7 shows the deflection magnitude as a function of energy for all particles. The red line serves to guide the eye along a curve following E −1 behavior, such as that expected by energy-dependent deflection in a uniform field. It is apparent at low energy that some events experience less deflection than would be expected, as these particles are lensed along magnetic field lines into the nearest structure seen in Figure 6.
In a second example, 10 4 protons are forward-tracked from an extragalactic isotropic source distribution through the same magnetic field model as above. The injection site for each particle is randomly selected within an injection disc of radius 1 kpc. The detector radius is set to 0.1 kpc. Figure 8 shows the sky distribution of observed events. Figure 9 shows normalized count histograms of the energies for detected and non-detected events.

Conclusion
We have presented a publicly available numerical tool for tracking cosmic ray events from various source distributions through Galactic magnetic field models. CRT incorporates components of basic models of the Galactic magnetic field and many possible source distributions from which to create and propagate cosmic ray events. The modularity of CRT allows for users to easily implement modules for new magnetic field models and sources. The integration accuracy of the program has been tested and deemed acceptable for most simulation configurations.