Cosmic-ray interactions with the Sun using the FLUKA code

The interactions of cosmic rays with the solar atmosphere produce secondary particle which can reach the Earth. In this work we present a comprehensive calculation of the yields of secondary particles as gamma-rays, electrons, positrons, neutrons and neutrinos performed with the FLUKA code. We also estimate the intensity at the Sun and the fluxes at the Earth of these secondary particles by folding their yields with the intensities of cosmic rays impinging on the solar surface. The results are sensitive on the assumptions on the magnetic field nearby the Sun and to the cosmic-ray transport in the magnetic field in the inner solar system.

The interactions of cosmic rays with the solar atmosphere produce secondary particle which can reach the Earth. In this work we present a comprehensive calculation of the yields of secondary particles as gamma-rays, electrons, positrons, neutrons and neutrinos performed with the FLUKA code. We also estimate the intensity at the Sun and the fluxes at the Earth of these secondary particles by folding their yields with the intensities of cosmic rays impinging on the solar surface. The results are sensitive on the assumptions on the magnetic field nearby the Sun and to the cosmic-ray transport in the magnetic field in the inner solar system.

I. INTRODUCTION
Cosmic rays entering in the Solar system after propagating for millions of years in the Galaxy can reach the planets and the Sun itself, producing emission of secondary particles, such as gamma rays and neutrinos, due to the interactions with the surfaces or the atmospheres of the celestial bodies.
The Moon [1,2] and the Earth [3] are both bright sources of gamma rays. Lunar and terrestrial gamma rays are originated from the hadronic interactions of cosmic-ray nuclei with the lunar surface and with the upper layers of the Earth's atmosphere, respectively. The Sun is also a bright source of high-energy gamma rays. While gamma rays from the Earth and from the Moon are originated from cosmic-ray nuclei, the solar gamma-ray emission consists of two components: the first one, called disk emission, is due to cosmic-ray nuclei interacting with the solar surface [4,5], and is localized around the solar disk; the second one, which is due to the inverse Compton scatterings of cosmic-ray electrons (and positrons) with the solar optical photons, is a diffuse component and extends up to tens of degrees from the Sun [4,[6][7][8].
Several attempts have already been done to calculate the secondary emission (e.g., gamma rays and neutrinos) due to the interactions of cosmic rays with the solar atmosphere (see for example [5,9]). In particular, the knowledge of such emission could be used to constrain exotic processes, such as the production of standard model particles in the annihilation of dark matter particles captured by the Sun [9][10][11][12][13]. * mazziotta@ba.infn.it Early predictions are based on semi-analytical calculations with the inclusion of solar magnetic field [5], while full numerical simulations for the production of neutrinos based on the Monte Carlo method have been performed in [14] and recently revisited and updated by Refs. [9]. However, in those Monte Carlo simulations, the effect of the magnetic field was neglected since the calculation was performed at high energies. The production of neutrinos is closely related to that of gamma rays in the solar disk, as both are originated from hadronic interactions of cosmic-ray nuclei.
In this work we have performed a full simulation with the FLUKA code to calculate the yields of secondary particles produced by the interactions of cosmic rays with the Sun. In particular, we have simulated the interactions of protons, helium nuclei and electrons impinging on the solar atmosphere in a wide energy range from 0.1 GeV/n to 100 TeV/n, while the energy of secondary particles has been simulated down to 100 keV. The low-energy region is extremely interesting for the proposed future gammaray telescopes [15][16][17], which aim to probe photon energy intervals extending well below the lower bound of that explored by the Fermi Large Area Telescope (a few tens of MeV) [18].
The present work is based on our previous ones, in which we evaluated, using FLUKA, the yields of secondary cosmic rays in the collisions of primary cosmic rays with the interstellar gas [19] and the lunar gamma-ray emission [2].
In simulating the interactions of cosmic rays with the solar atmosphere there are a number of important effects to be considered. First, the interplanetary magnetic field affects the spectra of cosmic rays reaching the Sun. Second, the strong heliospheric magnetic field nearby the Sun also affects the trajectories of charged particles in the solar atmosphere: in particular, the total path length increases with the intensity of the magnetic field, and this corresponds to an increase of the interaction probability and consequently to an increase of the cascades of secondary particles. Finally, the profile of the solar atmosphere needs to be accounted in detail, since the cascades usually develop from a low-density medium toward a denser medium; in addition, the yield of secondary particles far away from the Sun is also affected by the grammage along the line of the sight from the production point to the outer space.

II. SIMULATION SET-UP
The propagation and the interactions of cosmic rays with the solar atmosphere have been simulated with the FLUKA code [20][21][22]. FLUKA is a general purpose Monte Carlo code for the simulation of hadronic and electromagnetic interactions, used in many applications. It can simulate with high accuracy the interactions and propagation in matter of about 60 different species of particles, including photons and electrons from 1 keV to thousands of TeV, neutrinos, muons of any energy, hadrons of energies up to 20 TeV (up to 10 PeV when it is interfaced with the DPMJET code [23]) and all the corresponding antiparticles, neutrons down to thermal energies and heavy ions.
Hadronic interactions are treated in FLUKA following a theory-driven approach. Below a few GeV, the hadron-nucleon interaction model is based on resonance production and decay of particles, while for higher energies the Dual Parton Model (DPM) is used, implying a treatment in terms of quark chain formation and hadronization.
The extension from hadronnucleon to hadron-nucleus interactions is done in the framework of the PreEquilibrium Approach to NUclear Thermalization model (PEANUT) [24,25], including the Gribov-Glauber multi-collision mechanism followed by the pre-equilibrium stage and eventually equilibrium processes (evaporation, fission, Fermi break-up and gamma deexcitation).
FLUKA can handle even very complex geometries, using an improved version of the well known Combinatorial Geometry (CG) package, that has been designed to track correctly both neutral and charged particles, even in the presence of magnetic fields.
To simulate the interactions of cosmic rays with the solar environment we have implemented a spherical geometry, which includes the radial profiles of the chemical composition, of the density, of the temperature and of the pressure. In addition, we have included a model of the magnetic field in the region close to the Sun (r < R S = 2.5R ), which affects the interactions of CRs with the Sun, and a model of the interplanetary magnetic field (r > R S ), which affects the propagation of CRs to the Sun.
In our code we use a spherical reference frame centered on the Sun with polar axis (i.e. z-axis) corresponding to the Sun's rotation axis. In our simulation we adopt the value of R = 6.9551 × 10 10 cm for the Sun radius.
To evaluate the yields of secondary particles from the Sun we have simulated several samples of protons, electrons and 4 He nuclei with different kinetic energies impinging a sphere surrounding the Sun with an isotropic and uniform distribution. The particles are generated on a sphere of radius r = R S , which is the boundary between the inner and outer magnetic field regions. We have then studied the effects of the solar magnetic field on the yields of secondaries by varying the inner magnetic field. The primary kinetic energy values are taken on a grid of 97 equally spaced values in a logarithmic scale, from 100 MeV/n up to 100 TeV/n.

A. Solar composition
For the radial composition profile, describing the distribution of mass in the Sun as a function of the distance from the centre, we use the Standard Solar Models (SSMs) for the interior of the Sun provided by Ref. [26] (hereafter Model gs98) 1 . Since most of cosmicray interactions will take place close to the surface, we have extrapolated this model outside the Sun. Figure 1 shows the mass fractions of the main components as a function of the distance from the center of the Sun for the gs98 model. The main components are the hydrogen and 4 He, while the abundances of heavier isotopes are below 1%.  For the radial density, temperature and pressure profiles we use the model provided by Ref. [27] (hereafter Model S), since it extends up to about 500 km above the solar surface. We then extrapolate this model to higher altitudes, up to about 1400 km. We have also verified that the Model gs98 is very similar to the Model S up to r = R . Figure 2 shows the radial density (top panel), pressure (middle panel) and temperature (bottom panel) profiles. The Model S is shown with black points, the Model gs98 is shown with blue lines, and the extrapolation is shown with red lines.
In the FLUKA simulation set-up we have implemented 100 layers (i.e. shells) with different densities, divided in three sets equally spaced on a logarithmic density scale: the external 40 layers from about 10 −13 g/cm 3 up to 10 −3 g/cm 3 , the middle 40 layers from 10 −3 g/cm 3 to 10 −1 g/cm 3 and the inner 20 layers with density > 10 −1 g/cm 3 . In each shell we define a compound mixture material according to the mass composition, density, temperature and pressure profiles described above. In particular, we also set the temperature profile for the neutron cross section for the main isotopes (i.e. H, 3 He, 4 He and 12 C). The temperature has an effect in the capture of neutrons that produce the gamma-ray line of 2.2 MeV.

B. Inner magnetic field
The magnetic field near the Sun is complex and strongly time-dependent, and the coronal magnetic field is usually extrapolated from the observed photospheric fields. A widely adopted model is the potential field source surface (PFSS) model [28,29], in which the field is purely radial at the source surface with radius R S (the source surface radius is approximately 2.5 solar radii).
We use the values of the three components of the coronal magnetic field B r , B θ , B φ at 51 heights between the photosphere (i.e. r = R ) and the source surface with radius R S = 2.5R , calculated using the photospheric magnetic field observations [30][31][32]  Magnetic field intensity as a function of the Carrington longitude and latitude angles at r = R for the CR 2111.
In the present work we assume that the magnetic field inside the Sun is always equal to that at r = R . The intensity of the magnetic field on the solar surface is shown in Fig. 3 for the Carrington rotation (CR) 2111, covering the period from 2011-06-05 17h to 2011-07-03 00h. We point out that in small regions of the solar surface the field intensity can even exceed 10 G.

C. Interplanetary magnetic field
The propagation of cosmic rays in the solar system is usually performed by means of the Parker model [39] for r > R S . The components of the interplanetary magnetic field (IMF) are given by: The angle ξ is defined as: where θ is the polar angle, ω S = 2.69 × 10 −6 rad/s is the angular velocity of the Sun (corresponding to a period of about 27 days) and v SW is the velocity of the solar wind (its typical value is 400 km/s). At the distance R S the B φ (and the B θ ) components are null to ensure the continuity with the PFSS model (see Sect. II B).
In the previous equations the intensity of the field B E is given by B E = B 0 / 1 + tan 2 ξ(R E , π/2), where B 0 is the intensity of the magnetic field at the Earth (its typical value is about 5 nT), and R E = 1 AU is the Sun-Earth distance. The constant f is given by: where H is the Heaviside function and the angle θ is the polar position of the heliospheric current sheet (HCS) defined as: where φ is the azimuth angle and we have indicated with α the tilt angle, i.e. the maximum latitude of the HCS; finally, the ± sign in Eq. 1 depends on the polarity of the magnetic field. Figure 4 shows the time evolution at the Earth of the magnetic field B 0 , of the solar wind velocity v SW and of the tilt angle α averaged in each CR number from 2008 to 2018. The values of the tilt angle α and its polarity are extracted from the Wilcox Solar Observatory public website [40], while the intensity of the magnetic field at the Earth B 0 and the velocity of the solar wind v SW are derived from the observations of the ACE satellite extracted from the NASA/GSFC's OMNI dataset [41,42].
In our simulation we have implemented the magnetic field configurations corresponding to a few Carrington rotations between 2011 and 2014, when the maximum of the solar cycle 24 occurred.

III. SIMULATION RESULTS
The yield of secondary particles from the i-th species of cosmic ray primaries (here i = p, e − and 4 He) Y s,i (E s |E k ) is calculated by using the multiplicity of the secondary particles escaping from the generation surface. The yield is calculated as: where N i (E k ) is the number of primaries of the ith species generated with kinetic energy E k (in units of GeV/n) and N s,i (E s | E k ) is the number of secondary species s with energies between E s and E s + ∆E s produced by the primaries of the type i. Fig. 5 shows the yields of gamma rays produced by protons (top panel), helium nuclei (middle panel) and electrons (bottom panel). The yields have been calculated as a function of the secondary energy bins with logarithmic step, by using fine bins below 10 MeV (i.e., 32 energy bins per decade), and above with coarse energy bins (i.e., 8 bins per decade). The differential intensity of secondary particles (in units of particles GeV −1 cm −2 sr −1 s −1 ) emitted from the Sun is given by: where I i (E k ) is the intensity of the i-th species of cosmicray primaries at the Sun.
The flux of secondaries observed by a detector at Earth (in units of particles GeV −1 cm −2 s −1 ) is given by: where F(E s ) is the fraction of secondaries with energy E s which are able to reach the Earth's orbit from the Sun.
In our simulation we assume that the Earth's orbit lays on a sphere centered on the Sun with radius r = R E . We point out here that not all secondaries emitted outwards from the Sun are able to reach the Earth. Charged particles are deflected by the interplanetary magnetic field and, depending on their energy and initial direction, can be sent back to the Sun without reaching the Earth's orbit. In addition, there are some species of unstable secondaries, such as neutrons, which can decay during their journey from the Sun to the Earth. In these cases, the fraction of secondaries reaching the Earth will be F(E s ) ≤ 1. On the other hand, for the secondary gamma-ray and neutrino we assume F(E s ) = 1. 2 Cosmic rays impinging on the solar atmosphere are those which can reach the Sun from the interplanetary space. Hence the intensities I i (E k ) of the various cosmicray primaries in eq. 6 are those at the surface of the generation sphere of radius R S , which differ from those measured at Earth, since not all cosmic rays reaching the Earth are able to continue their journey to the Sun.
To evaluate the cosmic-ray intensities at the Sun we have used the custom code HelioProp [43,44] 3 , which describes the transport of cosmic rays in the Solar System. We have simulated sets of pseudo-particles injected on the surface of the generation sphere with an isotropic and uniform distribution. The pseudo-particles are followed backwards in time during their propagation until they reach a sphere of radius R E [45][46][47]. Their survival probabilities are used to scale the measured intensities of cosmic rays at Earth in order to properly set the intensities I i (E k ) in the right-hand side of eq. 6.
In our simulations we assume that the intensity of cosmic rays measured at the Earth is the same across a sphere or radius R E = 1 AU. At low energies (< 10 GeV) this assumption could be not valid because of a possible dependence on the charge sign of the propagation of cosmic rays from the outer space to 1 AU [43]. However, this effect is not expected to produce significant changes in our results, since only a small fraction of low-energy CRs are able to reach the Sun.
We use the cosmic-ray intensities at Earth measured by AMS-02: the proton intensity is taken from Ref. [48], the helium intensity is taken from Refs. [49,50] and the electron 4 intensity is taken from Ref. [51]. We also use the AMS-02 spectra measured for different Bartels' rotations [52,53]. 5 Since the AMS-02 spectra are available starting from about 0.4 GeV/n, we have extrapolated the data down to 0.1 GeV/n by fitting the measured intensities with a function given by [54]: For the proton and helium we fit the data points up to the break energy around 200 GeV/n; then for larger energies we include a smooth break with a harder spectral index, as indicated by Refs. [48,49]. In the case of electrons we also take into account the DAMPE data [55] including a break at about 900 GeV.  Figure 6 shows the results of the fitting procedure with the experimental data points corresponding to the CR 2111, covering the period from 2011, June 5 th to 2011, July 3 rd . In Fig. 6 we also show the modulated spectra at the Sun, evaluated from those at the Earth with Helioprop. Figure 7 shows the gamma-ray fluxes at the Earth evaluated with our simulation set-up for four different Carrington rotations (2111, 2125, 2138 and 2152) spanning the period from June 2011 to June 2014, that covers the AMS-02 measurements. The calculated gamma-ray fluxes are slightly different at low energies (< 1 GeV), due to the effect of the heliospheric magnetic field that affects both the cosmic-ray intensity at the Sun and the secondary yields. Finally, fig. 8 shows the gamma-ray flux at Earth obtained by averaging the fluxes calculated in the four different Carrington rotations.
In the Figs. 7 and 8 we show the total gamma-ray flux at Earth and the contributions of photons produced by the interaction of protons, Helium and electron primaries separately. The tipical contributions of protons, helium and electron primaries to the total gamma-ray fluxes are of about 74%, 24% and 2%, respectively.
The gamma-ray flux at the Earth exhibits two sharp peaks at 511 keV and at about 2.2 MeV, due to the positron annihilation and to the neutron capture (in the hadronic interactions) respectively. These two lines could be used as reference to calibrate the low energy gamma-ray telescope proposed for the next decade, such as ASTROGAM [15,16] and AMEGO [17]. At energies above tens of GeV the fluxes exhibit some fluctuations that are due to the limited statistic in the simulated data sets. 6 In figs. 7 and 8 we also show the experimental results for the disk component with 1.5 years [4] and 9 years [56] datasets collected by the Fermi Large Area Telescope (LAT) [18]. We stress here that these measurements have been performed in different time windows from the one covered by our simulation.
We have also cross-checked these results by backpropagating from the Sun to the Earth each particle simulated with FLUKA. Given a cosmic-ray primary at the Sun, a particle with opposite charge and with opposite direction is back-propagated from the generation sphere of radius R s to the sphere of radius R E . If this particle is able to reach the Earth, the primary particle assigned a survival probability P surv = 1, otherwise it is assigned P surv = 0. With this procedure, the secondary yield can be calculated as: where N i (E k ) is the number of primaries of the i-th species generated with kinetic energy E k and N s,i (E s | E k , P surv = 1) is the number of secondaries of the i-th species with energies between E s and E s + ∆E s produced by the primaries of the type i with P surv = 1. In this way, the secondary spectra at the Earth can be calculated inserting in the right-hand side of eq. 6 the proton, Helium and electron intensities measured at the Earth. Using this procedure we find the same results as when we evaluate the intensities of primary CRs at the Sun with HelioProp.
In Tab.I we show the integral fluxes of gamma rays at Earth above 100 MeV, 1 GeV and 10 GeV respectively, for the four CRs considered in this section.
The integral fluxes decrease with the CR number, as the Sun approaches to its maximum activity.
The secondary productions at high energies occur close the solar surface, where the secondary are emitted in a low-density medium in the forward direction with respect  [5]; gray points: 1.5 years Fermi-LAT data [4]; dark red points: 9-years Fermi-LAT data [56].  Fig. 7. Color lines and data points have the same meanings as those in Fig. 7.
to the high-energy primary particles. However, the combination of the solar magnetic field with the solar atmosphere density profile can affect the emission, even at high energies. Figure 9 shows the gamma-ray flux seen at the Earth as a function of the angle of sight for two different energy bins, i.e. [0.1, 0.133] GeV and [1, 1.33] GeV. We also show the corresponding spatial emission maps centered on the Sun and built with the HEALPix pixelization [57] 7 . The emission is mainly located nearby the solar surface and it becomes much narrow at higher energies. Figure 10 show the intensity of different species of secondaries (muon neutrinos and antineutrinos, electron neutrinos and antineutrinos, neutrons, electrons and positrons) produced by the interactions of cosmic rays with the Sun, evaluated on the generation sphere. The values of the intensities are obtained by averaging the results in the four Carrington rotations mentioned above. To calculate the fluxes at the Earth, the decays of unstable particles during their journey from the Sun to the Earth should be taken into account, as well as the propagation of charged particles in the IMF. We also remark here that in the calculation of the neutrino and antineutrino fluxes we did not include their interactions 2 in the Sun and their possible oscillations [9]. The neutrino intensity is similar to the calculation of Ref. [9], and its value is similar or larger than the intensity of neutrinos produced in cosmic-ray showers in the Earth's atmosphere (see for example [58]).

IV. EFFECT OF THE MAGNETIC FIELD ON THE SECONDARY YIELDS
The secondary emissivity of the Sun is strongly dependent on the intensity of the magnetic field close to the solar surface. To study this effect we have  3. enhanced B field configuration near the Sun (r/R < 1.01) following the BIFROST model [59][60][61], i.e. we increase the the original PFSS maps near the Sun to follow the BIFROST profile 8 . Figure 11 shows the gamma-ray fluxes at the Earth with the four different configurations of the PFSS magnetic field, i.e. the nominal model and the three alternative models illustrated above.
The gamma-ray flux without magnetic field is significantly enhanced at low energies with respect to the flux in presence of magnetic field, while for gammaray energies above 10 GeV the flux increases as the magnetic field increases. If the solar magnetic field is suppressed, low-energy cosmic rays can reach the Sun surface inducing a shower which produce secondary particles in the outer space. The presence of a solar magnetic field reduces the probability that low-energy cosmic rays can reach the Sun, but increases the probability of interaction for high-energy cosmic rays, 8 The BIFROST simulation is available for a limited region of the Sun. The enhanced factor is about 25 at the solar surface.
since they move along curved trajectories in a strong and non-uniform magnetic field and their path length increases as the magnetic field increases. This effect is well visible when comparing the gamma-ray fluxes with B = 0 (top left panel in Fig. 11) with the one with enhanced B field configuration (bottom right panel in Fig. 11)

V. DISCUSSION AND CONCLUSIONS
We have implemented a full simulation with the FLUKA code to calculate the yields of secondary particles produced by the interactions of primary cosmic rays with the solar atmosphere. Our simulation includes the current state-of-art models and data available to describe the solar atmosphere, the magnetic field nearby the Sun and in the interplanetary space.
The FLUKA toolkit provides a detailed simulation of hadronic and electromagnetic interactions in the matter in a wide energy range, with complex geometries and even in presence of magnetic fields. The geometry used in the present work is quite flexible, and it can be used for any other configuration.
The solar atmosphere and its chemical composition have been taken from the SSM gs98 and from the model S, with some extrapolation in the chromosphere region. However, the average density should drop below 10 −13 g/cm 3 at an altitude of about 1400 km from the solar surface, where the interaction probability should be negligible. We have also used the model ags09 9 and we found very similar results. νe +νe. The prediction of neutrinos by Ref. [9] is also shown. In addition, we include results from gamma-ray observations from Fermi-LAT [4,56], gamma-ray HAWC's 95% limit [62] and neutrinos IceCube's 90% limit [63].
The magnetic field adopted nearby the Sun is the one predicted by the PFSS model by using the synoptic map from HMI 720sline-of-sight magnetograms collected over 27-day solar rotations in a high-resolution Carrington coordinate grid. We have studied the effect of the magnetic field on the secondary yields by changing the original values of the PFSS maps. Indeed, we found that the yields are strong affected by the intensity of the magnetic field, even in the high energy region of the emission in the outer space. Recent developments of numerical solutions of a magneto-hydrodynamical (MHD) model together with the current observation of the Parker Solar Probe [64,65] could provide new insights to get a realistic description of the plasma dynamics and of the magnetic field nearby the Sun.
The calculated solar gamma-ray flux at the Earth has been compared with the Fermi-LAT data on the disk emission in different time windows with respect to those used in the present simulations. The detected gammaray emission from the solar disk above 1 GeV shows a harder spectrum (∼ E −2.2 γ ) than the cosmic-ray spectrum 9 http://www.ice.csic.es/personal/aldos/Solar_Data_files/ struct_b16_agss09.dat (∼ E −2.7 p,He ). This behaviour would require a high-intensity magnetic field configuration nearby the Sun, up to a factor 20 larger than the field predicted by the PFSS model used in the present simulation.
In Fig. 12 we show the predictions of our simulation about the fluxes at Earth of gamma rays and neutrinos in the CR 2111. The simulations have been performed with the standard PFSS solar magnetic field and with the enhanced one according to the BIFROST profile. The gamma-ray production is higher in the case of the more intense magnetic field, while the effect of the magnetic field on the neutrino production seems negligible. This could be a signature that the interactions in high magnetic field occur in the higher layers of the lowdensity solar atmosphere, resulting in an enhanced high energy gamma-ray emission. Anyway, above 100 GeV the predicted gamma-ray flux is below the HAWC's limit [62], even in case of the enhanced solar magnetic field. The predicted neutrino flux is lower than the calculation by Ref. [9], in particular below 1 TeV, due to the effect of the nearby solar magnetic field and is well below IceCube's 90% limit [63].
The magnetic field should also affect the inverse Compton gamma-ray emission close to the Sun, since the electrons should move along curved trajectories, whose lengths determine the interaction probability with the intense optical photon field. In this way the inverse Compton emission should also be peaked close to the solar surface, and might be not well separated by the disk emission due to the interaction of cosmic rays with the solar atmosphere. To get a complete picture of the solar gamma-ray emission the inverse Compton scattering needs to be calculated in presence of strong and irregular magnetic field. This simulation is beyond the scope of the current paper and it deserves a dedicated work.