Small-scale Magnetic Fields Are Critical to Shaping Solar Gamma-Ray Emission

The Sun is a bright gamma-ray source due to hadronic cosmic-ray interactions with solar gas. While it is known that incoming cosmic rays must generally first be reflected by solar magnetic fields to produce outgoing gamma rays, theoretical models have yet to reproduce the observed spectra. We introduce a simplified model of the solar magnetic fields that captures the main elements relevant to gamma-ray production. These are a flux tube, representing the network elements, and a flux sheet, representing the intergranular sheets. Both the tube and sheet have a horizontal size of order 100 km and serve as sites where cosmic rays are reflected and gamma rays are produced. While our simplified double-structure model does not capture all the complexities of the solar-surface magnetic fields, such as Alfvén turbulence from wave interactions or magnetic fluctuations from convection motions, it improves on previous models by reasonably producing both the hard spectrum seen by Fermi Large Area Telescope at 1–200 GeV and the considerably softer spectrum seen by the High Altitude Water Cherenkov Observatory (HAWC) at near 103 GeV. We show that lower-energy (≲10 GeV) gamma rays are primarily produced in the network elements and higher-energy (≳few × 10 GeV) gamma rays in the intergranular sheets. Notably, the spectrum softening observed by HAWC results from the limited effectiveness of capturing and reflecting ∼104 GeV cosmic rays by the finite-sized intergranular sheets. Our study is important for understanding cosmic-ray transport in the solar atmosphere and will lead to insights into small-scale magnetic fields at the photosphere.


INTRODUCTION
The Sun is a bright and time-steady source of GeV to TeV gamma rays, with the intensity showing a modest anti-correlation with the solar cycle.This emission is proposed to arise from two distinct channels, originating from different locations.The first is emission from the solar halo, arising from electron galactic cosmic-ray (GCR) interactions with solar photons via the inverse Compton process near the Sun.Theoretical predictions (Moskalenko et al. 2006;Orlando & Strong 2007, 2021) reasonably match observations (Abdo et al. 2011).The second source of emission arises from the solar disk, where the decay of π 0 and other mesons produced through nucleon-nucleon collisions occurs as a result of the bombardment of hadronic GCRs on the solar Corresponding author: Jung-Tsung Li li.12638@osu.edusurface.Theoretical predictions for the solar disk emission do not match the observed data.Consequently, in this paper, we focus on understanding this emission.
Observations with the Fermi Large Area Telescope (Fermi-LAT) demonstrate that the solar disk emits gamma rays in the 0.1-200 GeV range with a spectrum, dN γ /dE γ , following an ∼ E −2.2 γ power law (Abdo et al. 2011;Ng et al. 2016;Tang et al. 2018;Linden et al. 2022).(For the earlier work with EGRET data, see Orlando & Strong (2008).)More recently, the High Altitude Water Cherenkov Observatory (HAWC) reported the first detection of gamma-ray fluxes near 10 3 GeV, finding a E −3.6 γ power law (Albert et al. 2023).Both instruments reveal an anti-correlation between gammaray fluxes and solar activity, with fluxes during the solar minimum being approximately a factor ∼ 2 greater than those during the solar maximum.Interestingly, while the entire solar disk emits gamma rays, indicating that the emission is not just a limb effect, the emission distri-Li et al.
bution reported in Linden et al. (2018) shows moderate non-uniformity across the solar surface.
The pioneering theoretical work on this subject is Seckel et al. (1991) and Seckel et al. (1992), which focus on the reflection of hadronic GCRs by canopy fields via the magnetic mirroring effect.Their work assumes a simple pressure balance following a P (z) ∝ B (z) 2 scaling relation, where P (z) and B (z) are the gas pressure and the magnetic field strength of the flux tube at the height z, respectively.They further assume that particle trajectory is governed by the adiabatic invariance, cos θ (z) = 1 − B (z) sin 2 θ 0 /B 0 , where θ (z) is particle pitch angle at z, θ 0 is the initial pitch angle at the top of the tube, and B 0 is the magnetic field strength at the top of the tube.A recent study conducted by Hudson et al. (2020) suggests a scaling relation of P ∝ B 3.5±0.1 may be more appropriate based on data from the Bifrost magnetohydrodynamics (MHD) simulation code (Gudiksen et al. 2011).Several other approaches and ideas have been suggested in the literature.The studies by Mazziotta et al. (2020) and by Li et al. (2020) have integrated the Potential Field Source Surface model of the solar coronal fields into the FLUKA code and implemented numerical simulations of the yields of secondary particles and gamma rays from hadronic GCR showers.Gutiérrez & Masip (2020) and Gutiérrez et al. (2022) use HAWC data on the cosmic ray shadow of the Sun to deduce the averaged hadronic GCR absorption fraction and the gamma-ray yields.Banik et al. (2023) presents a contrasting perspective, arguing that the acoustic-like shock waves in the chromosphere can accelerate protons in the solar gas up to 10 3 GeV and produce gamma rays that match the HAWC observation.Zhou et al. (2017) argues that the solar magnetic fields are negligible for gamma-ray production at high enough energies, in which case the disk emission would be solely from the thin ring of the solar limb; the HAWC data show that this must be at gamma-ray energies above about 10 3 GeV.Overall, despite extensive theoretical efforts in this domain, no prediction can simultaneously account for the overall gamma-ray flux, the morphology of its emission, and its time variability.
In this paper, we aim to identify the magnetic field structures at the solar surface that are crucial for reflecting hadronic GCRs and to model the gamma-ray emission from these sites.We construct a double-structure solar surface magnetic field model comprising one vertical flux tube, representing the magnetic flux tubes that form the network elements, and one vertical flux sheet, representing the magnetic flux sheets at the downflow lanes between granules.Both the tube and sheet, which have horizontal sizes of ∼ 100 km, maintain magnetohy-drostatic equilibrium with the surrounding gas.Their field geometries are calculated following the numerical approach in Steiner et al. (1986).We inject proton GCRs at the top of the flux tube and sheet, numerically trace their trajectories, and calculate their gamma-ray emission.By doing so, we can determine the fraction of GCRs reflected in the flux tubes of the network elements and the fraction that plunges through to the flux sheets of the intergranule lanes.Our model helps us understand the relative roles of these components in capturing and reflecting GCRs across various energy ranges.Comparison to observations helps probe the solar magnetic environment, even below the surface of the photosphere (z = 0 km, defined by where the optical depth τ 5000 = 1 for light of wavelength 5000 Å).In our calculations of the gamma-ray emission, we aim for a precision of a factor of two, which is appropriate given the large dynamic ranges of the variables and the substantial uncertainties.
Last, it is important to note that while solar flares are known to generate gamma rays through reconnection and non-thermal particle acceleration, the emission is generally episodic and does not exceed a few GeV in maximum gamma-ray energy (Murphy et al. 1987;Schneid et al. 1996;Ackermann et al. 2014;Ajello et al. 2014;Pesce-Rollins et al. 2015;Omodei et al. 2018;Share et al. 2018).Moreover, such events show transient and burst-like behavior and are easily excluded from longduration gamma-ray observational datasets.As a result, solar flares are not relevant to the time-steady, multi-GeV gamma-ray emission discussed in this paper.
The remainder of this paper is organized as follows.We begin in Section 2 with an introduction to the solar magnetic field structure.Section 3 introduces the double-structure model of gamma-ray emission for the quiet photosphere.In Section 3, we also discuss the potential effect of GCR scattering by magnetic fluctuations.Section 4 presents the magnetohydrostatic solutions for the flux tube and flux sheet.Section 5 explains our simulation setup for hadronic GCR transport in the double-structure model.Section 6 presents the simulation results.We conclude in Section 7. In Appendix A, we show the most probable injection polar angles of GCRs for producing the observed gamma rays.In Appendix B, we show the average gamma-ray emission angle and depth.In Appendix C, we validate the assumption of collinearity in gamma-ray production.

SIMPLIFIED FRAMEWORK FOR SOLAR SURFACE MAGNETIC FIELDS
In this section, we provide a brief overview of photospheric convection, solar-surface magnetic fields, and open field lines.The reason for this focus will later be explained in Section 3.
Figure 1 is a schematic diagram of the quietphotosphere magnetic network structures that form the open field lines.Here we follow the canonical picture outlined in Cranmer & van Ballegooijen (2005) and Wedemeyer-Böhm et al. (2009).We begin our discussion with convection in the photosphere, as it is the driving force behind the formation of magnetic flux tubes.Convection near the photospheric surface manifests across various scales.The smallest observable scale is granulation, with a typical horizontal size of about 1000-1500 km (white polygons in Figure 1(a)).Each granule cell hosts a convective flow where hot material rises from the cell center, moves horizontally at the top, and descends at the cell edge.Between the granule cells are the intergranule lanes, where the flow is compressed and pushed downward (gray bands in Figure 1(a)).At a larger scale is supergranulation, with the typical horizontal size of about 10,000-30,000 km (largest arcs in Figure 1(b)).The convective flow responsible for supergranulation occurs on larger scales and at greater depths than the convective flow driving granulation.
Following the variety of scales in convection, the magnetic fields also show a variety of scales.In the case of network fields, the smallest-scale feature is vertical magnetic flux tubes at the photosphere.These tubes, with diameters of approximately 150 km and field strengths of around 1500 Gauss, are situated in the intergranule lanes at the edges of supergranule cells (Figure 1(a)).Each flux tube maintains magnetohydrostatic equilibrium with the surrounding gas.As the gas pressure decreases with increasing height, the flux tubes expand horizontally.A collection of nearby flux tubes eventually merges at a height of around 600-1000 km, forming a network element.At the merging height, the horizontal scale of each network element is about 2000-6000 km.Above the merging height, the magnetic fields transition into a thicker flux tube and continue to expand with height (upper-half of Figure 1(b)).A second merging occurs when the fields of two adjacent network elements, separated by a horizontal distance of about a supergranule cell, coalesce at the merging height of a few 1000 km above the surface.The magnetic fields continue to extend into the coronal region, becoming parts of the large-scale coronal-hole network fields.These network fields eventually extend into interplanetary space as open magnetic field lines (Figure 1(c)), which have their other footprint in the outer structure of the heliosphere.
Magnetic fields between network elements in the quiet photosphere are internetwork fields.Observations in-dicate that these regions comprise small-scale, mixedpolarity magnetic fields that show loop-like structures a few hundred km above the solar surface (dashed arcs in Figure 1(b)) (Orozco Suárez et al. 2007;Lites et al. 2008).These fields are situated above granule cells, with their polarities connecting to the intergranule lanes.Below the photospheric surface, strong vertical magnetic fields are formed in the intergranule lanes as a result of flux expulsion (Parker 1963;Weiss 1966;Galloway & Weiss 1981) and field amplification (Parker 1978;Webb & Roberts 1978;Spruit & Zweibel 1979) by the granule convective flows.In the near-surface environment, they have sheet-like structures in the intergranule lanes between granules and tube-like structures in the micropores located at the vertices where multiple intergranule lanes converge.Because direct observations below the surface are unavailable, the extent of these structures into the upper convection zone is uncertain.However, valuable insights have been obtained from state-ofthe-art numerical simulations of magnetoconvection in the solar photosphere (Weiss et al. 1996;Stein & Nordlund 1998;Steiner et al. 1998;Cattaneo 1999;Cattaneo et al. 2003;Schaffenberger et al. 2005;Vögler et al. 2005;Gudiksen et al. 2011;Freytag et al. 2012).These simulations demonstrate that the coherent structure of flux sheets at the intergranule lanes extends a few hundred km below the surface before being disrupted by turbulent flows.
The discussion thus far restricts to the quiet region of the Sun where magnetic structures near the solar surface are primarily small-scale and driven by granule and supergranule convection cells.Other regions on the solar surface include active regions with more complex magnetic structures.Active regions host strong, large-scale magnetic features, such as sunspots, coronal loops, solar flares, and coronal mass ejections.The activities within the active regions and the solar cycle are thought to stem from the magnetic fields generated by the large-scale solar dynamo emerging from the tachocline at the base of the convection zone.

SETUP OF THE GAMMA-RAY EMISSION MODEL IN THE QUIET PHOTOSPHERE
In this section, we provide the setup for our model for the time-steady solar gamma-ray emission of the quiet photosphere.Given the extensive range of scales discussed in Section 2, simulating GCR propagation and gamma-ray emission in the solar surface environment poses significant challenges.In the following, we explain the approach of our model, which optimizes the prediction of the gamma-ray emission based on appropriate approximations.We review solar modulation and hadronic interactions, then estimate the relevant depths in the photosphere for GCR interactions, outline our key assumptions, define the key magnetic field structures that we model, and discuss the potential effect of the magnetic fluctuations on particle trajectories.

Solar Modulation
As the GCR flux has not yet been measured near the Sun, it must be calculated based on the flux near Earth.When hadronic GCRs in interplanetary space propagate toward the Sun, their fluxes decrease due to the interactions with magnetic turbulence in the solar wind, a phenomenon known as solar modulation (Parker 1965;Gleeson & Axford 1968).In our previous work in Li et al. (2022), we utilized recent magnetic power spectral density measurements from the Parker Solar Probe to calculate the parallel diffusion coefficients.We then showed that the proton GCR flux is reduced by ≲ 15% for proton kinetic energy, E k p , in the range of 0.1 GeV < E k p < 1 GeV when transported from a heliocentric distance of 1 au to 0.1 au.For the higher energy range of GCRs considered here, the effects are even smaller, O (≲ 1%) for E k p ≳ 10 GeV.The GCR transport from 0.1 au to the solar surface (≈ 0.005 au) remains uncertain.Accurate evaluation of this transport requires theoretical models of magnetic turbulence, solar wind acceleration, and coronal loops.This region was not considered in our previous work in Li et al. (2022) and is beyond the scope of the present work.(See e.g., Petrosian et al. (2023) for modeling of electron GCR propagation down to the solar surface.)However, we expect that the modulation from 0.1 au to the solar surface should be no more than a few percent for E k p ≳ 10 GeV based on the model extrapolation in Li et al. (2022).Given that the focus of gamma rays in this work is above 1 GeV, hence proton GCRs above 10 GeV (see below), the modulation effects on our predictions should be no more than ≃ 10%, and are hence neglected.
In this work, we use the GCR spectra measured at 1 au from the Sun, which takes into account the significant modulation effects experienced as GCRs propagate from the heliopause to 1 au, for GCR injection into the solar surface magnetic fields.We have further assumed an isotropic GCR distribution.We use the proton GCR measurements from AMS-02 (Aguilar et al. 2015) for 0.5 GeV ≤ E k p ≤ 1.46 × 10 3 GeV and from CREAM (Choi et al. 2022) for 2.10 × 10 3 GeV ≤ E k p ≤ 5.33 × 10 5 GeV, linearly interpolating the fluxes between the data points.In Section 5.5, we describe how we take into account the moderate effect of helium in the GCRs and the Sun.

Hadronic GCR Interactions
The interactions of hadronic GCRs in the solar atmosphere can be understood based on well-established results on how these interact in Earth's atmosphere (Gaisser et al. 2016).Because the GCR energies considered are so large, the details of the atomic states (e.g., ionization) are wholly irrelevant, and even the details of the nuclear states (e.g., hydrogen vs. nitrogen) are largely so.
The dominant loss process for proton GCRs is inelastic nucleon-nucleon collisions ("pp interactions" below) that produce pions and other mesons.In our calcula-tions, we use the results of Kelner et al. (2006), plus validations with the FLUKA code (Battistoni et al. 2015;Ahdida et al. 2022).Though multiple mesons may be produced, the basic physics is well summarized by p + p → p + p + π 0 , followed by π 0 → γ + γ.The total inelastic hadronic cross section for proton-proton interactions is σ pp ≃ 3 × 10 −26 cm 2 , with weak energy dependence for E k p ≳ 1 GeV (Workman et al. 2022).The proton mean free path, expressed as a mass column density, is thus ≃ 55 g cm −2 (about half a meter of water equivalent), meaning that gamma rays are only produced when the GCRs encounter substantial material.The differential cross section can be roughly approximated by assuming that the most important gamma ray has E γ ∼ 0.1 E k p .Last, these pp interactions may break nuclei, either directly or through induced showers, and this may lead to nuclear de-excitation gamma rays; as those are below about 10 MeV, we ignore them here.
The second-most important -but ignorable -loss process is continuous losses due to GCR collisions with electrons.For neutral material, the loss rate due to atomic excitation and ionization is approximately 2 MeV/ g cm −2 for all materials except hydrogen and 4 MeV/ g cm −2 for hydrogen, the difference is due to the greater number of electrons per unit mass (Workman et al. 2022).The solar photosphere is largely neutral, but we note that even the losses in fully ionized material (so-called Coulomb losses) are only a few times larger (Strong & Moskalenko 1998).These loss rates have only weak dependence on energy in our range of interest.The continuous energy loss of a proton over one mean free path for inelastic nucleon-nucleon interactions is thus ≃ 0.2 GeV.Due to our focus on gamma rays above 1 GeV, and hence protons above about 10 GeV, this continuous energy loss is negligible.

Relevant Photospheric Depths
Using the proton mean free path, we estimate the possible depth range for interactions, which informs which magnetic field structures are most relevant.For simplicity, here we only assume vertical trajectories, ignoring helical motion.Interactions typically occur when σ pp n H (z) dz ≳ 1, where n H (z) represents the hydrogen number density of the gas at height z.
Using the mass density ρ of the Sun, which will later be given in Figure 3, we estimate that incoming proton GCRs typically penetrate no deeper than 800 km below the surface.Because the mass density varies exponentially with depth, this is a firm bound.As a result, proton GCR absorption and gamma-ray emission primarily occur within the photosphere and uppermost convection zone.
In Appendix B, we show that a full calculation with our model reveals that emission primarily occurs within a height ranging from −100 km to 500 km.This is shallower than noted in our estimate because our 3D numerical simulation of particle trajectory takes into account full helical motion, hence a longer particle trajectory at around the reflection height compared to the case for the vertical trajectories.Moreover, to produce outgoing gamma rays, proton GCRs must first be reflected from incoming to outgoing.

Key Assumptions of our Model
The Sun's magnetic environment is complex, with distinctive features ranging from the smallest scales to the size of the solar system.However, when considering how incoming GCRs produce outgoing gamma rays, there are basic facts that allow significant simplifications.First, hadronic GCRs only produce gamma rays when the accumulated mass column density they traverse is very large.Because the matter density in the solar surface varies exponentially, only a narrow range of heights is relevant, as discussed in Section 3.3.This suggests that, as a primary effect, we should consider magnetic fields from network elements and intergranule lanes within this specific height range.Second, there must be magnetic fields strong enough to reflect the very energetic GCRs from incoming to outgoing so that the gamma rays can escape over the whole surface of the Sun, as observed.In comparison, we can reasonably neglect deflections from other field structures, though they may matter at second order through how they change the angular distribution of incoming GCRs.
In our model, we make the following simplifications: 1.That the solar surface is spherically symmetric, and thus the model presented in this work applies to every small patch on the surface of the Sun.
2. That all GCRs reach the flux tube of the network elements in the quiet photosphere by propagating along the open field lines and network fields.
3. That higher-energy GCRs that penetrate the flux tubes of the network element plunge onto internetwork regions, which consist of vertical flux sheets.
The first simplification allows us to focus on the main parts of the emission process; however, it means we will not be considering the different latitudes of the coronal holes and open field lines.The second simplification aids in understanding how emission occurs in the quiet photosphere, even though we will not explore the relationship between gamma-ray flux, active regions, and solar activity in this paper.The third simplification helps us separate the trajectories of low-and high-energy GCRs and gives us a method to evaluate the fraction of highenergy GCRs entering the internetwork regions.

Details of the Double-Structure Model
Figure 2 shows the schematics of the double-structure magnetic fields we use to model GCR propagations.It consists of one isolated flux tube and one isolated flux sheet.The tube represents those forming the network elements, as depicted in Figure 1(a).The sheet represents those formed in the intergranule lanes, as depicted in Figure 1(b).In our model, this sheet is symmetric along the long horizontal axis and extends indefinitely.The magnetic structures of the flux tube and flux sheet are essential to determining the 3D trajectories of the GCRs.In Section 4, we demonstrate the numerical solutions for these structures under the magnetostatic equilibrium with the surrounding gas.
Driven by our assumption that hadronic GCRs propagate towards the Sun along open field lines and network fields, eventually reaching the network elements, we start the GCR trajectories at the top of the flux tube.Hadronic GCRs with an isotropic distribution are uniformly injected across the cross-section surface of the tube at the merging height, set at z = 1600 km, as shown in Section 4. Lower-energy GCRs -those with Larmor radii, r L , much smaller than field scale height, (∇ ln|B|) −1 -are tightly confined by the field lines within the flux tube.As they spiral downward, the magnetic field strength increases due to compression from the surrounding gas, causing the particle pitch angle to approach 90 • .Eventually, the radial component of the field imparts a "kick" to the particles, initiating the upward spiral, the process known as the magnetic bottle effect.
On the other hand, higher-energy GCRs -those with r L ≫ (∇ ln|B|) −1 -do not spiral inside the tube; instead, they plunge through the tube with a slight deflection.After they exit the tube, we maintain their polar angles (relative to the vertical axis) while isotropizing their azimuthal angles.These GCRs are then injected into the flux sheet, onto a horizontal surface at z = 600 km, as depicted in Figure 1(b) and Figure 2. Let L be the separation between the two adjacent granule lanes, and W sh be the width of the flux sheet at z = 600 km.Then the ratio W sh /L is the fraction of those higher-energy GCRs exiting the tube that should be injected into the sheet.For the GCRs that enter the sheet, we calculate their particle trajectories and gamma-ray yields.
At the very highest energies, GCRs are not captured by solar-surface magnetic fields.They will plunge straight through the flux sheets and interact with the Sun but do not produce outgoing gamma rays.

Magnetic Fluctuations in the Solar Atmosphere
The magnetic field structure described in Section 3.5 is a magnetostatic configuration, which does not contain any magnetic fluctuations.In this subsection, we discuss the potential effect of magnetic fluctuations on particle trajectories.
We distinguish the impacts of magnetic fluctuations originating from three distinct regions within the solar atmosphere in the open field scenario.The first region is the solar corona, where Alfvénic turbulence occurs.This turbulence is thought to contribute to the heating of the plasma in the coronal region and the acceleration of the solar wind (Tu 1987(Tu , 1988;;Matthaeus et al. 1999;Cranmer et al. 2007;Verdini et al. 2010;van Ballegooijen et al. 2011;van Ballegooijen & Asgari-Targhi 2016).The energy source of the turbulence comes from the Alfvén waves launched from the photospheric surface due to the buffeting of the footprints by solar granulation flows (van Ballegooijen et al. 1998;Nisenson et al. 2003;Cranmer & van Ballegooijen 2005).As these waves propagate outward from the photospheric surface, they interact with the inward-propagating waves arising from wave reflection in the lower coronal region (Heinemann & Olbert 1980;Velli et al. 1989;Velli 1993;Matthaeus et al. 1999;Dmitruk et al. 2002;Dmitruk & Matthaeus 2003;van Ballegooijen et al. 2011;van Ballegooijen & Asgari-Targhi 2016;Meyrand et al. 2023).The energy cascade then happens due to the nonlinear interactions between the counter-propagating Alfvén waves, producing microscopic Alfvénic turbulence at smaller scales (i.e., higher wavenumbers) (Iroshnikov 1964;Kraichnan 1967).The turbulence generated this way in the solar corona is transported away from the Sun by the solar wind into the heliosphere.We point out that the origin of the solar modulation discussed in Section 3.1 is due to GCR interaction with this kind of turbulence from the solar wind.The reduction of GCR intensity due to solar modulation within 1 au is small (Li et al. 2022) and thus neglected in this work.Because the gamma-ray flux produced in the corona is much smaller than that produced near the photospheric surface, as discussed in Section 3.3, the magnetic fluctuations in this region are not the primary concern of this work.
The second region is the photosphere and the lower chromosphere, where the gamma rays are most likely to be produced, as shown in Section 3.3.Specifically, we refer to the region from z = 0 km to z = 1600 km, above which the density is not high enough to produce significant gamma rays.The magnetic fluctuations in this region are macroscopic MHD transverse waves (fluctuations) driven by photospheric motions and the waves reflected from the lower coronal region.(Specifically, the waves below the merging heights of the flux tubes at the height of 600-1000 km are classified as kink waves, while those above the merging height are identified as Alfvén waves.Both types are referenced as MHD transverse waves in our discussion.)Cranmer & van Ballegooijen (2005) highlights that the energy cascade efficiency is low in this region in the open field scenario, which is also supported by the findings in Cranmer et al. (2007).This low efficiency is evident in Figure 14 of Cranmer & van Ballegooijen (2005), which illustrates that the Alfvén wave reflection time, t ref , is shorter than the nonlinear outer-scale eddy cascade time, t eddy .Following the hierarchy of the energy cascade by Dmitruk & Matthaeus (2003), a weak energy cascade occurs when t ref < t eddy , as Alfvén waves propagate away before the turbulence has sufficient time to develop and heat the plasma.In contrast, the strong energy cascade occurs when t eddy < t ref , as there is sufficient time for the Alfvén turbulence to develop.Therefore, the condition in Cranmer & van Ballegooijen (2005) of t ref < t eddy in the open field regions of the photosphere and chromosphere suggests a weak magnetic power spectrum of the Alfvén turbulence for scales smaller than the outer scale of the driving eddies, i.e., an inefficient process to produce smaller-scale sideway displacements on the magnetic field lines.
This weak turbulence suggests that particle scattering with Alfvén turbulence is likely also weak, though whether it is negligible in the study of solar gammaray emission requires a careful evaluation against pitch-angle scattering and magnetic mirroring times.Similar analysis has been conducted in Malyshkin & Kulsrud (2001) and Effenberger & Petrosian (2018) on the topics of thermal conduction within stochastic magnetic fields and their application to energetic particle transport in solar flares.The two references indicate that if the escape time in the magnetic mirroring region is less than the pitch-angle scattering time, then magnetic mirroring due to the converging fields has a stronger effect than the scattering of the particles on the Alfvén turbulence, and the opposite is true if the scattering time is less than the escape time.However, a precise quantitative analysis of particle pitch-angle scattering requires a detailed magnetic power spectrum of Alfvén turbulence as functions of wavenumber and height in the photosphere and chromosphere.Such information has not been shown in the existing literature.The derivation of this magnetic power spectrum would entail modeling macroscopic MHD and acoustic waves originating from photospheric motions, their reflection at the coronal base and lower corona, nonlinear wave interactions for direct and inverse energy cascade processes, and a plasma heating model that aligns with observational data.These requirements exceed the scope of the present study.As such, our current flux tube and sheet model does not consider microscopic Alfvén turbulence, though this aspect remains a focus for future investigations.
The third region is the upper-most convection zone, situated within a few hundred km below the photospheric surface.The macroscopic fluctuations of the magnetic fields and the fluids are caused by the convective nature of the granule cells.Even in this convective environment, a coherent structure in tubes and sheets can emerge due to flux expulsion (Parker 1963;Weiss 1966;Galloway & Weiss 1981) and field amplification (Parker 1978;Webb & Roberts 1978;Spruit & Zweibel 1979), and it is this coherent structure that we investigate in this study.In contrast, the macroscopic fluctuating parts of the fields cannot be easily obtained through theoretical approaches.They require 2D or 3D magnetoconvection simulations, which exceed the current paper's scope but will be pursued rigorously in future research.In this study, we omit the macroscopic fluctuating component and the microscopic Alfvén turbulence of the magnetic fields, focusing solely on the coherent component of the field configuration in the uppermost convection zone, photosphere, and chromosphere.It is noteworthy that our results, as presented herein, reasonably agree with the observational gamma-ray data, signifying that the double-structure approach in a magnetostatic configuration setting captures key character-

Li et al.
istics of the hadronic GCR reflection and the resulting solar gamma-ray emission.

MAGNETOHYDRODYNAMIC SOLUTIONS OF FLUX TUBE AND FLUX SHEET
In this section, we describe the magnetic flux tube and flux sheet structures employed in our model.We first present the equations governing the magnetohydrostatic equilibrium with the surrounding gas, then develop the numerical solutions, for which we follow the approach in Steiner et al. (1986).We emphasize that the coordinate variable names used in this section should not be conflated with those used in Section 5 for particle trajectories.

Basic Equations
Flux tube: The flux tube representing the network elements is approximated as a vertical, axisymmetric, untwisted magnetic flux tube that is in magnetohydrostatic equilibrium with the surrounding gas.We describe the flux tube in the cylindrical coordinates (r, ϕ, and z).The magnetic field B of this geometry has azimuthal symmetry, i.e., ∂B/∂ϕ = 0.By introducing the stream function Ψ ≡ rA ϕ where A ϕ is the ϕ component of the vector potential, the field components can be described by The stream function Ψ satisfies the quasi-linear (Grad-Shafranov) equation, where J is the ϕ-directional current density (with the dimension of [current/length 2 ]) at the surface between the flux tube and the surrounding gas.
Flux sheet: The flux sheet representing the intergranule sheets is approximated as a vertical, untwisted magnetic flux sheet in a magnetohydrostatic equilibrium with the surrounding gas.We describe the flux sheet in Cartesian coordinates (x, y, and z).The magnetic field of this geometry has translational symmetry, i.e., ∂B/∂x = 0.The stream function is the y component of the vector potential, i.e., Ψ = A y .The field components are described by where Ψ satisfies the quasi-linear equation, with J the x-directional current density at the surface between the flux sheet and the surrounding gas.

Common features:
The magnetic fields at the surfaces of the flux tube and flux sheet are discontinuous and thus form current sheets with surface current density J ⋆ (with the dimension of [current/length]).In this paper, we refer to this current sheet as the "flux boundary".For magnetohydrostatic equilibrium, the total pressure, P +B 2 /8π, must be continuous across the flux boundary, i.e., where the subscripts i and e refer to the internal and external boundary surfaces of the flux tube or flux sheet.Using Ampere's Law and the pressure continuity in Equation ( 5), J ⋆ can be expressed as To solve the quasi-linear equations in Equations ( 2) and ( 4), we make several simplifications about the surrounding gas.First, we consider an idealized case where the gas temperature T at any given height z is in equilibrium in the horizontal direction, i.e., T i (z) = T e (z) = T (z).Thus, the interior and exterior gas pressure scale heights are equal and are given by where R ⊙ is the solar radius, k b is the Boltzmann constant, G is Newton's constant, µ is the mean molecular weight, and M ⊙ is the solar mass.The interior and exterior gas pressures are then expressed as where P i (z b ) and P e (z b ) are the internal and external gas pressures immediately adjacent to the flux boundary, and z b is the base of the computation box.Additionally, we assume that at z = z b , the vertical magnetic field components inside the flux tube or flux sheet have a uniform distribution, while outside the flux tube or flux sheet, the magnetic fields are zero.As a result, the difference between the interior and exterior gas pressures at the base is expressed as

Numerical Solutions
We follow the iterative procedure from Steiner et al. (1986) to solve for Ψ from the quasilinear equations in Equations ( 2) and (4).Because of the symmetries of the flux tube and flux sheet, we only solve half of the tube domain (the rz plane) in Equation ( 2) and half of the sheet domain (the yz plane with y ≥ 0) in Equation (4).In both cases, the boundaries at the top of the computational domains follow the Neumann condition where ∂Ψ/∂z = 0 assuming B r = 0. To obtain stable solutions, we have used the implicit underrelaxation method introduced in Patankar (1980).
To evaluate the pressure scale height in Equation ( 7), we must also specify T and µ as a function of z.For this purpose, we use the HSRASP model from Chapman (1979).This model combines the HSRA model from Gingerich et al. (1971), describing the gas properties in the chromosphere and photosphere, with the upper convection zone model from Spruit (1974).
Figure 3 shows the HSRASP model for the upper convection zone, photosphere, and chromosphere.The photosphere layer is located in 0 km < z < few × 100 km, with the chromosphere located above and the convection zone below.Because the HSRASP model does not provide ρ, P , and T above z = 1000 km and µ above z = 0 km, we make assumptions to extrapolate these data.First, we set T = T (1000 km) = 6070 K for 1000 km ≤ z ≤ 1600 km.This temperature aligns with the model Last, we assume that Sun's atmosphere is comprised of only hydrogen (H) and helium (He) atoms.Heavier elements are ignored since they only account for ≲ 2% of the total mass.As a result, the mass fraction of the hydrogen atom is given as where A H = 1 and A He = 4 are the mass numbers of hydrogen and helium, respectively.

Flux Tube Solution
For the numerical solution of the flux tube, the computational domain is a 2800 km × 500 km box that is discretized on a rectangular mesh of 65×65 mesh points.The base of the box is set at z b = −1200 km.We select an initial radius for the flux tube, R ⋆ , and vertical magnetic field strength, B z (z b ) at the base of the computa-

Li et al.
tional domain, with values of 36 km and 10, 000 Gauss, respectively.This choice of initial conditions gives the radius of the tube at z = 0 km as 80 km, the filling factor as 3.6%, and the axial field strength at r = z = 0 km as 1580 Gauss.The numerical value of the tube radius at z = 0 km is close to the photometric measurements reported in Muller & Keil (1983), which shows that the typical flux tube diameter at z = 0 km for the facular points forming the network fields in the quiet photosphere is 150 km.The numerical values of the filling factor and the axial field strength agree with the network properties revealed by the Stokes I and V lines of Fe-I reported in Solanki & Stenflo (1984), which show that the vertical field component of the network field ranges from 1400 G to 1700 G, and the filling factor is between 3% to 4%.
Figure 4(a) shows the numerical solution for the flux tube, which is invariant in the azimuthal direction.The radius of the tube at the top of the simulation box (z = 1600 km) is R top = 430 km.The figure is divided into groups of subfigures to show the variations for B z and B r at different heights.At each height, we find that as r increases, B r increases while B z decreases.
Lastly, we emphasize that the magnetic field structures discussed here are located in the photosphere and chromosphere.They are finite-sized and the building blocks of the large-scale coronal and interplanetary magnetic fields.These flux structures located at the photosphere have much stronger field strength (> 1000 G) than those in the coronal region, which have a field strength of ∼ few G at 1 R ⊙ from the photospheric surface (Patzold et al. 1987).The large field strength of the structures we consider, plus their location in dense regions, is essential to redirecting GCRs from inward to outward below the column density of material needed to produce gamma rays.

Flux Sheet Solution
For the numerical solution of the flux sheet, the computational domain is an 1800 km × 500 km box that is discretized on a rectangular mesh of 65 × 65 mesh points.The base of the computational box is set at z b = −1200 km, as for the tube.We choose the initial half-width of the sheet and the vertical magnetic field strength B z (z b ) at the base of the computation box to be 30 km and 10, 000 Gauss, respectively.This choice of the initial conditions gives the width of the sheet at z = 0 km as 330 km, the filling factor as 55%, and the vertical field strength at r = z = 0 km as 1608 Gauss.The half-width of the sheet at the top of the simulation box, z = 600 km, is given as y top = 330 km.Consequently, W sh = 2y top = 660 km.The numerical values of the widths at z = 0 km and at 600 km, as well as the vertical field strength at z = 0 km, are similar to the results obtained from the 3D magnetoconvection simulation presented in Freytag et al. (2012).
Figure 4(b) is the numerical solution for the flux sheet, which is invariant along the x direction.(The sheet structure is redundant in the x direction.)The interpretation is the same as Figure 4(a), except for the blue dashed line denoting B y variation in this model.

SIMULATION APPROACH FOR GCR PROPAGATION AND GAMMA-RAY EMISSION
In this section, we explain our numerical simulation of solar disk gamma-ray emission.Our simulation approach is factorized into four separate stages, taking advantage of the fact that we need to calculate only the effects on the average GCR proton and gamma-ray distributions, as opposed to needing to follow the details of every particle.First, we simulate proton GCR trajectories in the flux tube and flux sheet, recording the accumulated optical depths of these trajectories, but ignoring interactions.Next, we calculate the energy spectrum of gamma rays resulting from pp interactions, based on the proton optical depths.We then integrate the gammaray fluxes along the trajectories of the proton GCRs to obtain the total emission flux from the solar disk due to pp interactions.Finally, we incorporate helium contributions to the GCR flux and the solar gas density by using the nuclear enhancement factor.

Proton GCR Propagation
Proton GCR motions follow the Lorentz force equation, where v is proton velocity, Γ is the Lorentz factor, q is proton charge, m p is proton rest mass, c is the speed of light, and B (r) is the local magnetic field at location r = (x, y, z).We do not include the electric field E in the Lorentz force equation as its magnitude is negligible compared to the v × B/c term.We motivate this within an ideal magnetohydrodynamics framework where E = −U × B/c, with U denoting the plasma flow velocity.
The typical convective flow speed in the granule cell is |U | ∼ 1 km/s, five orders of magnitude smaller than the hadronic GCR speed, |v| ≈ c.

Injection into the Flux Tube
We uniformly inject protons across the horizontal cross-section surface of the flux tube at z = 1600 km.Due to the azimuthal symmetry of the tube, we only inject protons along one of the horizontal axes.We begin the injection process at an axial distance of r 0 = 20 km and progressively increase r 0 in increments of 40 km until it reaches the boundary of the flux tube.
At each r 0 , we inject proton GCRs based on the following procedure, with a focus on those propagating into (downward) the flux tube while disregarding those moving away (upward) from it.We define θ 0 and ϕ 0 as the initial polar and azimuthal angles, respectively, relative to the vertical direction.We consider a range for θ 0 spanning from 90 • to 180 • (all downward directions) with increments of ∆θ 0 = 1 • , and a range for ϕ 0 spanning from 0 • to 360 • with increments of ∆ϕ 0 = 45 • .For every combination of values (r 0 , θ 0 , ϕ 0 ), we evaluate E k p ranging from 1 GeV to 10 5 GeV, dividing this range into eight equally spaced logarithmic E k p bins per decade.We numerically simulate 3D proton trajectories in the flux tube.
For each combination of values r 0 , θ 0 , ϕ 0 , E k p , one proton GCR is injected.The simulation of each proton trajectory is terminated when the proton satisfies one of the following conditions: (1) z < −1200 km, (2) z > 1600 km, or (3) exiting the surface of the tube.In (1), protons are not magnetically reflected before being absorbed by the Sun.We select a minimum height of the simulation box at −1200 km where the pp interaction probability exceeds 95%.This is to ensure a sufficiently high probability of pp interactions at the bottom of the simulation box.In (2), protons have been reflected and exit the tube from the top cross-section surface, above which the gas density is too low to produce significant gamma rays.In (3), protons pass through the edge of the tube; those with downward velocities later enter the internetwork regions.

Injection into the Flux Sheet
We apply a similar injection process at the top crosssection of the flux sheet at z = 600 km (see Figure 1(b) and Figure 2).We only inject proton GCRs along the y direction, as B is invariant in the x direction.Starting at y 0 = 20 km, we increment y 0 by 40 km until y 0 reaches the boundary of the flux sheet.
We maintain the same parameter values for θ 0 , ϕ 0 , and E k p as in the flux tube.We use θ0 and φ0 to denote the initial polar and azimuthal angles at the injection site of the flux sheet.For each combination of values y 0 , θ0 , φ0 , E k p , one proton is injected.The simulation of each proton trajectory is terminated when the proton satisfies one of the following conditions: (1) z < −1200 km, (2) z > 600 km, or (3) exiting the surface of the sheet.
We note that the simulation procedure for the flux sheet is independent of the injected GCR spectrum.However, the injected GCR spectrum for the flux sheet is anticipated to be different from the injected GCR spectrum for the flux tube.This is because low-energy proton GCRs are confined within the flux tube, whereas high-energy GCRs traverse the tube with a slight angular deflection, as described in the double-structure model explanation found in Section 3.5.Therefore, only the high-energy GCRs can exit the flux tube and subsequently enter the flux sheet.To address this effect, we separately compute the distribution for GCR injection into the flux sheet in Section 5.1.3.

Angular and Energy Efficiency of GCRs Injected into the Flux Sheet
In the simulation of proton GCR injection into the flux tube presented in Section 5.1.1,we record the ones that pass through the flux tube surface at heights above z = 0 km.For each r 0 and E k p , we record the number of proton GCRs, Q, and their polar angles θ ′ (relative to the ẑ direction) upon exiting the tube surface.We disregard the azimuthal angular dependence in Q, assuming each proton enters the flux sheet isotropically in the azimuthal direction at the injection site.To account for the higher number of protons injected at larger r 0 , we weight Q over the ring area with the radius r 0 and the width dr 0 , i.e., where A tot = πR top 2 .Note that if there were no magnetic field within the flux tube, we would have ⟨f θ ′ , E k p ⟩ → 1 for all values of θ ′ and E k p .Here, we clarify the reason for recording the proton GCRs that pass through the flux tube surface at heights above z = 0 km instead of above z = 600 km to obtain ⟨f θ ′ , E k p ⟩. First, due to the comparatively weaker magnetic field strength in the network field within the lower coronal and chromosphere regions, it is anticipated that the majority of high-energy proton GCRs will already traverse the network fields at heights far above z = 1600 km and enter the internetwork regions (Figure 1(b)).In other words, we anticipate that ⟨f θ ′ , E k p ⟩ → 1 for E k p above a certain threshold that needs to be determined.However, in our model, the injection of proton GCRs into the flux tube occurs at z = 1600 km, which is already close to the photospheric surface.At this height, the flux tube area is A top = 5.89 × 10 5 km 2 .Now, if we choose z = 600 km as the criteria for recording the proton GCRs passing the flux tube surface, the flux tube area at this height is 1.96 × 10 5 km 2 ≈ 33% × A top .This choice would not result in ⟨f θ ′ , E k p ⟩ → 1 because proton GCRs with initial injection polar angle θ 0 ≳ 155 • have not yet passed through the flux tube surface for such a short vertical distance, as from z = 1600 km to z = 600 km.(Note that in the ideal scenario where the injection started at heights within the chromosphere or lower coronal region, these proton GCRs with θ 0 ≳ 155 • would have already passed through the network field and entered the internetwork regions.)To circumvent this issue, we choose z = 0 as the criterion for recording the proton GCRs passing the tube surface.At z = 0 km, the flux tube radius is 2.01 × 10 4 km 2 ≈ 3.4% × A top .For this choice, the proton GCRs injected at θ 0 close to 180 • would have enough vertical distance, from z = 1600 km to z = 0 km, to escape the flux tube, resulting in ⟨f θ ′ , E k p ⟩ → 1 for large E k p . Figure 5 shows the numerical results of ⟨f θ ′ , E k p ⟩.Each thin colored line represents the polar angle θ ′ of a proton GCR upon exiting the flux tube surface.At E k p ≲ 10 2 GeV, we find ⟨f θ ′ , E k p ⟩ → 0 because protons injected into the tube are magnetically reflected and do not pass through the tube surface.At E k p ≳ 10 5 GeV, we find ⟨f θ ′ , E k p ⟩ → 1 for all θ ′ because protons have such large momenta that the trajectories are not affected by the magnetic fields in the flux tube.At 10 2 GeV ≲ E k p ≲ 10 5 GeV, protons are deflected by the magnetic fields and exit the tube surface with smaller polar angles.Because of the deflection, a fraction of proton GCRs entering the tube with different r 0 , θ 0 , and ϕ 0 can exit the tube surface with nearly the same θ ′ , leading to certain lines in Figure 5 showing ⟨f θ ′ , E k p ⟩ > 1.In Figure 5, the red thick line represents the average of ⟨f θ ′ , E k p ⟩ over θ ′ ranging from 90 • to 180 which follows ⟨f E k p ⟩ ≤ 1.We can understand ⟨f E k p ⟩ as the fraction of the total numbers of injected protons with E k p capable of passing through the flux tube surface.It reveals that within the range of 10 2 GeV ≲ E k p ≲ 10 3 GeV, which corresponds to the rising part of the red line, a fraction of the injected protons traverse through the flux tube, albeit with a slight angular deflection.For E k p ≳ 10 3 GeV, protons injected into the flux tube are not captured by the flux tube, but all pass through the tube surface.
Last, we assume that proton GCRs exiting the tube surface follow straight-line trajectories (i.e., no magnetic fields) until they reach the top of the flux sheet.As a result, we have ⟨f θ ′ , E k p ⟩ → ⟨f θ0 , E k p ⟩.

Gamma-ray Energy Spectra from pp Collisions
In this subsection, we follow the methodology of Kelner et al. ( 2006) for calculating the gamma-ray energy spectra from pp collisions.They provide an analytical form for the high-energy regime (E γ ≥ 100 GeV) and a δ-function approximation approach for the low-energy regime (1 GeV ≲ E γ ≤ 100 GeV), where they note that this approximation is more accurate in this energy range.

Higher-energy Regime
For E γ ≥ 100 GeV, we utilize the analytical expression of F γ , the gamma-ray energy spectrum resulting from pp interactions, as provided by Kelner et al. (2006) in their Equation (58).First, we consider the interaction of a proton GCR with a proton in the solar gas, where the proton GCR has total energy E p = E k p + m p c 2 .The number of gamma-ray photons per pp interaction within the energy interval [E γ , E γ + dE γ ] is given by The analytical expression of F γ provided in Kelner et al. (2006) is a parameterization of the gamma-ray energy spectrum obtained from the numerical simulation performed with the SIBYLL code (Fletcher et al. 1994), which takes into account inelastic pp interactions and the subsequent decays of the secondary π 0 and other mesons into gamma rays.They have shown that their result maintains an accuracy within a few percent under conditions where E γ /E p ≥ 10 −3 and E p > 100 GeV.Because of the second condition, we will use F γ in their Equation ( 58) to calculate the solar disk gamma-ray emission spectrum, dN γ /dE γ , from the solar disk for E γ ≥ 100 GeV.

Lower-energy Regime
For 1 GeV ≤ E γ ≤ 100 GeV, we utilize the δ-function approximation provided in Kelner et al. (2006) in their Equations ( 77) and ( 78).We show their methodology of the δ-function approximation in the same form as in Equation ( 15), where G γ , denoting the gamma-ray energy spectrum produced in pp interactions under the assumption of the δ-function approximation, is presented as Here, the parameter n f is selected as a free parameter to match dN γ /dE γ from the δ-function approximation for 1 GeV ≤ E γ ≤ 100 GeV and dN γ /dE γ from the analytical expression for E γ ≥ 100 GeV at E γ = 100 GeV.The parameter K π represents the mean fraction of E k p transferred to the secondary π 0 and other mesons.The value of K π = 0.17 provides a satisfactory agreement with Monte Carlo simulations based on laboratory data, as discussed in Kelner et al. (2006) and Aharonian & Atoyan (2000).The δ-function approximation is valid when σ pp is nearly constant, which occurs for E p ≳ 3E th = 3.66 GeV where E th = m p + m π + m 2 π /2m p c 2 = 1.22 GeV is the threshold energy of π 0 production with m π being the rest mass of π 0 (Kelner et al. 2006).In this work, while we present our numerical results for dN γ /dE γ down to E γ = 1 GeV, we consider the validity of the results to be limited to E γ ≳ 3.66 GeV.
Last, F γ from the analytical expression and G γ from the δ-function approximation method given in Kelner et al. (2006) do not cover the angular profile of secondary gamma rays from pp interactions.As a result, F γ and G γ should be used under the assumption that secondary gamma rays are collinear with the primary proton.This assumption is applied for calculating solar disk gamma-ray emission in Section 5.3.In Appendix C, we demonstrate that considering a nonzero emission angle would only enhance the total solar disk gamma-ray spectrum by ≈ 3% at E γ ∼ 1 GeV.Therefore, the assumption of collinearity is adequate for this work.

Gamma-ray Emission for Higher-energy Regime
In this subsection, we calculate solar disk gamma-ray spectra from the flux tube and the flux sheet for the higher-energy regime (E γ ≥ 100 GeV).
Flux tube: For proton GCRs injected into the tube ("tb") at an axial distance r 0 and the injection height z = 1600 km, the resulting gamma-ray flux in the gamma-ray energy interval [E γ , E γ + dE γ ] is given by p , E k p + dE k p ] at the injection site of the flux tube.We use the GCR flux at 1 au measured by the AMS-02 and CREAM (see Section 3.1).The factor cos θ 0 accounts for the effective flux of proton GCRs entering the horizontal tube cross-section surface at the injection height (z = 1600 km), and dΩ 0 = sin θ 0 dθ 0 dϕ 0 is the differential solid angle for injection into the flux tube.The factor dE p /E p is a product of 1/E p from Equation ( 15) and the differential integration variable dE p used to integrate Φ p .For the upper limit on E p , we use 100 E γ .The contribution from E p > 100 E γ to the gamma-ray flux is less than 1% and can be neglected.
The function S p represents the integrated absorption probability along the particle trajectory for a proton GCR injected into the tube with initial conditions r 0 , θ 0 , ϕ 0 , E k p , and is furthermore weighted by the gamma-ray transmission factor, ζ (r), for a gamma ray produced at location r, thereby taking into account the fraction of gamma rays transmitted outward from the Sun.We express S p as (19) Here, χ p = χ p (r) is the integrated column density of a proton GCR along its trajectory from the injection site r inj to the location r in hydrogen gas, where m H is the mass of a hydrogen atom.The upper bound of the integral, χ max p , in Equation ( 19) denotes the total integrated column density throughout the particle trajectory, from the injection site r inj to the exit point r exit , i.e., χ max p = χ p (r exit ).The function P abs χ p , E k p in Equation ( 19) is the absorption probability of a proton GCR with the accumulated column Li et al.
where σ inel (E p ) = σ inel E k p + m p c 2 is the cross-section of the inelastic pp interaction between the proton GCR and a hydrogen atom in the solar atmosphere.We use the numerical fit of σ inel (E p ) presented in Kelner et al. (2006) in their Equation ( 73), where D = ln E p /10 3 GeV and 1 mb = 10 −27 cm 2 .The gamma-ray transmission factor ζ (r) introduced in Equation ( 19) denotes the probability for a gamma ray photon produced at r from the pp interaction being transmitted through the solar gas.It is expressed as where λ γ is the photon mass attenuation length.For E γ ≳ 1 GeV, λ γ is approximately 80 g/cm 2 for both hydrogen and helium gases (Workman et al. 2022).Here, t γ (r) is the mass column density of the gamma-ray photon produced at r, which propagates to a far distance, r infty , from the Sun, i.e., In our numerical calculations of t γ , we have incorporated the curvature of the Sun and the z-dependence of ρ.We assume that the gamma-ray momentum is collinear with the primary proton GCR momentum during each pp interaction.Consequently, the gamma-ray trajectory used in t γ (r) in Equation ( 24) is a straight line beginning at location r and is collinear with the velocity vector of the primary proton GCR at location r.In Appendix C, we discuss the validity of the collinearity assumption.
The gamma-ray flux in Equation ( 18) only represents the case where proton GCRs are injected at r 0 .To factor in the cross-section surface area at the injection height, Equation ( 18) is weighted over the ring area 2πr 0 dr 0 /A tot .It is expressed as where A tot = πR top 2 .Flux sheet: For proton GCRs injected into the sheet ("sh") at distance y 0 at the injection height z = 600 km, the resulting gamma-ray flux in the energy interval where much is as above in Equation ( 18).The component ⟨f θ0 , E k p ⟩ is the angular and energy efficiency of protons injected into the flux sheet, as shown in Figure 5.It accounts for the effect that those high-energy protons not confined by the flux tube magnetic fields can traverse through the flux tube and subsequently enter the flux sheet.The component W sh /L takes into account the ratio of the flux sheet cross-section area at z = 600 km to the area between the two adjacent granule lanes, as discussed in the model schematic diagram in Section 3.5.We take the separation of the granule lanes to be the mean size of a granule, L = 1200 km.
To factor in the cross-section surface of the flux sheet at the injection height, Equation ( 26) is weighted over dy 0 /y top .It is expressed as Finally, the averaged gamma-ray flux evaluated at the solar surface is the sum of gamma-ray fluxes from the flux tube in Equation ( 25) and from the flux sheet in Equation ( 27).By "averaged," we mean that all patches on the surface are considered equivalent due to our assumption of the spherical symmetry of the solar surface.As a result, the averaged gamma-ray flux evaluated at 1 au from the Sun is given as

Gamma-ray Emission for the Lower-energy Regime
In this subsection, we present the formulae for solar disk gamma-ray spectra from the flux tube and the flux sheet in the low-energy regime (1 GeV ≲ E γ ≤ 100 GeV).
Flux tube: For proton GCRs injected into the flux tube, the gamma-ray flux is given by where Kelner et al. 2006).The free parameter n f in G r in Equation ( 29) for the lowenergy regime needs to be adjusted to 2.71 to match dN γ,tb (r 0 )/dE γ | R⊙ in Equation ( 18) for the high-energy regime.
Flux sheet: For proton GCRs injected into the flux sheet, the gamma-ray flux is given by The free parameter n f in G r in Equation ( 30) for the low-energy regime needs to be adjusted to 1.23 to match dN γ,sh (y 0 )/dE γ | R⊙ in Equation ( 26) for the high-energy regime.
The methodology for calculating the averaged gammaray flux at 1 au from the Sun for the low-energy regime follows the same approach as outlined in Section 5.3 for the high-energy regime.As a result, we do not repeat the formulae here.

Nuclear Enhancement Factor
Gamma-ray production is not solely dependent on protons; helium nuclei, or alpha particles, also play a relevant role.In fact, they make up ≃ 10% of the number abundance in both GCRs and photospheric gas.To factor in the effect of helium nuclei, we utilize the result of nuclear enhancement factor ε M from Kachelriess et al. (2014), which gives results for a variety of cases.Following Zhou et al. (2017), which makes choices based on recent cosmic-ray data, we take ε M = 1.8.This means that the total gamma-ray spectrum should be that much larger than the one calculated using only the proton densities in the Sun and the GCRs.
6. PREDICTED GAMMA-RAY SPECTRUM Figure 6 shows our numerical results of gamma-ray energy spectra, E γ 2 dN γ /dE γ , over an energy range of 1 GeV to 2 × 10 3 GeV.The green line is our model prediction of the gamma-ray energy spectrum from the flux tube, and the magenta is from the flux sheet.Both spectra result from proton GCR interactions with hydrogen gas.For E γ ≲ 10 GeV, the gamma rays from the flux tube dominate, while the gamma rays from the flux sheet dominate at E γ ≳ few × 10 GeV.This shift arises because lower-energy proton GCRs (E k p ≲ 10 2 GeV) are magnetically confined within the flux tube.In contrast, the higher-energy proton GCRs (E k p ≳ few × 10 2 GeV) are capable of traversing through the flux tube and entering the flux sheet, as indicated from the angular and energy efficiency in Figure 5. Furthermore, due to the higher magnetic field strength at the injection height of the flux sheet compared to the flux tube, particle reflection is more effective, consequently leading to increased gamma-ray yields in the flux sheet.
In Figure 6, the gray line is the combined gamma-ray energy spectrum summing the proton interactions in the flux tube and flux sheet.To account for the helium GCR interactions discussed in Section 5.5, we multiply the gray line by the nuclear enhancement factor, ε M = 1.8, resulting in the black line.Consequently, the black line represents the total gamma-ray energy spectrum from the solar disk.The black line shows a modest spectral slope, dN γ /dE γ ∼ E −2.4 γ , for 1 GeV ≲ E γ ≲ 10 2 GeV, aligning close to the spectral shape from Fermi-LAT observations over the entire solar cycle.This moderate spectral slope is a product of two combined emission sources: the gamma-ray emission from the flux tube for E γ ≲ 10 GeV and the emission from the flux sheet for 10 GeV ≲ E γ ≲ 10 2 GeV.At much higher E γ , the spectrum steepens to dN γ /dE γ ∼ E −3.6 γ , agreeing with the spectral shape from the HAWC observation.This steep slope is due to the limited effectiveness of GCR capture and reflection within the finite-sized flux sheet -an inherent characteristic of the intermittent field distributions generated by granule convective flows.Additionally, the normalization of the black line is reasonably consistent with the observational data, remaining within a factor of about two, ranging from 1 GeV to 10 3 GeV.This agreement strongly suggests that the observed gamma-ray spectra are influenced and shaped by the flux tubes in the network elements and the flux sheets in the intergranule lanes.
Based on the discussion above, important conclusions and insights can be drawn from Figure 6.
1.In a realistic solar environment, the observed gamma-ray spectra in the ∼ 1-100 GeV range are likely shaped by the combination of multiple finite-sized magnetic flux structures at the photosphere.We note that a simple P to B scaling relation would not reproduce the observed spectra over a wide range of E γ .
2. Taking into account the finite-sized magnetic flux structures is key to the considerably softer gammaray spectra near 10 3 GeV from HAWC.This is due to the large fraction of high-energy proton GCRs that plunge through both the flux tube and flux sheet.Figure 6.Gamma-ray spectrum of the solar disk.Green and magenta lines are our calculated gamma-ray spectra for the flux tube and sheet due to proton GCR interactions, with the gray line denoting the combined flux.Black line is an 80% increase over the gray line, factoring in helium GCR interactions.Vertical dashed line at 3.66 GeV is the lowest Eγ for this work to be valid.Also presented are Fermi-LAT data from the solar minimum (2008-2010) (Abdo et al. 2011) and over the full cycle (2008-2020) (Linden et al. 2022), in addition to HAWC data from the solar maximum ("max", 2014-2017) and minimum ("min", 2018-2021) (Albert et al. 2023).Blue line and the shaded band are the best-fit and statistical uncertainties of HAWC's 6. 1-year (2014-2021) data.The band represents a single energy bin of 0.5-2.6TeV; the blue data points indicate how its height varies over the solar cycle.
3. While the data can be reasonably interpreted through these two basic structures, incorporating additional flux structures and field features could potentially enhance the accuracy of the model.This could include variability in size and average magnetic field strength of flux tubes and flux sheets, non-vertical flux structures, and the magnetic turbulence caused by the convective flow of the granules.Taking these into account could lead to a larger predicted flux in the 10 3 GeV range.
Building on prior theoretical work on the solar-disk gamma-ray emission (Seckel et al. 1991(Seckel et al. , 1992;;Zhou et al. 2017;Hudson et al. 2020;Mazziotta et al. 2020;Li et al. 2020;Gutiérrez & Masip 2020;Gutiérrez et al. 2022), our model goes much further by considering magnetic fields motivated by solar physics data and magnetoconvection simulations.Our approach introduces the concept of finite-sized flux structures into the problem, differentiating between higher-and lower-energy proton GCR behaviors by considering two separate flux structures, calculating 3D particle trajectories, and comparing theoretical predictions to data across three orders of magnitude in gamma-ray energy.
How do our results compare to those of previous work?First, the pioneering work of Seckel et al. (1991), which focused on magnetic mirroring in the magnetic canopy fields, produced predictions that disagreed with measurements by a factor of about 5-10.Moreover, their spectrum was predicted only for gamma-ray energies up to 5 GeV.Several alternative approaches have been proposed in recent years.For instance, Mazziotta et al. (2020) and Li et al. (2020) focus on the effects of coronal fields, which extend towards the photosphere.However, their results are about three times smaller than the observational data at around 10 GeV and 10 times smaller at 100 GeV.Most importantly, these models do not align with the new HAWC data in terms of normalization and spectral shape at about 10 3 GeV.Nevertheless, their findings point to the potential significance of solar coronal fields.
For future work, in Figure 6, while we show our theory predictions down to E γ = 1 GeV, we note that they are considered valid only for E γ ≳ 3.66 GeV.This is because the δ-function approximation provided in Kelner et al. (2006) is only valid for E p ≳ 3E th = 3.66 GeV where σ pp is nearly constant.Although the method for E γ < 3.66 GeV is currently absent in the literature, this task could potentially be performed using publicly avail-able Monte Carlo codes, for instance, FLUKA (Battistoni et al. 2015;Ahdida et al. 2022) and GEANT4 (Agostinelli et al. 2003).Further investigation in this aspect is an interesting direction for future work.
Last, in Appendix A, we present the most probable injection polar angles for GCRs that produce the escaping gamma rays from the solar disk as ≈ 171 • ± 5 • for the flux tube and ≈ 135 • ± 15 • for the flux sheet.In Appendix B, we perform the averaged emission angle and height of the gamma-ray flux.We find that the emission primarily occurs within a height ranging from −100 km to 500 km, encompassing the photosphere and extending ∼ 100 km into the uppermost convection zone.Moreover, we find an emission angle of ≈ 65 • ± 15 • relative to the normal direction of the solar surface, indicating that the gamma rays predicted by our model predominantly originate from the outer regions of the solar disk.This range of predicted angles is a result of our exclusive focus on vertical flux structures and the omission of magnetic turbulence in our model.In a more realistic solar surface environment, the presence of non-vertical flux structures with magnetic turbulence would have an impact on particle trajectories, causing deviations from the magnetic bottle effect.This can potentially lead to smaller emission angles and a more uniform distribution of gamma-ray emission across the solar disk.We will explore this in future work.

DISCUSSION AND CONCLUSIONS
In this paper, we present a simple model of solarsurface magnetic fields, aiming to understand the broad outline of how hadronic GCRs shape gamma-ray emission.We have investigated GCR propagation and gamma-ray emission within a simplified solar-surface magnetic field model comprising a flux tube and a flux sheet.In the ∼ 1-100 GeV range, where emission is due to both the tube and the sheet, our model produces a hard spectrum (dN γ /dE γ ∼ E −2.4 γ ).In the ∼ 100-1000 GeV range, where the tail of the sheet emission dominates, our model produces a considerably softer spectrum (dN γ /dE γ ∼ E −3.6 γ ) due to the limited effectiveness of GCR capture and reflection by the finite-sized flux sheet.Despite the fact that the model presented in this work does not capture all the complexities of the solar-surface magnetic fields, the reasonable agreement between our model predictions and data (within a factor of 2 over about three orders of magnitude in gamma-ray energy) is encouraging and highlights the critical roles of finite-sized field structures of flux tubes and flux sheets in shaping solar-disk gammaray emission.
We emphasize that the model, as presented in this study, remains a work in progress and does not yet offer a conclusive solution to the long-standing issue of solardisk gamma-ray emission.Given the multi-scale nature of the solar magnetic field, which spans from large-scale heliospheric magnetic fields to small-scale Alfvén wave turbulence and dissipation, incorporating a realistic solar magnetic field model into the current problem is a challenging task and necessitates a stepwise approach.The reasonable agreement with data suggests that the double-structure model with finite-sized magnetic flux structure captures key characteristics of the hadronic GCR reflection at the solar surface and the resulting solar gamma-ray emission.This simplified approach sets the stage for further refined modeling and advanced computational efforts, details of which are elaborated as follows.
In future work, it will be important to go beyond the assumptions here of a simple double-structure model in the quiet photosphere with an assumption of the spherical symmetry.Aspects to investigate include the influences of active regions and coronal-hole open field distributions in GCR transport and gamma-ray emission.Ultimately, we need to evaluate the impact of active regions on GCR transport toward the solar surface.We also need to quantify the relationship between the latitudinal dependence of coronal-hole distribution and the resulting gamma-ray emission.Both directions hold great potential in unraveling the observed anti-correlation between solar activity and the gamma-ray flux (Ng et al. 2016;Tang et al. 2018;Linden et al. 2018Linden et al. , 2022;;Albert et al. 2023).
Future studies should also explore GCR propagation and solar gamma-ray emission within the framework of magnetoconvection in the photosphere and uppermost convection zone.In this work, we focused on the vertical magnetic flux tube and flux sheet in magnetohydrostatic equilibrium with the surrounding gas.The perpendicular field components leading to particle reflection are caused by the compression exerted by the surrounding gas pressure.Eventually, we need to understand whether and how magnetic turbulence resulting from the convection of the granule cells can cause a similar effect and affect GCR propagation.Investigating the impact of turbulence may further elucidate additional magnetic field structures at the photosphere that contribute to the observed gamma-ray spectra.
Solar gamma-ray observations have presented an intriguing new opportunity to investigate the characteristics of magnetic fields in the photosphere.The origin of these magnetic fields, whether arising from the decay of active regions (Spruit et al. 1987) or local fast Li et al. dynamo actions in the uppermost convection zone (Cattaneo 1999;Cattaneo et al. 2003;Vögler & Schüssler 2007), remains uncertain.By refining theoretical models of solar magnetic fields and the transport of GCRs, we have the potential to gain valuable insights into this unresolved matter.In this regard, future research should focus on exploring these issues through magnetoconvection simulations and utilizing high-resolution measurements of small-scale magnetic structures at the solar surface, such as those provided by the Daniel K. Inouye Solar Telescope (Rimmele et al. 2020).
We are grateful for stimulating discussions with Ofer Cohen, Federico Fraschetti, Hugh Hudson, Mikhail Malkov, and Eleonora Puzzoni.This work was supported by NASA grant Nos.80NSSC20K1354 and 80NSSC22K0040.J.F.B. was additionally supported by National Science Foundation grant No. PHY-2012955.J.T.L. thanks the Lunar and Planetary Laboratory for their hospitality while parts of this manuscript were completed.

A. MOST PROBABLE GCR POLAR INJECTION ANGLES
In this appendix, we show the most probable injection polar angles for GCRs that produce the gamma rays that escape from the solar disk.
For proton GCR injection into the flux tube, we average S p in Equation ( 19) over the injection area, 2πr 0 dr 0 , and the azimuthal injection angle, ϕ 0 , (A1) The term ⟨S p ⟩ tb is interpreted as the averaged absorption probability for a primary proton GCR along the section of its trajectory where the produced gamma rays from pp interactions successfully escape without being absorbed by the Sun.
For proton GCR injection into the flux sheet, S p is averaged over the injection width dy 0 and the azimuthal injection angle dϕ 0 .It is expressed as ) where ⟨f θ0 , E k p ⟩ accounts for the angular and energy efficiency of the injected particles into the flux sheet.Note that Equation A2 accounts for monoenergetic protons with a kinetic energy of E k p , without incorporating any weighting based on the proton GCR spectrum.
Figure 7(a) shows our calculations of ⟨S p ⟩ tb in Equation (A1) for proton GCR injection into the flux tube.The three lowest E k p lines indicate that the most probable polar injection angle is θ the injected protons are reflected too high from the solar surface to accumulate significant column density.For θ 0 ≳ 175 • , the injected protons are not yet reflected by the magnetic fields by the time they reach the bottom of the simulation box; thus, they do not contribute outward-directed gamma rays.
In Figure 7(a), the three lowest E k p have approximately equivalent ⟨S p ⟩ tb .In contrast, for the highest E k p , then ⟨S p ⟩ tb has substantially smaller values across all most probable θ0 .This difference arises because protons with E k p ≲ 100 GeV are nearly all reflected within the flux tube via the magnetic bottle effect.In contrast, protons with E k p ≳ 10 3 GeV do not spiral within the flux tube but pass through the tube surface with only small angular deflection.This result echos the findings from the angular and energy efficiency of GCRs injected into the flux sheet, ⟨f θ0 , E k p ⟩, as presented in Figure 5. p has lower ⟨S p ⟩ sh .This is because the majority of the proton GCRs are magnetically reflected by the flux tube, while only a small fraction of protons successfully penetrate through the flux tube and enter the flux sheet.The highest E k p also has lower ⟨S p ⟩ sh .This is because the efficiency of capturing and reflecting proton GCRs by flux sheet magnetic fields starts diminishing for E k p ≳ few × 10 3 GeV.Finally, for E k p ≳ 10 4 GeV, ⟨S p ⟩ sh ≈ 0 across all θ0 , and are not shown in Figure 7.
In Figure 7(b), all four energies show ⟨S p ⟩ sh ≈ 0 for θ0 ≳ 160 • .This result implies that protons injected into the flux sheet with θ0 ≳ 160 • do not experience magnetic reflection before being absorbed by the Sun.Therefore, any gamma rays generated from these protons do not contribute to the outward-directed gamma rays.Note that the range of these polar injection angles is significantly broader compared to the flux tube scenario depicted in Figure 7(a).

B. AVERAGE EMISSION ANGLE AND HEIGHT
In this appendix, we present the average emission angle and height of the solar-disk gamma rays.We only present the methodology for the higher-energy regime (E γ ≥ 10 2 GeV).The methodology for the lower-energy regime (1 GeV ≲ E γ ≤ 10 2 GeV), which has been discussed in Section 5.4, is not repeated here.

B.1. Average Emission Angle
Flux tube: First, we calculate the average gamma-ray emission angle in the case of a flux tube.We define θ p (r) as the angle of the proton's velocity vector at location r relative to the normal direction of the surface.Thus, for proton GCRs injecting into the flux tube at axial distance r 0 at the injection height z = 1600 km, the θ p (r)-weighted gamma-ray flux is given as (B3) Here, S θp p is defined as S p in Equation ( 19) weighted by θ p (r), i.e., As a result, the average emission angle ⟨θ p ⟩ for the gamma-ray flux at E γ is given as where dN γ,tb /dE γ | R⊙ is provided in the calculation of solar-disk gamma-ray flux in Equation ( 25).Next, we calculate the root mean square (RMS) emission angle, θ p (r) − ⟨θ p ⟩, for the gamma-ray flux at E γ .This time, the gamma-ray flux is weighted by the mean squared polar angle, (θ p (r) − ⟨θ p ⟩) 2 , of each proton's moving direction at location r.It is given as  (B10) Flux sheet: The methodology for the average and RMS emission angles in the case of the flux sheet is the same as in the flux tube.The only difference is the integrants in Equations ( B3) and (B7) should be further multiplied by the angular and energy efficiency ⟨f θ0 , E k p ⟩ shown in Figure 5. Consequently, we do not repeat the formulae here.
Figure 8 shows the average (blue solid lines) and RMS (blue shaded bands) emission angles in the cases for the flux tube (top) and flux sheet (bottom) for gamma rays successfully transmitted from the Sun.Our result indicates that for each small patch at the surface of the Sun, the average emission angle is ≈ 65 • , relative to the normal direction, with a corresponding RMS emission angle of approximately 15 • .Finally, gamma rays are fully absorbed by the Sun whenever θ p (r) ≥ θ p, crit ≈ 93 • , which are highlighted in red shaded bands in Figure 8.

B.2. Average Emission Height
The methodology employed to calculate the average and RMS emission heights (denoted as ⟨z p ⟩ and ⟨∆z p ⟩, respectively) are consistent with the approach described in Appendix B.1 for computing emission angles.Therefore, we do not reiterate the formulas here.
Figure 9 shows ⟨z p ⟩ and ⟨∆z p ⟩ for the gamma rays produced in the flux tube (top) and flux sheet (bottom).Our findings indicate that emission typically occurs at the height of z ≈ 0 km with an RMS of ≈ 300 km for the flux tube and with an RMS of ≈ 150 km for the flux sheet.Consequently, this suggests that the layers responsible for producing the gamma rays successfully transmitted from the Sun occur in the photosphere and a thin layer extending ∼ 100 km into the upper convection zone.

C. IMPACT OF NON-COLLINEAR EMISSION ON GAMMA-RAY FLUX
In this appendix, we estimate the errors incurred by assuming that gamma rays are collinear with their parent protons.
In this simplified model, we consider a gamma-ray emission cone with a half-angle denoted as ∆Θ γ .Accordingly, the solid angle, Ω γ , of the emission cone is expressed as 2π (1 − cos ∆Θ γ ).We assume that the number of gamma-ray photons distribute uniformly in the angular direction in Ω γ .Under this assumption, we conduct a numerical calculation of S p in Equation ( 19) to evaluate the impact of a nonzero emission cone angle on the solar gamma-ray flux.(To remind the reader, S p in Equation ( 19) is the total proton GCR absorption probability along a particle's trajectory and is further weighted by the gamma-ray transmission probability, as discussed in Section 5.3.) We consider a proton GCR injecting into the flux tube at θ 0 = 171 • because this is where ⟨S p ⟩ tb peaks (see Figure 7  and E k p = 15.4GeV has nearly the same ⟨S p ⟩ tb .Thus, the choice of E k p = 15.4GeV is sufficient for our purpose.Last, when E k p is as low as 15.4 GeV, the selections of ϕ 0 and r 0 have minimal influence on ⟨S p ⟩ tb as they yield similar particle trajectories.Thus, we opt for ϕ 0 = 180 • , and r 0 = 20 km.Our chosen set of parameters r 0 , θ 0 , ϕ 0 , E k p for proton GCR injection into flux tube should aptly represent the peak S p .Therefore, the changes of S p observed for this particular set of injection parameters can be regarded as an upper limit on the true impact on the solar gamma-ray emission.
Figure 10 presents our results for S p with different angular sizes of emission cones.The x axis is the polar angle of the primary proton GCR along its trajectory.The y axis is S p as a function of bins of θ p and is further normalized by ΣS p, collinear , which is the summation of S p over all bins of θ p in the case of collinearity.The black line represents S p /ΣS p, collinear in the case of collinearity.Therefore, the area under the black line equates to 100%, as displayed in the top-right corner of the diagram.The black line shows that the gamma rays successfully transmitted from the Sun are produced at θ p < 90 • , echoing the findings in the average angular emission plot in Figure 8.
In Figure 10, color lines show emission cones with ∆Θ γ of 10 • , 20 • , and 30 • .All three color lines indicate an increase in S p when θ p > θ p,crit ≈ 93 • .This increase is due to a fraction of gamma rays leaving the Sun before θ p reaches θ p,crit .In contrast, when θ p < θ p,crit , we see a drop in S p compared to the black line.This decrease is due to a fraction of gamma rays remaining fully absorbed by the solar gas at this stage.
In Figure 10, the three colored percentages (100.8%,102.9%, and 106.3%) correspond to the sum of S p across all θ p bins, normalized by ΣS p, collinear , for ∆Θ γ values of 10 • , 20 • , and 30 • respectively.These percentages are over 100%, which implies that considering a finite-sized emission cone leads to a slightly higher total gamma-ray flux from the solar disk than the collinearity scenario.However, this S p increase is about a few percent and should only increase the solar-disk gamma-ray flux by about the same order of magnitude.
Last, to estimate the range of gamma-ray emission angles resulting from pp interactions, we utilize the particle shower Monte Carlo simulation code FLUKA (version 4-2.2) (Battistoni et al. 2015;Ahdida et al. 2022) and its graphical user interface Flair (version 3.2-4.5)(Vlachoudis 2009).We inject a proton into a slab consisting of equal numbers of protons and electrons.We calculate the required ∆Θ γ that contains either 90% (solid line) or 68% (dashed line) of the total number of gamma rays with E γ ≥ 1 GeV.
Figure 11 shows the result from our FLUKA simulation.The result suggests that 90% of the gamma rays are contained within ∆Θ γ = 20 • for E k p = 10 GeV, roughly equivalent to peak gamma-ray production for E γ ∼ 1 GeV.Based on these findings and the S p result of ∆Θ γ = 20 • presented in Figure 10, we conclude that the collinearity assumption should not cause the actual flux to deviate by more than ≈ 3%.

Figure 1 .
Figure 1.Schematic diagrams (not to scale) of the magnetic network fields and coronal-hole open field lines, with the viewing area progressively expanding from (a) to (c).(a) Bird's eye view of solar surface: Vertical magnetic flux tubes arise from granules (white polygons) that reside in the intergranule lanes (gray lanes) along the edges of the supergranule cells.A collection of flux tubes merges at a height of approximately 600-1000 km, forming a network element.(b) Vertical 2D slice of the solar surface: The network fields expand laterally and merge with others.The regions between the network elements are filled with granules (small solid arcs with arrows), with magnetic loops (dashed arcs) formed near each granule.(c) Solar corona and interplanetary magnetic fields: Lines with one footpoint are the network fields in the coronal-hole regions, forming the open magnetic field lines that extend into interplanetary space.Lines with two footpoints are the closed loops.(a) and (b) are redrawn illustrations based on Cranmer & van Ballegooijen (2005) and Wedemeyer-Böhm et al. (2009), respectively.

Figure 2 .
Figure 2. Schematics of the double-structure model.Left:Lower-energy hadronic GCRs (black helix) are magnetically mirrored inside the flux tube, producing outgoing gamma rays (red wavy line).Higher-energy hadronic GCRs (blue curve) plunge through the flux tube with a slight deflection.Right: Higher-energy hadronic GCRs (blue helix) then enter the flux sheet, with a fraction being magnetically reflected, producing outgoing gamma rays.

Figure 3 .
Figure 3. Adopted model conditions for the upper convection zone, photosphere, and chromosphere.Solid lines are data from the HSRASP model.Dashed lines are our extrapolation.

Figure 4 .
Figure 4. (a) For the flux tube case: vertical cross-section plane and magnetic field strength variations.The black line denotes the edge of the flux tube.At each height (labeled by a horizontal gray line), the r-direction variations of Bz and Br, normalized to the value of Bz at the axis (with their numerical values labeled in red), are shown in the red solid and blue dashed lines, respectively.(b) For the flux sheet case.(Note the vertical scale in this representation has been compressed by factors of 2.8 and 2.25 in relation to the horizontal scales depicted in (a) and (b), respectively.The actual dimensions of the flux tube and sheet are more slender than the depicted structures here.)

Figure 5 .
Figure 5. Angular and energy efficiency, ⟨f θ ′ , E k p ⟩, of proton GCRs passing through the flux tube surface.The color bar corresponds to the polar angle, θ ′ , at which the proton GCR exits the flux tube surface.
18)where the subscript R ⊙ denotes this flux is evaluated at the solar surface.The function Φ p E k p represents the differential proton GCR flux per steradian per proton kinetic energy interval[E k

Figure 7 .
Figure 7.Most probable GCR polar injection angles.The higher the values of ⟨Sp⟩, the greater the production of outward-directed gamma rays.(a) ⟨Sp⟩ tb for the flux tube shown in Equation (A1).(b) ⟨Sp⟩ sh for the flux sheet shown in Equation (A2).
Figure 7(b) shows our calculations of ⟨S p ⟩ sh in Equation (A2) for proton GCR injection into the flux sheet.Note that the lines with E k p = 1.54 GeV and 1.54 × 10 1 GeV shown earlier in Figure 7(a) are not shown in Figure 7(b) because the majority of proton GCRs at such low energy do not penetrate through the flux tube.In Figure 7(b), all four energies indicate that the most probable injection angle is 120 • ≲ θ0 ≲ 150 • .The two intermediate E k p show high ⟨S p ⟩ sh , indicating that those protons are effectively reflected and produce outwarddirected gamma rays.Notably, these two intermediate energies lie within the rising part of ⟨f E k p ⟩ in Figure 5.In contrast, the lowest E k B4)The θ p -weighted gamma-ray flux in Equation (B3) is further weighted over the area along the cross-section surface of the flux tube at the injection height,

Figure 8 .
Figure 8.Average (blue lines) and RMS (blue bands) emission angles for gamma rays successfully transmitted from the Sun.
(a)).As for injecting E k p , we choose 15.4 GeV because proton GCR at this kinetic energy should produce most of the gamma-ray flux at E γ ∼ 1 GeV.(Note: gamma rays of energy E γ are roughly produced by proton GCRs of kinetic energy E k p ∼ 10 × E γ .)In addition, Figure7(a) has also shown that E k p = 1.54 GeV

Figure 9 .
Figure 9. Average (black lines) and RMS (blue bands) emission heights for gamma rays successfully transmitted from the Sun.
Figure 10.Proton GCR absorption probability in the flux tube weighted by gamma-ray transmission probability under non-collinear gamma-ray emission.The injection polar angle θ0 is 171 • and proton GCR E k p is 15.4 GeV.