# Solar disk gamma-rays emission via synthetic magnetic field from photosphere to low corona

Eleonora Puzzoni, Federico Fraschetti, Jozsef Kota, Joe Giacalone

► To cite this version:

Eleonora Puzzoni, Federico Fraschetti, Jozsef Kota, Joe Giacalone. Solar disk gamma-rays emission via synthetic magnetic field from photosphere to low corona. The Astrophysical Journal Letters, 2025, 989 (1), pp.L3. 10.3847/2041-8213/adf2a6. hal-05193572

HAL Id: hal-05193572

https://hal.science/hal-05193572v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Solar Disk Gamma-Ray Emission via Synthetic Magnetic Field from Photosphere to Low Corona

Eleonora Puzzoni<fnref n="1" />,<fnref n="2" />✉, Federico Fraschetti<fnref n="1" />,<fnref n="3" />✉, József Kóta<fnref n="1" />✉, and Joe Giacalone<fnref n="1" />✉

Received 2025 February 25; revised 2025 July 10; accepted 2025 July 19; published 2025 August 5

## Abstract

Gamma-ray emission in the GeV–TeV range from the solar disk is likely to arise from collisions of galactic cosmic rays (GCRs) with solar atmospheric plasma. In a previous study, we demonstrated that closed turbulent magnetic arcades trap GCRs efficiently, leading to a gamma-ray flux consistent with the Fermi-HAWC observations (from ~0.1 GeV to ~1 TeV). Here, we model a synthetic magnetic field with a static, laminar structure of open field lines in the chromosphere, increasingly braided near the solar surface, with a scale height of ~$10^{-2}$ R⊙. The height-dependent increase in magnetic field line braiding is modulated by an exponential scalar function, mimicking the bending of the photo- and chromospheric magnetic field revealed by polarimetric observations and reproduced by MHD simulations. Employing 3D test-particle numerical simulations, we investigate how distorted magnetic field lines affect the gamma-ray production by injecting GeV–TeV protons into both magnetically laminar and braided regions. We find that with the chosen spatial resolution, this synthetic magnetic field can account for the >10 GeV gamma-ray spectrum observed by Fermi/HAWC. A rebrightening between approximately 30 and 100 GeV (following a ~30 GeV spectral dip) suggests an enhanced confinement within the photo-/chromospheric layer by stronger braiding.

*Unified Astronomy Thesaurus concepts: Solar gamma-ray emission (1497); Galactic cosmic rays (567); Solar magnetic fields (1503); Solar atmosphere (1477)*

## 1. Introduction

Interest in high-energy gamma-ray emission from the solar disk has grown in recent years following the foundational work of D. Seckel et al. (1991), hereafter SSG91, which introduced the first theoretical model based on gamma-ray production from galactic cosmic rays (GCRs) mirroring within magnetic flux tubes and emitting gamma-rays in the interaction with ambient plasma of the solar atmosphere. However, observations from Fermi-Large Area Telescope (LAT, an imaging gamma-ray detector spanning energies from ~20 MeV to ~300 GeV; see W. B. Atwood et al. 2009) revealed several limitations in the model: the predicted ~1 GeV flux was a factor of 7 lower than the observations (see A. A. Abdo et al. 2011). Second, while the SSG91 model predicted an abrupt cutoff at ~5 GeV, Fermi-LAT observed gamma rays from the solar disk up to 100 GeV (T. Linden et al. 2018), and the High-Altitude Water Cherenkov Gamma-ray Observatory (HAWC), which detects gamma rays from 100 GeV to 100 TeV (A. U. Abeysekara et al. 2013), recorded gamma rays up to ~2.6 TeV (A. Albert et al. 2023). Other notable features in the solar disk gamma-ray emission include its anticorrelation with the solar cycle, initially detected by K. C. Y. Ng et al. (2016) and subsequently confirmed by T. Linden et al. (2018, 2022) and A. Acharyya et al. (2025). It was observed that the gamma-ray flux can be up to a factor of 2 lower during solar maximum compared to solar minimum. T. Linden et al. (2022), utilizing Fermi-LAT data spanning the entire 11 yr solar cycle from 2008 to 2020, and B. Arsioli &

E. Orlando (2024), employing Fermi-LAT data from 2008 to 2022, both identified energy-dependent anisotropies in the gamma-ray emission from the solar disk. In particular, T. Linden et al. (2022) found that emission at 10–50 GeV is nearly uniform across the solar disk, whereas above 50 GeV, it is concentrated along the equator with little polar emission. B. Arsioli & E. Orlando (2024) observed that during the solar maximum, higher-energy emission (20–150 GeV) is predominantly localized at the solar south pole, while lower-energy emission (5–20 GeV) is concentrated at the north pole. Additionally, the SSG91 model did not account for the observed spectral dip at ~30 GeV (T. Linden et al. 2018, 2022; Q.-W. Tang et al. 2018), which was not detected by the CALorimetric Electron Telescope (CALET; N. W. Cannady 2022).

Several modeling attempts have been made to go beyond the SSG91 model and explain the new observations. The solar atmospheric magnetic field is known to significantly influence gamma-ray emission (see, e.g., B. Zhou et al. 2017). In order to address the gamma-ray excess from the disk, J.-T. Li et al. (2024a) proposed a two-component model with a flux tube (an open magnetic field line) and a flux sheet. Given the scarcity of studies focusing on closed photo-/chromospheric magnetic structures, E. Puzzoni et al. (2024) investigated the impact of closed magnetic arcades on gamma-ray emission from the solar disk. Smaller-scale magnetic fluctuations were found essential to explain the <30 GeV Fermi-LAT flat spectrum. Additionally, K. C. Y. Ng et al. (2024) estimated gamma-ray production semianalytically by examining how the horizontal component of the magnetic field below the photosphere affects the flux in the TeV energy range, measured by HAWC. R. Alfaro et al. (2024) emphasized the significance of the photospheric magnetic field in modeling the GCR deficit, and the resulting excess of solar disk gamma rays, at the Sun’s

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Lunar and Planetary Laboratory, University of Arizona, 1629 E University Boulevard, Tucson, AZ 85721, USA; epuzzoni@arizona.edu</footnote>

<footnote n="2">²Observatoire de la Côte d'Azur, Laboratoire Lagrange, Bd. de l'Observatoire, CS 34229, 06304 Nice cedex 4, France</footnote>

<footnote n="3">³Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>
location in the sky. This effect, known as the Sun shadow (see, e.g., M. Amenomori et al. 2018; J. Becker Tjus et al. 2020), varies with energy, as higher-energy GCRs are less deflected by the magnetic field and are therefore more likely to be blocked by the Sun, creating the shadow. M. Gutiérrez et al. (2022) analyzed the energy-dependent shadow of the Sun measured by HAWC to model the excess of gamma rays, neutrinos, and neutrons. M. N. Mazziotta et al. (2020) used FLUKA simulations to model the yields of secondary particles, including gamma rays from the solar disk, produced by GCR interactions in the solar atmosphere.

The present study extends the analysis from E. Puzzoni et al. (2024) to a synthetic 3D model incorporating static open magnetic field lines that are increasingly distorted near the Sun’s surface. Specifically, our model captures the transition of the magnetic field near the Sun’s surface from predominantly radial to braided within the photosphere. This synthetic field is not derived from a physical driver but is designed to mimic the braiding expected from processes such as granular flows, MHD turbulence, or magnetic reconnection, as shown by observations (see, e.g., R. J. Morton et al. 2023). In this sense, the model is conceptually aligned with 3D-MHD-based reconstructions (see, e.g., H. S. Hudson et al. 2020). The simulated low-photosphere field strengths range from 100 G to 1 kG, in agreement with both MHD models (e.g., N. Shchukina & J. Trujillo Bueno 2011; M. Rempel 2014; T. del Pino Alemán et al. 2018) and with observations (e.g., S. Danilovic et al. 2010; J. O. Stenflo 2013; R. Ishikawa et al. 2021). By adopting this model, we aim to explore in greater detail the migration of GCRs from the open magnetic field lines, which channel particles into the solar atmosphere, to the nearly closed photospheric magnetic structures, an aspect that remained an open question in E. Puzzoni et al. (2024). The particle trajectory is tracked with the PLUTO code (A. Mignone et al. 2007) for distinct magnetic braiding.

The Letter is organized as follows. Section 2 outlines the magnetic field model. Section 3 details the numerical setup and analysis methodology. The results are presented in Section 4, with Section 4.1 focusing on how the GCRs interaction with ambient protons depends on the magnetic structure, and Section 4.2 examining the direction of outgoing gamma rays. Section 4.3 compares the computed gamma-ray flux with Fermi-LAT and HAWC observations. Finally, Section 5 summarizes the findings.

## 2. Theoretical Model

### 2.1. Magnetic Field Structure

The upper chromosphere/low corona magnetic field $B_0$ is assumed to be stationary on the particle motion timescale and uniform. We span the range $B_0 = 5–50$ Gauss along the z-direction (see Figure 1 at $z \gtrsim 0.02 R_\odot$ up to the corona), consistent with recent spectropolarimetric measurements from the Daniel K. Inouye Solar Telescope (DKIST; see T. A. Schad et al. 2024) and with 3D-MHD resistive simulations of jets adjacent to arcades or more generally closed structures (J. J. González-Avilés et al. 2018). These values are also in agreement with the magnetic field strength at the top of the magnetic arcade as assumed in E. Puzzoni et al. (2024) and the mean field in H. S. Hudson et al. (2020). Such values of $B_0 = 5–50$ G in the low corona correspond to a photospheric magnetic field strength of

<figure><img src="image_4.png" /><figcaption>Figure 1. Magnetic field magnitude in the 3D domain (with logarithmic color scale in Gauss) together with the decomposition of the particle velocity vector defined by the angles $\theta$ and $\phi$. The xy-plane is tangent to the solar surface, and the z-axis is parallel to the local radial direction.</figcaption></figure>

$B_{max} \sim 100 – 1000$ G, consistent with MHD simulations as well as remote-sensing observations (see, e.g., S. Danilovic et al. 2010; J. O. Stenflo 2013).

As GCRs approach the photosphere, the magnetic field acquires a 3D structure. Therefore, we use a synthetic 3D total magnetic field to mimic the complex 3D photospheric field that can be expressed as

$$ \mathbf{B}(x, y, z) = B_0 \hat{k} + \nabla \times (f(z) \mathbf{S}(x, y, z)), \quad (1) $$

satisfying the divergence-free condition $\nabla \cdot \mathbf{B}(x, y, z) = 0$.

Expanding Equation (1), we obtain

$$ \mathbf{B}(x, y, z) = B_0 \hat{k} + f(z) \nabla \times \mathbf{S}(x, y, z) + \nabla f(z) \times \mathbf{S}(x, y, z), \quad (2) $$

where $f(z)$ is an arbitrary scalar function modulating the onset of the transverse components and $\mathbf{S}(x, y, z)$ represents the vector potential associated with the synthetic magnetic field $\delta\mathbf{B}(x, y, z)$, i.e., $\delta\mathbf{B}(x, y, z) = \nabla \times \mathbf{S}(x, y, z)$, derived in J. Giacalone (2021) as

$$ \mathbf{S}(x, y, z) = \sum_{n=1}^{N_m} \frac{A_n}{k_n} (\sin \alpha_n \hat{x}_n' + i \cos \alpha_n \hat{y}_n') e^{ik_n z_n' + i\beta_n}. \quad (3) $$

In our model, we take $f(z) = e^{-(z/\Lambda)^2}$ (with $\Lambda \simeq 1.41 \times 10^{-2} R_\odot$ being the scale height). The total magnetic field can thus be rewritten as

$$ \mathbf{B}(x, y, z) = B_0 \hat{k} + f(z) \delta \mathbf{B}(x, y, z) + \nabla f(z) \times \mathbf{S}(x, y, z). \quad (4) $$

The spatial variation of $|\mathbf{B}|$ in the 3D domain of size $L^3$, where $L = 0.06 R_\odot$, is depicted in Figure 1: the magnetic field remains relatively uniform except for the low chromosphere and below, where distortions build up. The plane $z = 0$
represents the base of the photosphere, i.e., the surface where the vertical optical depth for radiation at a wavelength of 500 nm reaches unity (R. G. Athay 1976). We note that $S(x, y, z)$ and $\delta B(x, y, z)$ have nonzero components along all three directions, while the third term on the right-hand side of Equation (4) is zero along the z-axis ($\partial_x f(z) = \partial_y f(z) = 0$). As a result, no component of the gradient-curvature drift velocity vanishes.

$\delta B(x, y, z)$ is defined in J. Giacalone & J. R. Jokipii (1999) as

$$ \delta B(x, y, z) = \sum_{n=1}^{N_m} A_n (\cos \alpha_n \hat{x}_n' + i \sin \alpha_n \hat{y}_n') e^{ik_n z_n' + i\beta_n}, \quad (5) $$

with

$$ \begin{bmatrix} x' \\ y' \\ z' \end{bmatrix} = \begin{bmatrix} \cos \tilde{\theta}_n \cos \tilde{\phi}_n & \cos \tilde{\theta}_n \sin \tilde{\phi}_n & -\sin \tilde{\theta}_n \\ -\sin \tilde{\phi}_n & \cos \tilde{\phi}_n & 0 \\ \sin \tilde{\theta}_n \cos \tilde{\phi}_n & \sin \tilde{\theta}_n \sin \tilde{\phi}_n & \cos \tilde{\theta}_n \end{bmatrix} \begin{bmatrix} x \\ y \\ z \end{bmatrix}, \quad (6) $$

by summing over a large number of plane waves with wave mode $n$ (for a total of $N_m = 556$ modes), each characterized by amplitude $A(k_n)$, wavenumber $k_n$, randomly oriented polarization $\alpha_n$, and phase $\beta_n$. The angles $\tilde{\theta}_n$ and $\tilde{\phi}_n$ represent the random direction of propagation of each mode along the $z'$-axis of the local (primed) reference frame. The random parameters are reshuffled every 100 particles, leading to 100 different realizations of the total magnetic field. We adopt

$$ A^2(k_n) = \sigma'^2 G(k_n) \left[ \sum_{n=1}^{N_m} G(k_n) \right]^{-1}, \quad (7) $$

where

$$ G(k_n) = \frac{\Delta V_n}{1 + (k_n L_s)^q} \quad (8) $$

with $q = 11/3$. Here, $\sigma'^2 = \langle \delta B^2(x, y, z) \rangle$, where $\langle \cdot \rangle$ indicates the average over the different magnetic field realizations, corresponds to the relative energy density of the magnetic components added to $B_0$; $\Delta V_n = 4\pi k_n^2 \Delta k_n$ is the normalization factor dependent on the logarithmic spacing $\Delta k_n$ between $k_n$ (meaning that $\Delta k_n/k_n$ is a constant); and $L_s$ represents the largest horizontal scale of the magnetic photospheric structures.

In this work, we set $L_s = 0.003 R_\odot \simeq 2$ Mm, comparable with the granular scale ($\sim 1$ Mm; see, e.g., J. Trujillo Bueno et al. 2004). $S(x, y, z)$ in Equation (3) spans a wavenumber interval from $L_{\min} = 10^{-3} R_\odot$ to $L$, corresponding to wave-numbers $k_{\max} = 2\pi/L_{\min}$ and $k_{\min} = 2\pi/L$, respectively. The assumed range in scales ensures resonant interactions across all GCR energies. We explore different values of $\sigma^2 = \sigma'^2/B_0^2 = 0.1, 1, 10$.

The upper plot of Figure 2 shows the average angle $\langle \theta \rangle_\perp = \arctan(\sqrt{B_x^2(x, y, z) + B_y^2(x, y, z)/B_0})$ between the magnetic field and the vertical direction $\hat{z}$ as a function of $z$, illustrating the deviation of the magnetic field from the vertical as it approaches the photosphere base. This definition of $\langle \theta \rangle_\perp$ is clearly an upper limit for $\sigma^2 \gg 1$, which includes longitudinal compressions. The $\langle \theta \rangle_\perp$ monotonically increases downward with a peak value for $\sigma^2 = 10$ at $50^\circ-60^\circ$ from the

vertical or $\sqrt{B_x^2 + B_y^2}/B_0 = 1.2-1.7$. This predominance of the horizontal component of the photospheric field is confirmed by high-resolution measurements, which report an average value $\sim 2.0-5.6$ (O. Steiner et al. 2008). As a further clarification, this distortion is not the result of turbulent cascade or wave-wave collisions and does not contradict observations that identify maximal turbulent fluctuations at an intermediate height within the photosphere (see, e.g., T. A. Schad et al. 2024). The general trend of vertical magnetic field dominance in the chromosphere and low corona, contrasted with horizontal field dominance in the photosphere, is also observed by S. Danilovic et al. (2010) and J. O. Stenflo (2013). We note that our synthetic model cannot reproduce the peak in $\langle \theta \rangle_\perp$ due to the small-scale dynamo in the photosphere found numerically by M. Rempel (2014), as the peak cannot be resolved by the grid resolution for the values of $B_0$ adopted here, comparable with the vertical field therein. The magnetic field lines and magnitude are shown in the lower left ($xz$-plane) and right ($xy$-plane) panels of Figure 2 for a single realization of the total magnetic field with $B_0 = 5$ G and $\sigma^2 = 10$. Note that $L_s < \Lambda$ to ensure that the largest magnetic structures are contained in the scale height $\Lambda$. The increasing distortion of the laminar field lines in the lower left panel is visible below $z < L/2 = 0.03 R_\odot$, corresponding to the height of the closed arcade in E. Puzzoni et al. (2024). Small-scale magnetic structures form at $z \lesssim 0.01 R_\odot$, consistent with the heights of similar structures reported by H. S. Hudson et al. (2020). The lower right panel shows magnetic structures near the bottom of the photosphere, with sizes and a maximum magnetic field strength of $\sim 100$ Gauss with a large filling factor that aligns with the values inferred by combining radiative transfer simulations with polarization observations (J. Trujillo Bueno et al. 2004), which ruled out earlier estimates of the magnetic field at least 2 orders of magnitude smaller in the photosphere (see, e.g., M. Faurobert-Scholl et al. 1995; M. Faurobert et al. 2001). Here, $B_0 = 50$ G with $\sigma^2 = 10$ leads to a maximum magnetic field strength of $\sim 1$ kG, consistent with observations (see, e.g., J. O. Stenflo 2013).

The synthetic field $\delta B(x, y, z)$ is periodic in both the x- and y-directions. The wavenumbers $k_x$ and $k_y$ are selected to be close to randomly chosen values but must be multiples of $L$. It is essential that these new wavenumbers maintain consistency with the specified $k$: if the new values of $k_x$ and $k_y$ lead to $k_x^2 + k_y^2 > k^2$, the values must be reselected, using a new random number.

## 2.2. Density Profile

Following E. Puzzoni et al. (2024), we incorporate the density profile from 1D hydrostatic models (J. M. Fontenla et al. 1993) for the low photosphere ($\rho_0 \sim 10^{-6}$ g cm$^{-3}$ at $z=0$) up to 4 Mm (R. J. Morton et al. 2023) and from J. J. González-Avilés et al. (2021) up to 30 Mm. The motion of GCRs in the total magnetic field $B(x, y, z)$ in Equation (1) is computed by solving the Lorentz equation. In the upper chromosphere, the uniform magnetic field $|B(x, y, z)| = B_0 = 5$ G, with a proton density $\rho \sim 10^{-16}$ g cm$^{-3}$ (from J. J. González-Avilés et al. 2021), results in a typical Alfvén speed $v_A = B_0/\sqrt{4\pi\rho} \simeq 1.4 \times 10^8$ cm s$^{-1}$. The ratio between this Alfvén velocity to the particle speed ($v \sim c$) is $v_A/c \simeq 4.7 \times 10^{-3}$, which validates the magnetostatic assumption also for the high fields and very dense low photosphere.
<figure><img src="image_5.png" /><figcaption>Figure 2. Top panel: average angle ($\theta$)$_{\perp}$ (in degrees) as a function of the altitude z with $\sigma^2 = 10$. Bottom panels: magnetic field lines (in black) and magnitude (color bar) in the xz-plane (lowest y = 0.0003 $R_{\odot}$ plane in the grid, left panel) and xy-plane (lowest z = 0.0003 $R_{\odot}$ plane in the grid, right panel) with $\sigma^2 = 10$.</figcaption></figure>

### 3. Numerical Setup

#### 3.1. Computational Domain

Test-particle proton 3D numerical simulations are conducted with the PLUTO code (A. Mignone et al. 2007, 2012). We adopt a 3D computational domain with dimensions $L^3$. In particular, the domain is defined within $-0.03 < x/R_{\odot} < 0.03$, $0 < y/R_{\odot} < 0.06$, and $0 < z/R_{\odot} < 0.06$. The 2D resolution of $1000^2$ used in E. Puzzoni et al. (2024) is extended to three dimensions in this work, employing a grid of $100^3$ cells with a spatial resolution of $6 \times 10^{-4} R_{\odot} \simeq 418$ km. As explained in the previous section, periodic boundary conditions are employed in the x- and y-directions. Our simulations show that the GCR excursion out of the box sides before either GCR-proton interaction or GCR escape both downward into the Sun and upward toward the corona is only a few times L; thus, $\simeq 0.1-0.2 R_{\odot}$, supporting the approximation of the Sun surface with its tangent xy-plane. Outflow boundary conditions are employed along the z-direction for particles precipitating into the Sun surface or leaving the chromosphere with no interaction.

#### 3.2. Particle Initial Conditions and Motion in the Magnetic Field

GCR protons are initialized with one particle per cell at two distinct heights: in the laminar (i.e., $z = z_{up} = 0.05 R_{\odot} \sim 35,000$ km) and braided field region ($z = z_{down} = 0.005 R_{\odot} \sim 35,000$ km).

In each simulation run, a total of $N_{inj} = 10,000$ protons are injected. GCRs are injected at both heights with an isotropic velocity distribution as in E. Puzzoni et al. (2024), supported by the observed isotropy of the anomalous cosmic rays (ACRs) up to a few MeV as close as $\sim 0.1$ au (J. S. Rankin 2024).

The particle orbits are integrated using the Boris algorithm, which is already implemented in the PLUTO code (A. Mignone et al. 2007, 2012). The total magnetic field is interpolated (C. K. Birdsall & A. B. Langdon 1991) at the particle’s position (see A. Mignone et al. 2018 for further details on the interpolation). Orbit integration continues until $t_{stop} = 25 R_{\odot}/c \sim 58$ s, which corresponds to the time when all GCRs have either interacted within the computational domain or escaped from it.

#### 3.3. Interacting Protons and Gamma-Ray Flux

As in E. Puzzoni et al. (2024), the interaction time $t_{int}$ of the GCR proton-solar atmosphere proton is taken into account by

$$ \Delta t/t_{\text{int}} = \int_{t_i}^{t_f} dt/t_{\text{int}}(z) = \sigma_{pp} v \int_{t_i}^{t_f} n[z(t)] dt, \quad (9) $$

where $\Delta t$ is the total elapsed time from $t=0$ up to the final time $t_f$, $dt$ is the time step of the code, $v \simeq c$ is the typical GCR speed, and $n[z(t)]$ represents the number density in the solar atmosphere of single proton (in the nearly fully ionized chromosphere) or neutral hydrogen (in the photosphere). The
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left panel: ratio $N_{int}/N_{inj}$ ($N_{esc}/N_{inj}$) for the interacting (escaping) GCRs marked by the solid (dashed) lines as a function of $\sigma^2$ for injection at $z = z_{up}$ and $z = z_{down}$ for 100 GeV (red and magenta lines, respectively) and 10 TeV (blue and black lines, respectively) GCRs. Right panel: same as Figure 2, lower left panel, with color scale showing the gyroradius of 10 TeV GCRs with magnetic field lines (in white).</figcaption></figure>

total inelastic cross section of the proton–proton interaction $\sigma_{pp}$ is obtained from Equation (1) of E. Kafexhiu et al. (2014). We are interested in photons produced within the computational domain ($\Delta t/t_{int} > 1$) and observed at Earth ($v_z > 0$). Indeed, gamma rays are assumed to be emitted along the same direction as the GCRs at $t = t_{int}$ due to momentum conservation (see Section 4.2).

The typical photon energy is $E_\gamma \sim 0.1 E_p$ via the $\pi^0$ meson primary channel (S. R. Kelner et al. 2006) for proton–proton interactions (i.e., $p + p \to \pi^0 \to \gamma + \gamma$; see J. Dörner et al. 2025 for a comparison with empirical hadronic interaction models). Following the method outlined in E. Puzzoni et al. (2024) and J.-T. Li et al. (2024a), we calculate the gamma-ray flux:

$$ \Phi_{\gamma}^{*}(E_{\gamma}) = \frac{dN_{\gamma}}{dE_{\gamma}} = cn_{p} \int_{E_{\gamma}}^{\infty} \sigma_{pp}(E_{p})J_{p}(E_{p})F_{\gamma} \\ \times \left(\frac{E_{\gamma}}{E_{p}}, E_{p}\right) \frac{dE_{p}}{E_{p}}, \qquad (10) $$

where $n_p$ is the proton number density integrated along the z-direction within the dense atmosphere; $F_\gamma(E_\gamma/E_p, E_p)$ is an empirical function accounting for the flux of outgoing $E_\gamma$ photon from an impinging $E_p$ proton (S. R. Kelner et al. 2006, Equation (58) therein); and $J_p(E_p) = N_{int}(E_p)/N_{inj} \times (E_p^2 dN/dE_p)_{obs} \times 1/(E_p^2 c)$, where $N_{int}/N_{inj}$ represents the ratio of interacting to injected protons, averaged over the injection altitudes $z_{up}$ and $z_{down}$. Here $(E_p^2 dN/dE_p)_{obs}$ denotes the observed, approximately isotropic, GCR proton energy spectrum from the Alpha Magnetic Spectrometer (AMS; M. Aguilar et al. 2021) for energies $\lesssim 1$ TeV and from the Cosmic Ray Energetics And Mass mission

for the International Space Station (ISS-CREAM) Experiment (G. H. Choi et al. 2022) at higher energies. The final gamma-ray flux expression is $\Phi_\gamma(E_\gamma) = 2\pi R_\odot^2/L^2 \times (R_\odot/R_{\text{au}})^2 \times \Phi_\gamma^*(E_\gamma)$, where $R_{\text{au}}$ is the distance between the Sun and Earth.

# 4. Results

## 4.1. Interacting versus Escaping GCRs

The left panel of Figure 3 shows $N_{int}$ (solid lines) and $N_{esc}$ (dashed lines) divided by $N_{inj}$ as a function of $\sigma^2$. The red (blue) and magenta (black) lines show the results for injection at $z = z_{up}$ and $z = z_{down}$, respectively, for 100 GeV (10 TeV) GCRs. Due to the periodic boundary conditions in the xy-plane, protons can escape only along the z-direction (i.e., move outward in the corona or precipitate into the Sun surface). The panel shows that when particles are injected at $z = z_{down}$, $N_{int}$ increases with $\sigma^2$, as a result of trapping within smaller-scale magnetic structures. Consequently, $N_{esc}$ decreases (see the corresponding dashed lines). For protons injected at $z = z_{up}$, $N_{int}$ remains relatively small and approximately constant across all values of $\sigma^2$. In this scenario, most GCRs escape the computational domain before interacting, as $N_{esc}/N_{inj} \sim 1$, freely streaming out from the top of the box due to positive vertical velocities ($v_z > 0$, in half of GCRs at injection) along the laminar magnetic field due to the large $r_g$ (see right panel of Figure 3).

In this scenario, field distortions in the xy-plane increase the residence time of GCRs in the dense photo- and chromosphere. This condition enhances trapping and gamma-ray emission for higher $\sigma^2$. In E. Puzzoni et al. (2024), a large-scale, nonuniform unperturbed field $B_0(x)$ is overlapped with ubiquitous small-scale fluctuations: $\lambda_\parallel$ is finite everywhere and scales inversely with fluctuation power (F. Fraschetti et al. 2022). In that magnetic field geometry, a larger cross-field
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left panel: angular components $\theta$ and $\phi$ of the velocity vector of the GCRs at $t_{\text{int}}$ for $E_p = 100$ GeV with $\sigma^2 = 10$ and injection at $z = z_{\text{down}}$. Right panel: same as left panel for $E_p = 10$ TeV.</figcaption></figure>

diffusion (F. Fraschetti & J. R. Jokipii 2011) enables horizontal migration to the arcades, and a smaller $\lambda_\|$ favors the GCR-proton collision (higher $N_{\text{int}}$).

## 4.2. Angular Distribution of Outgoing Gamma Rays

Figure 4 shows the $0^\circ < \theta < 90^\circ$ and $-180^\circ < \phi < 180^\circ$ angles of the GCR velocity vector at time $t_{\text{int}}$ for the 100 GeV (left panel) and 10 TeV (right panel) cases for $\sigma^2 = 10$ and for injection at $z = z_{\text{down}}$. Due to momentum conservation, the direction of the outgoing GCR is nearly aligned with the direction of the outgoing photons in the solar plasma frame. Z. Li et al. (2024b) demonstrated that the angular distribution broadens for $E_p \lesssim 100$ GeV (or $E_\gamma \lesssim 10$ GeV); such an effect is neglected herein as it introduces only a correction smaller by a factor $\sim 4-5$ for $E_p < 10$ GeV and even smaller at $E_p > 10$ GeV.

Figure 4 also shows, for $E_p = 10$ TeV (right panel), that most GCRs interact across the entire $\phi$ range, within the narrow range $70^\circ \lesssim \theta \lesssim 90^\circ$ nearly aligned with the tangent plane (i.e., xy-plane; see Figure 1), whereas for $E_p = 100$ GeV (left panel), the $\theta$ range broadens down to $\theta \approx 50^\circ$. This implies that as GCR energy increases, gamma rays are emitted close to the plane tangent to the Sun surface, as also found in E. Puzzoni et al. (2024). A similar anisotropic gamma-ray emission pattern on the solar disk is also suggested by Fermi-LAT observations (T. Linden et al. 2022).

## 4.3. Gamma-Ray Flux

Figure 5 presents the solar disk gamma-ray flux derived from our simulations for $B_0 = 5$ G and different values of $\sigma^2$: 0.1 (green), 1 (brown), and 10 (orange diamonds, leading to a low-photosphere $B_{\text{max}} \sim 100$ G). The orange stars represent results obtained with $\sigma^2 = 10$ and $B_0 = 50$ G, leading to $B_{\text{max}} \sim 1$ kG. These results are compared with the Fermi-LAT data (T. Linden et al. 2022) with and without the solar flares (purple and black points, respectively), with the HAWC observations (in blue; A. Albert et al. 2023) and the estimated gamma-ray flux at 1 au assuming 100% efficiency in converting protons to gamma rays (solid black line). The latter is calculated using the expression $\phi_\gamma(1 \text{ au}) \approx (E_p^2 dN/dE_p)_{\text{obs}}(R_\odot/R_{\text{au}})^2(E_\gamma/E_p)$,

where $E_\gamma/E_p \sim 1/10$. Results at energies $E_\gamma < 10$ GeV (corresponding to $E_p < 100$ GeV) are omitted, as the gyromotion of particles in this energy range cannot be resolved given the spatial resolution of $\sim 418$ km.

The gamma-ray flux increases with $\sigma^2$ across the entire energy range analyzed due to the increase in $N_{\text{int}}$ (see left panel of Figure 3). Small-scale magnetic structures in the low photosphere effectively trap GCRs, extending their residence time and consequently increasing the probability of interaction. This finding aligns with expectations and our earlier work (E. Puzzoni et al. 2024). A notable difference from the arcade trapping of GCRs analyzed in E. Puzzoni et al. (2024) is that, at higher energies ($E_\gamma \gtrsim 100$ GeV), the gamma-ray flux from an arcade does not depend on the amplitude of turbulent fluctuations. That result indicates that the gamma-ray flux is due to the migration from open to closed field lines: once a GCR is trapped in the arcade, the flux is independent of the turbulence strength and drops as a result of the drop of the incoming GCRs flux unable to enter the arcade.

With the specific magnetic field geometry in the present Letter, we observe a novel trend in the “dip” energy range (from $\sim 30$ up to 100 GeV) as $\sigma^2$ increases, as highlighted in the zoomed-in inset at the top right of Figure 5. For $\sigma^2 \lesssim 0.1$, the gamma-ray flux decreases in this region, at a rate comparable to the maximum efficiency spectrum (solid black line). As $\sigma^2$ increases, however, the gamma-ray flux begins to flatten within this range, even for $B_0$ 10 times stronger. As elaborated in Section 5, this effect is not numerical but genuinely physical.

# 5. Discussion and Conclusions

This Letter extends the previous work on the role of closed magnetic arcades in the solar disk gamma-ray emission presented in E. Puzzoni et al. (2024) by using a new model for the magnetic field near the solar surface, which has the property of being open to space at the top and increasingly horizontal near the surface. In our new model, the GCRs follow the open field lines from the top toward the solar surface where they encounter an irregular field (see J. Giacalone 2021).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Solar disk average gamma-ray spectrum observed by Fermi-LAT with (purple points) and without (black points) solar flares, from T. Linden et al. (2022), by HAWC (blue line, from A. Albert et al. 2023), and obtained from our simulations, with $\sigma^2 = 10$ (orange diamonds), $\sigma^2 = 1$ (in brown), and $\sigma^2 = 0.1$ (in green). The orange stars represent simulations with $\sigma^2 = 10$ but with a total magnetic field strength 10 times higher ($B_0 = 50$ G leading to low-photosphere $B_{max} \sim 1$ kG). A detailed comparison is shown in the zoomed-in inset at the top right. The solid black line represents the estimated (see text for further details) gamma-ray flux considering a 100% efficiency of conversion from GCRs to gamma rays.</figcaption></figure>

In agreement with J.-T. Li et al. (2024a) and E. Puzzoni et al. (2024), GCR protons predominantly interact with solar atmospheric protons in the lower dense solar photosphere (see, e.g., R. J. Morton et al. 2023). $N_{\text{int}}$ increases as expected with the field irregularity ($\sigma^2$) if GCRs are injected within the nonvertical field layer. However, $N_{\text{int}}$ is approximately independent on $\sigma^2$ for injection farther up in the laminar field region. In this case, the majority of high-energy GCRs ($E_p \gtrsim 10$ TeV) escape upward into the lower corona due to their large $r_g$ (as sideway escape is treated via periodic boundary conditions in the xy-plane, unlike in E. Puzzoni et al. 2024). To clarify this, we injected GCRs at lower altitudes while remaining within the laminar field region (specifically, at $z = 0.025 R_\odot$ instead of $z = 0.05 R_\odot$). This adjustment resulted in an increase in $N_{\text{int}}$ though it had a minimal impact on the intensity of the gamma-ray flux.

It is worthwhile to note that our simulations are consistent with the GCR pitch-angle isotropy (for ACRs, see J. S. Rankin 2024): if initially only $v_z < 0$ is allowed, the resulting gamma-ray flux would be up to 5 times higher than our calculated low-energy flux, while the spectral slope would remain approximately unchanged. An isotropic distribution is also assumed in a study focusing on TeV energy range only (K. C. Y. Ng et al. 2024).

Gamma-ray emission spans the full range in the $\phi$ angular component of the outgoing GCR at all energies (see Figure 4). At higher energies, a gamma-ray “trapping cone”, i.e., angular range from the vertical to the Sun surface where proton escape is prohibited, extends down to $\theta \approx 90^\circ$ (along the xy-plane); thus, gamma-ray emission is squeezed to the tangent plane. At lower energy, the trapping cone is narrower, gamma rays are emitted at smaller $\theta$, and the solar disk emission appears more

isotropic. This energy-dependent emission pattern on the solar disk is roughly consistent with Fermi-LAT observations (see T. Linden et al. 2022) and numerical results reported in E. Puzzoni et al. (2024).

The new synthetic magnetic field model also reflects a realistic solar atmosphere, characterized by nearly closed magnetic structures extending to altitudes on the order of Mm, as observed by Solar Orbiter (P. Antolin et al. 2023), with maximum magnetic field strengths ranging from 100 G up to 1 kG and profile of the horizontal field both consistent with remote-sensing observations (see, e.g., S. Danilovic et al. 2010; J. O. Stenflo 2013) and MHD simulations (see, e.g., M. Rempel 2014 and references therein). Rather than deriving the field from a specific physical mechanism, we constructed it to mimic a geometry consistent with observations and likely produced by processes such as granular flows, MHD turbulence, or reconnection-driven braiding, which are known to distort field lines in the low solar atmosphere (see, e.g., R. J. Morton et al. 2023). The model is consistent with 3D MHD-based reconstructions (see, e.g., J. J. González-Avilés et al. 2018; H. S. Hudson et al. 2020).

A highly irregular photospheric field was essential for revealing features of the gamma-ray flux within the dip region. In E. Puzzoni et al. (2024), we argued that the dip/rebrightening of the gamma-ray flux between approximately 30 and 100 GeV can be due to an energy cutoff in the migration of the GCRs into magnetic arcades (for a given turbulence strength, only GCRs with large enough $E_p$ can migrate into the arcades and emit gamma rays). Here, we show that as the field becomes more irregular, the gamma-ray flux flattens between 30 and 100 GeV. It is noteworthy that this flattening saturates at the shape shown in Figure 5, rather than
forming a local maximum for $\sigma^2 \gtrsim 10$. This has been demonstrated by running the $\sigma^2 = 50$ case as well as the case with a stronger magnetic field ($B_{max} \sim 1$ kG reached with $B_0 = 50$ G). Notably, this spectral dip was not captured in earlier modeling attempts (see, e.g., M. Gutiérrez et al. 2022; J.-T. Li et al. 2024a, with the single-structure model).

We argue that the dip/rebrightening is associated with the condition $r_g(300 \text{ GeV}) \simeq L_s$ in a $B_0 = 5$ G field (see Figure 2). To investigate this, we break such a condition by increasing $L_s$ to $3L_s \sim \Lambda$, with unchanged $\Lambda$ and $\sigma^2$, thereby increasing the length scale that maintains coherent the horizontal field. Horizontal GCR motion for a longer time reduces the escape along the vertical direction, leading to an enhanced flux. Our findings show that for $\sigma^2 = 0.1$, the slope of the gamma-ray flux in the dip region remains unchanged, while its intensity increases. The flux does not change significantly for higher $\sigma^2$ (or higher $\langle\theta\rangle_\perp$). GCR-proton interactions (in the dip energy range) have a mean-free path $\lambda_{pp} \simeq (\sigma_{pp} \rho_0/m_p)^{-1} \lesssim L_s$ for photospheric densities $\rho_0 \sim 10^{-7}-10^{-6} \text{ g cm}^{-3}$ (see Section 2.2); thus, the flux at the dip is unaffected by a greater $L_s$.

In addition, we have examined the effect of $\Lambda$ on the gamma-ray flux. Doubling $\Lambda$ for all the values of $\sigma^2$ with an unchanged density vertical profile does alter neither the intensity nor the slope of the resulting gamma-ray spectrum at $E_\gamma \gtrsim 30$ GeV. The reason is that, in our synthetic field, modifying the height of deformation of the vertical field $\Lambda$ does not affect the underlying density profile: the denser region where the majority of GCR-proton collisions take place remains unchanged, as well as the gamma-ray production. Therefore, we conclude that the dip/rebrightening is not a numerical artifact resulting from the chosen computational setup but rather arises from a physical effect associated with the small-scale magnetic field structures.

The present model has limitations: we cannot include closed magnetic field lines, which efficiently trap GCRs enhancing the gamma-ray flux, and the time dependence of the latitudinal emission (equatorial versus polar) over the solar cycle; in addition, smaller-scale structures that contribute low-energy flux ($\lesssim 10$ GeV) are not resolved. It is important to underline that this model provides a geometric framework, consistent with observed/simulated irregularity, that explains the gamma-ray production; the magnetic field structure is not evolved via a specific MHD model. Identifying which physical processes are responsible for the braiding remains an open question for future work. We also assume that the field is static over the GCR propagation timescale (tens of seconds), which is much shorter than the characteristic timescale of quiescent photo- and chromospheric field evolution (e.g., tens of minutes for the active regions loops, F. Reale 2014). This justifies neglecting temporal variations in the field while focusing on its spatial structure. With the inclusion of open magnetic field lines in our modeling, a natural progression is to investigate the effect of interchange reconnection, which could be addressed in future work by including a multiscale field model with both open and closed topologies, informed by high-resolution MHD simulations. Despite these limitations, this work supports the findings of our previous work (E. Puzzoni et al. 2024), highlighting the critical role of the photospheric magnetic field in shaping the gamma-ray flux from the solar disk. Our inclusion of an irregular field provides a crucial new perspective, helping to address previously unexplained

spectral features, such as the dip observed in the $\sim 30-100$ GeV range, which had remained unexplained in earlier models.

## Acknowledgments

The authors are grateful to the referee for very constructive feedback. The authors express gratitude for the insightful feedback from Drs. M. Asgari-Targhi, J. Beacom, S. S. Cerri, O. Cohen, C. Fong, H. S. Hudson, J.-T. Li, T. Linden, K. C. Y. Ng, and A. Peter. This work was partially supported by NASA under grant 80NSSC22K0040. F.F. was partially supported by NASA under grants 80NSSC21K0119 and 80NSSC21K1766. E.P., J.K., and J.G. were partially supported by NASA under grant 18-DRIVE18_2-0029. J.G. acknowledges partial support from NSF under grant No. 1931252.

*Facilities:* Fermi, HAWC.

*Software:* PLUTO code (A. Mignone et al. 2007, 2012).

## ORCID iDs

Eleonora Puzzoni https://orcid.org/0009-0009-5314-348X

Federico Fraschetti https://orcid.org/0000-0002-5456-4771

József Kóta https://orcid.org/0000-0002-3715-0358

Joe Giacalone https://orcid.org/0000-0002-0850-4233

## References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJ, 736, 131

Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2013, ApH, 50, 26

Achariya, A., Adelfio, A., Ajello, M., et al. 2025, arXiv:2505.06348

Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2021, PhR, 894, 1

Albert, A., Alfaro, R., Alvarez, C., et al. 2023, PhRvL, 131, 051201

Alfaro, R., Alvarez, C., Arteaga-Velázquez, J. C., et al. 2024, ApJ, 966, 67

Amenomori, M., Bi, X. J., Chen, D., et al. 2018, PhRvL, 120, 031101

Antolin, P., Dolliou, A., Auchère, F., et al. 2023, A&A, 676, A112

Arsioli, B., & Orlando, E. 2024, ApJ, 962, 52

Athay, R. G. 1976, The Solar Chromosphere and Corona: Quiet Sun (Dordrecht: Reidel)

Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071

Becker Tjus, J., Desiati, P., Döpper, N., et al. 2020, A&A, 633, A83

Birdsall, C. K., & Langdon, A. B. 1991, Plasma Physics via Computer Simulation (Bristol: Adam Hilger)

Cannady, N. W. 2022, ICRC (Berlin), 37, 604

Choi, G. H., Seo, E. S., Aggarwal, S., et al. 2022, ApJ, 940, 107

Danilovic, S., Schüssler, M., & Solanki, S. K. 2010, A&A, 513, A1

del Pino Alemán, T., Trujillo Bueno, J., Štěpán, J., & Shchukina, N. 2018, ApJ, 863, 164

Dörner, J., Morejon, L., Kampert, K.-H., & Becker Tjus, J. 2025, JCAP, 2025, 043

Faurobert, M., Arnaud, J., Vigneau, J., & Frisch, H. 2001, A&A, 378, 627

Faurobert-Scholl, M., Feautrier, N., Machefert, F., Petrovay, K., & Spielfiedel, A. 1995, A&A, 298, 289

Fontenla, J. M., Avrett, E. H., & Loeser, R. 1993, ApJ, 406, 319

Fraschetti, F., Alvarado-Gómez, J. D., Drake, J. J., Cohen, O., & Garraffo, C. 2022, ApJ, 937, 126

Fraschetti, F., & Jokipii, J. R. 2011, ApJ, 734, 83

Giacalone, J. 2021, ApJ, 912, 83

Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204

González-Avilés, J. J., Guzmán, F. S., Fedun, V., et al. 2018, ApJ, 856, 176

González-Avilés, J. J., Murawski, K., Srivastava, A. K., Zaqarashvili, T. V., & González-Esparza, J. A. 2021, MNRAS, 505, 50

Gutiérrez, M., Masip, M., & Muñoz, S. 2022, ApJ, 941, 86

Hudson, H. S., MacKinnon, A., Szydlarski, M., & Carlsson, M. 2020, MNRAS, 491, 4852

Ishikawa, R., Trujillo Bueno, J., del Pino Alemán, T., et al. 2021, SciA, 7, eabe8406

Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, PhRvD, 90, 123014

Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, PhRvD, 74, 034018

Li, J.-T., Beacom, J. F., Griffith, S., & Peter, A. H. G. 2024a, ApJ, 961, 167

Li, Z., Ng, K. C. Y., Chen, S., Nan, Y., & He, H. 2024b, ChPhC, 48, 045101

Linden, T., Beacom, J. F., Peter, A. H. G., et al. 2022, PhRvD, 105, 063013

Linden, T., Zhou, B., Beacom, J. F., et al. 2018, PhRvL, 121, 131103
Mazziotta, M. N., Luque, P. D. L. T., Di Venere, L., et al. 2020, PhRvD, 101, 083011

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228

Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13

Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7

Morton, R. J., Sharma, R., Tajfirouze, E., & Miriyala, H. 2023, RvMPP, 7, 17

Ng, K. C. Y., Beacom, J. F., Peter, A. H. G., & Rott, C. 2016, PhRvD, 94, 023004

Ng, K. C. Y., Hillier, A., & Ando, S. 2024, arXiv:2405.17549

Puzzoni, E., Fraschetti, F., Kóta, J., & Giacalone, J. 2024, ApJ, 973, 118

Rankin, J. S. 2024, in AGU Fall Meeting 2024, What's Next for Science (Washington, D.C.: AGU), SH21C-2832

Reale, F. 2014, LRSP, 11, 4

Rempel, M. 2014, ApJ, 789, 132

Schad, T. A., Petrie, G., Kuhn, J., et al. 2024, SciA, 10, eadq1604

Seckel, D., Stanev, T., & Gaisser, T. K. 1991, ApJ, 382, 652

Shchukina, N., & Trujillo Bueno, J. 2011, ApJL, 731, L21

Steiner, O., Rezaei, R., Schaffenberger, W., & Wedemeyer-Böhm, S. 2008, ApJL, 680, L85

Stenflo, J. O. 2013, A&A, 555, A132

Tang, Q.-W., Ng, K. C. Y., Linden, T., et al. 2018, PhRvD, 98, 063019

Trujillo Bueno, J., Shchukina, N., & Asensio Ramos, A. 2004, Natur, 430, 326

Zhou, B., Ng, K. C. Y., Beacom, J. F., & Peter, A. H. G. 2017, PhRvD, 96, 023015