Cosmic Ray transport in mixed magnetic fields and their role on the observed anisotropies
Margot Fitz Axen, Julia Speicher, Aimee Hungerford, Chris L. Fryer
MMNRAS , 1–14 (2020) Preprint 20 January 2021 Compiled using MNRAS L A TEX style file v3.0
Cosmic Ray Transport in Mixed Magnetic Fields and their role onthe Observed Anisotropies
Margot Fitz Axen, , (cid:63) Julia Speicher, Aimee Hungerford and Chris L. Fryer Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, Texas 78712-1205, USA American Astronomical Society, 2000 Florida Ave., NW, Suite 300, Washington, DC 20009-1231, USA Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos, NM 87544, USA
Accepted 2020 November 5. Received 2020 October 29; in original form 2020 May 13
ABSTRACT
There is a growing set of observational data demonstrating that cosmic rays exhibit small-scaleanisotropies (5-30 ◦ ) with amplitude deviations lying between 0.01-0.1% that of the averagecosmic ray flux. A broad range of models have been proposed to explain these anisotropiesranging from finite-scale magnetic field structures to dark matter annihilation. The standarddiffusion transport methods used in cosmic ray propagation do not capture the transport physicsin a medium with finite-scale or coherent magnetic field structures. Here we present a MonteCarlo transport method, applying it to a series of finite-scale magnetic field structures todetermine the requirements of such fields in explaining the observed cosmic-ray, small-scaleanisotropies. Key words: astroparticle physics — ISM: magnetic fields — ISM: cosmic rays
Observations of cosmic rays in the TeV-PeV energy range havedemonstrated both small- and large-scale anisotropies. The large( > ◦ ) and small-scale anisotropies have been observed by a largenumber of instruments (Amenomori et al. 2005, 2006; Guillianet al. 2007; Abdo et al. 2008; Aglietta et al. 2009; Abdo et al.2009; Amenomori et al. 2010; Munakata et al. 2010; Abbasi et al.2010; Cui 2011; Abbasi et al. 2011; Aartsen et al. 2013). The large-scale dipole anisotropy is reasonably well fit by the asymmetries innearby sources smoothed by the subsequent diffusion of these cos-mic rays (Erlykin & Wolfendale 2006; Blasi & Amato 2012; Pohl& Eichler 2013; Sveshnikova et al. 2013). It is possible that mag-netic fields in the heliosphere could affect the anisotropy (Desiati &Lazarian 2013; Schwadron et al. 2014).The small-scale anisotropies are much more difficult to ex-plain. These small-scale anisotropies typically have size scalesbetween − ◦ and amplitudes between 0.01-0.1% of the back-ground cosmic ray flux (Blasi & Amato 2012). A number of mod-els have been proposed to explain such anisotropies. For example,anisotropies could arise from a nearby, as yet undetected, supernovaremnant (Salvati & Sacco 2008), perhaps mediated by a local, co-herent magnetic field or asymmetry in the propagation (Drury &Aharonian 2008; Malkov et al. 2010; Biermann et al. 2013). Anotherset of proposals argue that properties of the heliosphere can drivethe observed anisotropies (Drury & Aharonian 2008; Lazarian &Desiati 2010; Desiati & Lazarian 2013). More exotic models havealso been proposed invoking strangelet production or dark matter (cid:63) E-mail: [email protected] annihilation (Harding 2013; Kotera et al. 2013). More recently, aseries of models have proposed that turbulent magnetic fields aresufficient to explain the anisotropies (Giacinti & Sigl 2012; Ahlers2014).Magnetic fields generated in turbulence are believed to existon all scales, from the smallest scales in the Kolmogorov spectrumto the largest eddy scales(Kulsrud & Zweibel 2008; Saveliev et al.2012; Krasheninnikov 2013; Kritsuk et al. 2017; Kim et al. 2020).In the Milky Way, this leads to magnetic field structures rangingfrom well below the parsec scale up to a kiloparsec(Kulsrud &Zweibel 2008; Saveliev et al. 2012; Krasheninnikov 2013). If thescale of the magnetic fields were limited to the smallest turbulencescales, particle transport within these fields could be treated in thediffusion limit where the magnetic fields are treated as a scatteringterm with a net energy loss as is done in codes like GalProp (Stronget al. 2007). This treatment does not accurately model any possiblelarger scale, coherent structures in the magnetic field.However, transport methods that can model both small andlarge scale fields face numerical challenges. For instance, themethod described in Fryer et al. (2007) and Harding et al. (2016)allows for transport in one of two extreme solutions for each spatialzone: either dominated by small scales (isotropic scattering limit)or dominated by coherent fields. In this paper, we generalize theMonte Carlo method from Fryer et al. (2007) and Harding et al.(2016) to allow for more general magnetic field profiles that includeboth small- and large-scale features. As in Harding et al. (2016),we focus on a magnetic field configuration that studies the inter-action of a single point source for cosmic ray production with acoherent magnetic field of varying strength relative to the small-scale field (sufficiently small with respect to the spatial scale that © a r X i v : . [ a s t r o - ph . H E ] J a n Fitz Axen et al. the interaction can be treated as isotropic). With this study, we hopeto determine the magnetic field properties needed to explain theobserved anisotropies in the cosmic ray flux more realistically.In section 2, we describe the properties of our grid, particlesand magnetic field configurations. In section 3, we describe themethods we use for particle propagation and their differences fromprevious work. In section 4 we describe the verification tests we usedon the code. In section 5, we study the propagation of cosmic raysin a box to apply this code to a cosmic-ray transport application. Insection 6, we study the observational implications of these transportcalculations. Finally, in section 7, we conclude with a discussion ofthe properties needed to produce cosmic ray anisotropies.
The setup of our grid is similar to that used by Harding et al. (2016).We define a three dimensional grid of cubic zones, with physicalproperties that are constant within each zone. The grid we use is a150x150x150 pc grid, divided into 50 zones of 3x3x3 pc. Thoughthe zone size is arbitrary, it is set to be much larger than the cosmicray scattering length and Larmor radius of particle motion. Wedefine the tally plane as being the upper x face of the grid, at x=150pc. The particles are propagated from a chosen starting locationthrough this grid until they hit one of the grid faces or until they areabsorbed into the ISM.The current code physics is currently limited to the propagationof protons, which are expected to be the dominant form of cosmicrays that create the observed anisotropies in the TeV-PeV energyrange (Harding et al. 2016). We assume a point source for cosmicrays placed in the center of the simulation grid, emitting particlesin random directions according to an isotropic distribution. All par-ticles are assumed to stay traveling at relativistic velocities (v ≈ c).We tested primarily particles with energy of 10 TeV, which is atthe lower end of the energies observed for cosmic ray anisotropies.Additionally, we did one study in which we varied the cosmic rayenergy to higher energy values, up to 100 PeV. Note that the protonmean free path is what sets the fundamental scale of the computedsolutions and there is a straightforward relationship between protonenergy and proton mean free path shown in equation (6). Varia-tion of the assumed proton energy has been used as our primarymeans to study the impact of scaling the proton mean free path.In principle, these particles lose energy due to Coulomb scattering,ionization, and other processes, as was studied in Harding et al.(2016). However, for the size-scale and conditions in their, and oursimulation grid, they found that energy losses of protons were min-imal (~1-10 MeV over the course of their entire propagation) and,for the calculations in this paper, we do not include energy losses. The code is designed so that the properties of the magnetic field canbe varied in each zone. This includes the small-scale magnetic fieldamplitude and the coherent magnetic field amplitude and direction.For the calculations in this paper, we focus on a bimodal distributionof field structures: small-scale fields with size-scales well below thezone size, and coherent magnetic fields that are equal to or greaterthan the zone size. In the extreme case of no coherent magneticfield, we can transport using a statistically sampled isotropic mag-netic field direction equivalent to a symmetric diffusion coefficient. Focusing on a single coherent field allows us to test the ability ofthe code to model small- and large-scale magnetic field structurescombined in a single zone and to study in detail the effects of theseglobal structures on the transport. This setup is a considerable im-provement over the grid design used by Harding et al. (2016), whosesetup only allowed for either a turbulent field or a coherent field ineach zone of the grid.For our simulations, we assume that the entire grid is filledwith a small-scale magnetic field component of constant amplitude B t µ G. This field is assumed to be produced in the smallest tur-bulence scales and is much smaller than our simulation grid-scale.With our grid size-scales, a typical cosmic ray undergoes manypitch-angle scatterings as it transports through the grid. The timefor energy evolution is much longer than the scattering timescalefor our high-energy protons (Schlickeiser & Miller 1998). As dis-cussed above, this long energy-evolution timescale and the relativelysmall simulation grid-size means that, although the proton under-goes many "scatterings", its energy does not vary considerably asthe particle transports through the grid (Harding et al. 2016). Underthese conditions, we can mimic the cosmic ray interaction with thesmall-scale fields as isotropic and assume the energy is constantduring this period. Therefore, we can model the turbulent magneticfield through a random vector sampled at every particle step throughthe simulation: cos θ χ − , (1) φ πχ , (2)where χ i are standard deviates between 0 and 1. Converting theseexpressions to Cartesian coordinates gives a three-component, ran-dom vector.In addition to this turbulent field, we define six grid zoneboundaries that bound a region that includes a coherent magneticfield component. This single coherent field is only speculative butis meant to approximate some of the features of the broad range ofmagnetic field scales that may exist in the ISM. Inside this region,the coherent field component has a constant magnitude and directiongiven by B g x , y , z B g , , . Clearly, modeling the coherent magneticfield in this way is an oversimplification of the net effects large scalemagnetic field structures in the ISM would have, but it allows usto test the importance of these possible large-scale coherent fields.We study a variety of different configurations for the box regioncontaining the coherent component to the magnetic field. We havechosen to keep the y and z coordinates of the box zone numbersconstant at 108 pc < y,z <
120 pc. The geometry of the coherentfield region is varied using the x zone value of the box, for whichwe vary the offset from the tally plane and the extent of the box.We tested box extents of 6 and 24 pc and varied the box offset fromthe tally plane from 0 to 18 pc. This is in contrast to Harding et al.(2016), whose coherent magnetic field region spanned the entire xlength of the simulation space. A comparison of this is shown inFigure 1.The amplitude of the coherent magnetic field is varied as afraction of the maximum amplitude of the small scale field. Theexact values for the three parameters studied are shown in Table 1.
In this section, we describe the propagation of particles throughour code. The methods we use are similar to those used in Harding
MNRAS000
MNRAS000 , 1–14 (2020) onte Carlo CR Transport Figure 1.
Comparison of a few of our simulation setups to that used by Harding et al. (2016). The source of cosmic rays is placed in the center of our simulationgrid, shown by the red star. Along a rectangular path, shown by the blue box, we place a coherent magnetic field on top of the small-scale field. The amplitudeof this coherent magnetic field is varied as a fraction of the amplitude of the small-scale field. We vary the length of the coherent magnetic field and its offsetfrom the tally plane, shown by the green region. The top left plot shows Harding’s setup, running the length of the simulation space, while the other three showsome variations of the box offset and extent which we use.
Table 1.
Shows the full set of test configurations we chose. Altogether therewere 42 different coherent magnetic field configurations.Variable Tested ValuesExtent [pc] 6, 24Offset [pc] 0, 3, 6, 9, 12, 15, 18 B g / B t et al. (2016). At the beginning of the particle’s lifetime, its startinglocation is determined by the input source location and its initialdirection is sampled randomly from an isotropic distribution. Itis then propagated through the grid using one of two methods, depending on whether it is in the coherent magnetic field region ornot. In Harding et al. (2016), the methods used included a directtransport Monte Carlo method if the particle was in the coherentmagnetic field region and, if not, a Discrete Diffusion Monte Carlo(DDMC) method (Evans et al. 2003). Given the simplicity of themagnetic field configuration, such an approach is reasonable. ThisDDMC approach reproduces the results of simple diffusion trans-port schemes such as GalProp Strong et al. (2007). While we wereable to use a DDMC algorithm very similar to that used by Hardinget al. (2016) for the majority of our grid, the direct transport methoddid not allow for the magnetic field configurations that include the MNRAS , 1–14 (2020)
Fitz Axen et al. isotropic, turbulent field component in addition to the coherent fieldcomponent inside the box.To study these hybrid regions, we modified the transportmethod used by Harding et al. (2016) to account for particle prop-agation inside the coherent field region. We tested two differentmethods for this hybrid transport/diffusion Monte Carlo algorithm.Outside the coherent field region, the particles were propagated us-ing a DDMC method similar to that used by Harding et al. (2016). In3.1 we first discuss the direct transport Monte Carlo scheme whichwas used in Harding et al. (2016) and which is the basis for our hy-brid transport method. In 3.2 we discuss the Monte Carlo diffusionapproximation which we employ for most of the grid. Finally in 3.3we discuss the two methods we use for our coherent field region inorder to account for both components to the magnetic field.
A particle transport Monte Carlo code tracks the particle motion asthe particle goes through a zone and passes into the next. Particleswithin a magnetic field of strength B propagate according to theequation of motion: dvdt cr L v × ˆ B , (3)where ˆ B is the direction of the net magnetic field, v is the directionof the particles velocity, and r L is the Larmor radius of the particlemotion. This causes particles to move in a circular path of radius r L around the magnetic field line, perpendicular to it. The componentof the particle’s initial velocity parallel to the magnetic field linedoes not change. Particles within a magnetic field of strength B ,energy E , and charge Ze have a Larmor radius of (Schlickeiser2002): r L . × − Z − (cid:18) E TeV (cid:19)(cid:32) B µ G (cid:33) − pc . (4)Solving equation (3) is computationally difficult due to thenumber of directional changes the particles make spiraling aroundthe magnetic field lines. Therefore, the approach that Harding et al.(2016) and many others take is to approximate the particle mo-tion as simply following the field line it experiences. Using thisapproximation, one step through the simulation is a step over whichthe amplitude and direction of the total magnetic field is constant.Effectively, the particle’s velocity vector is either parallel or antipar-allel to the magnetic field line, depending on its initial direction inapproaching it. The particle’s position can then be reset as: xt x v · BB ˆ Bt , (5)where x is the particle’s previous position, v is the initial directionof the particle’s velocity when approaching the field line, ˆ B is themagnetic field direction, and t is the time the particle took for thestep. The full path length traversed by the particle is ct .Equation (5) is assumed to be a valid approximation if theLarmor radius is on par with or smaller than the path length ofthe particle, under the assumption that solving equation (3) directlywould be unlikely to move the particle into a region with a differentmagnetic field. We note that it is possible that solving equation (3)explicitly can alter the particle motion, and should be studied infuture work.The time t that a particle is expected to take for one stepis dependant on the mean free path of the particle motion. Bydescribing the turbulent magnetic field interaction as a scattering term, we can alse describe the distance to scattering interaction. Thescattering mean free path describes how far the particle is expectedto travel before encountering a change in the turbulent componentof the magnetic field, and is determined by the particle energy andthe strength of the magnetic field. It is given by (Schlickeiser &Miller 1998; Fryer et al. 2007) λ sc × (cid:18) λ max cm (cid:19) (cid:18) B . mG (cid:19) − (cid:18) E TeV (cid:19) cm , (6)where λ max is the scale length of the turbulent magnetic field, B isthe amplitude of the total magnetic field (turbulent and coherent)and E is the proton energy. We use λ max cm, as was donein Harding et al. (2016). For the cosmic ray energies we consider,magnetic field scattering is the only significant particle interaction,and other interactions such as proton-proton scattering are negligi-ble. Therefore, we take the total mean free path to be equal to thescattering mean free path: λ λ sc . Particles with a mean free path λ follow an exponential probability distribution for their motion, Px e − x λ , where x is the distance traveled in one step. Solving thisdistribution for its cumulative distribution function (CDF) gives theexpression χ − e − x λ , where χ is a random variable sampled be-tween 0 and 1. This can be rearranged to obtain a distance traveledfor that step, along with the time it took for the particle to go thatdistance: t xc − ln χ (cid:18) λ c (cid:19) . (7) In regions with only a turbulent magnetic field, we can treat the trans-port in the diffusion approximation and we use a DDMC methodwhich is very similar to that used by Harding et al. (2016). In theseregions, particle motion changes direction too quickly for it to becomputationally feasible to track the particle motion directly. How-ever, particles in an isotropic magnetic field exhibit “random walk"motion, with their displacement proportional to the square root oftheir travel time. DDMC methods combine a number of smallerrandom walk steps that the particle would take into a larger stepbased on this principle. Rather than picking a distance to collisionusing equation (7), we instead pick a travel time and then sample aparticle distance traveled in that time.For a three dimensional random walk, the probability of mov-ing distance R after N steps is (Rycroft & Bazant 2005; Hardinget al. 2016) PR π R (cid:32) π Na (cid:33) exp (cid:32) − R Na (cid:33) , (8)where a is the expectation value for displacement in a single step. Toapproximate the transport motion, we relate the number of steps, N ,to the total distance traveled by the particle and the average step size: N ct λ . The expected distance traveled in one step a √︁ (cid:104) x (cid:105) − (cid:104) x (cid:105) √︁ (cid:104) x (cid:105) is calculated from the transport equation probability functionto be a √ λ . Putting in these values gives the expression: PR π R (cid:32) π ct λ (cid:33) exp (cid:32) − R ct λ (cid:33) . (9)This expression does not lend itself to sampling directly via aninversion technique, so we instead cast it in a hybrid form using bothinversion sampling and rejection sampling techniques. Specifically,we break the full probability distribution PR into two components, gR and hR . Defining the constant C ct λ , it can be seen that: PR ∝ Re − C R Re − C R gRhR . (10) MNRAS000
A particle transport Monte Carlo code tracks the particle motion asthe particle goes through a zone and passes into the next. Particleswithin a magnetic field of strength B propagate according to theequation of motion: dvdt cr L v × ˆ B , (3)where ˆ B is the direction of the net magnetic field, v is the directionof the particles velocity, and r L is the Larmor radius of the particlemotion. This causes particles to move in a circular path of radius r L around the magnetic field line, perpendicular to it. The componentof the particle’s initial velocity parallel to the magnetic field linedoes not change. Particles within a magnetic field of strength B ,energy E , and charge Ze have a Larmor radius of (Schlickeiser2002): r L . × − Z − (cid:18) E TeV (cid:19)(cid:32) B µ G (cid:33) − pc . (4)Solving equation (3) is computationally difficult due to thenumber of directional changes the particles make spiraling aroundthe magnetic field lines. Therefore, the approach that Harding et al.(2016) and many others take is to approximate the particle mo-tion as simply following the field line it experiences. Using thisapproximation, one step through the simulation is a step over whichthe amplitude and direction of the total magnetic field is constant.Effectively, the particle’s velocity vector is either parallel or antipar-allel to the magnetic field line, depending on its initial direction inapproaching it. The particle’s position can then be reset as: xt x v · BB ˆ Bt , (5)where x is the particle’s previous position, v is the initial directionof the particle’s velocity when approaching the field line, ˆ B is themagnetic field direction, and t is the time the particle took for thestep. The full path length traversed by the particle is ct .Equation (5) is assumed to be a valid approximation if theLarmor radius is on par with or smaller than the path length ofthe particle, under the assumption that solving equation (3) directlywould be unlikely to move the particle into a region with a differentmagnetic field. We note that it is possible that solving equation (3)explicitly can alter the particle motion, and should be studied infuture work.The time t that a particle is expected to take for one stepis dependant on the mean free path of the particle motion. Bydescribing the turbulent magnetic field interaction as a scattering term, we can alse describe the distance to scattering interaction. Thescattering mean free path describes how far the particle is expectedto travel before encountering a change in the turbulent componentof the magnetic field, and is determined by the particle energy andthe strength of the magnetic field. It is given by (Schlickeiser &Miller 1998; Fryer et al. 2007) λ sc × (cid:18) λ max cm (cid:19) (cid:18) B . mG (cid:19) − (cid:18) E TeV (cid:19) cm , (6)where λ max is the scale length of the turbulent magnetic field, B isthe amplitude of the total magnetic field (turbulent and coherent)and E is the proton energy. We use λ max cm, as was donein Harding et al. (2016). For the cosmic ray energies we consider,magnetic field scattering is the only significant particle interaction,and other interactions such as proton-proton scattering are negligi-ble. Therefore, we take the total mean free path to be equal to thescattering mean free path: λ λ sc . Particles with a mean free path λ follow an exponential probability distribution for their motion, Px e − x λ , where x is the distance traveled in one step. Solving thisdistribution for its cumulative distribution function (CDF) gives theexpression χ − e − x λ , where χ is a random variable sampled be-tween 0 and 1. This can be rearranged to obtain a distance traveledfor that step, along with the time it took for the particle to go thatdistance: t xc − ln χ (cid:18) λ c (cid:19) . (7) In regions with only a turbulent magnetic field, we can treat the trans-port in the diffusion approximation and we use a DDMC methodwhich is very similar to that used by Harding et al. (2016). In theseregions, particle motion changes direction too quickly for it to becomputationally feasible to track the particle motion directly. How-ever, particles in an isotropic magnetic field exhibit “random walk"motion, with their displacement proportional to the square root oftheir travel time. DDMC methods combine a number of smallerrandom walk steps that the particle would take into a larger stepbased on this principle. Rather than picking a distance to collisionusing equation (7), we instead pick a travel time and then sample aparticle distance traveled in that time.For a three dimensional random walk, the probability of mov-ing distance R after N steps is (Rycroft & Bazant 2005; Hardinget al. 2016) PR π R (cid:32) π Na (cid:33) exp (cid:32) − R Na (cid:33) , (8)where a is the expectation value for displacement in a single step. Toapproximate the transport motion, we relate the number of steps, N ,to the total distance traveled by the particle and the average step size: N ct λ . The expected distance traveled in one step a √︁ (cid:104) x (cid:105) − (cid:104) x (cid:105) √︁ (cid:104) x (cid:105) is calculated from the transport equation probability functionto be a √ λ . Putting in these values gives the expression: PR π R (cid:32) π ct λ (cid:33) exp (cid:32) − R ct λ (cid:33) . (9)This expression does not lend itself to sampling directly via aninversion technique, so we instead cast it in a hybrid form using bothinversion sampling and rejection sampling techniques. Specifically,we break the full probability distribution PR into two components, gR and hR . Defining the constant C ct λ , it can be seen that: PR ∝ Re − C R Re − C R gRhR . (10) MNRAS000 , 1–14 (2020) onte Carlo CR Transport Figure 2.
Shows the normalized probability density functions using C=1 forthe full 3D random walk distribution PR (Equation 9), the 2D approximationwe use gR (Equation 11), and the 3D correction hR . The red points show theaverage value of R sampled for PR and gR . The function g(R) is the 2D random walk probability distribution,while h(R) has the correction for the full 3D solution (Rycroft &Bazant 2005).To sample the complete function, we first sample a value R samp via inversion of the CDF for gR , and then choose to accept R samp based on a rejection sample of the function hR . The efficiency ofthe rejection step is 80 % , and does not dramatically increase thecomputational cost of each run. However, for algorithm simplicity,most of our runs employ an approximate sampling technique wherethe rejection step is removed. To justify this, we compared selecttest cases using the complete function and discovered that, withinour error range, there was no noticeable difference in results whenthe rejection step was removed. This is largely because gR is aclose approximation to PR in terms of mean diffusion distanceproperties. The primary effect of removing the rejection check wasto slightly decrease the diffusive component’s mean free path. Thisresults in a slight enhancement of the influence of the diffusive fieldcomponent relative to the coherent field component. As such, usingthe 2D approximation to the full 3D solution provides a conservativeestimate of the influence of coherent magnetic fields on the cosmicray flux. These results can be seen intuitively from Figure 2, whichshows the three functions PR , gR , and hR for C=1. Plotted over PR and gR in the red points is the average value of R sampled for thetwo distributions. The function gR has a lower peak and is skewedfarther towards lower distances than PR , leading to a slightly loweraverage value of R sampled.With this approximation, the employed probability distributionfunction is given by: PR Re − C R , (11)which yields a sampled distance given by: R √︁ − ln χ ct λ. (12)The values we use for t are t years and t years, dependingon the particle’s situation. We discuss this more in Section 4.If the particle’s travel takes it to the edge of a zone wall beforethe end of its step, it stops and checks whether it is going into theregion with a coherent magnetic field. If it is, it goes to the edge of the zone and exits the diffusion loop, with a travel time of xc .Otherwise, it continues its journey along the path until it either hitsanother zone wall or finishes its timestep. Therefore all particles willkeep going until they either reach the end of the sampled distanceor hit the edge of a zone in which the transport mechanism isnecessary. For each substep along its sampled distance that it hitsa zone wall, the particle’s remaining time along its current path isrecalculated as the difference between its originally sampled timeand its substep time along that path. Once the particle reaches theend of its timestep, its direction is resampled.Note that, although equations (1) and (2) are used in both atransport method and DDMC method to approximate the turbulentmagnetic field, they have a different meaning for the two. Becausethe transport algorithm approximates the particle motion directly,the randomly sampled direction represents the direction of the tur-bulent magnetic field influence for that particle step. In the diffusionalgorithm, this vector represents a direction of particle travel overmany steps, through many different magnetic field lines. For the region inside the coherent magnetic field, there is no clear an-alytic solution to describe the particle motion, arising from the factthat we can neither assume the diffusion approximation or assumethe particle flows along field lines. Instead, the solution lies some-where between these two extremes. For the particle transport MonteCarlo method, the magnetic field direction is sampled, and the nextdirection of particle motion is based on the particle’s initial velocitydirection with respect to this vector. In contrast, for the diffusionapproximation, the particles direction change after many steps issampled; hence the particle’s previous direction of travel doesn’tmatter. Therefore, for the region containing the coherent magneticfield, we explored two methods for modeling the particle motionwhich attempt to “combine" the transport and diffusive methods indifferent ways. Both of these methods give solutions for particletransport that are somewhere between the isotropic magnetic fieldsolution and the transport solution.The first method we explored (“Method 1") forces particlesto follow either transport motion or diffusive motion inside thecoherent magnetic field box. For every step, we sample a randomnumber r between 0 and 1 and compare that to the coherent fieldamplitude ratio ( g B g B g B t ). If r < g , then the particle uses transportmotion and follows the coherent magnetic field with one of twodirections, based on the velocity of its previous step. If r > g , thenthe travel direction is sampled randomly. Thus, for example, if theamplitudes of the coherent and turbulent field are the same ( B g B t . ), then g . , and a particle that travels inside the box will followtransport motion about half of the time and random walk motionthe other half.The second method we explored (“Method 2") uses a directionof travel which is the vector sum of the coherent magnetic fieldvector and the turbulent magnetic field vector. The turbulent mag-netic field vector is randomly sampled. The coherent field vector,however, in this case, represents the direction the field line wouldgive the particle rather than the field line direction itself. It is one oftwo directions; either parallel to the field or antiparallel to it, basedon the particle’s previous direction of travel.A consequence of using a timestep t years outside thecoherent magnetic field box is that it was inadequate to simply usethe transport solution inside the coherent field box and the diffusionapproximation outside. The boundary conditions of the transportsolution require that when particles leave the transport box they do MNRAS , 1–14 (2020)
Fitz Axen et al. not resample their direction; which means they are biased to moveaway from the box. The distance which particles have to travel torecover random walk motion can be relatively large in our case,since with a timestep t years the average distance a particletravels in one diffusion step is Rt . pc. This causes a deficit ofparticles inside the coherent magnetic field region.In order to rectify this situation, we made two modificationsto stop the particles that entered the coherent field region fromquickly propagating away from it. First, we put a “sheath" of onegrid zone around the coherent magnetic field in which the particlesfollowed the transport solution, if they had just exited the fieldregion. Second, once the particles had exited the sheath region intothe diffusion regime, they used a timestep of t years rather than t years. These modifications allowed us to get the correctprecision for particles that entered the field region while runningthe particles that never entered the field region quickly. One testable property of an isotropic magnetic field is that the parti-cle motion follows a certain distance distribution. Because the meanfree path of the particle motion is the same everywhere, the particlesshould always follow a root-mean square (rms) distance distributionfor their travel time following equation (12) of Rt √ ct λ . To con-firm our code reproduces this property, we extended the simulationbox to calculate a broad particle distribution.In order to make this simulation run faster with a larger trans-port box, we set the mean free path in every grid zone to λ . pc,which is approximately the mean free path of a 100 TeV- 1 PeVproton in a 3 µ G magnetic field (our run simulations, with 10 TeVprotons, have a mean free path of approximately λ . pc. Wespecified a “box region" given by 42 pc < y,z <
102 pc and 90 pc < x <
150 pc, in which the packets must follow the transport solution;just as would be done if there was a coherent magnetic field in thisregion. With this setup, about 65 percent of the total number ofparticles run entered this region at some point during their lifetime,and they spent on average 12 percent of their lifetime in the boxregion.In total, we did 100 runs of the code, each with 10000 particles.Each particle’s position was recorded every 20000 years. Figure 3shows the difference between the rms distance traversed for the10000 particles of a run at each time and the analytic solution atthat time. Each point on the plot represents a separate run. The linerepresents the average of these 100 points at each time. This testwas done for both just DDMC and the IMC/DDMC hybrid.As can be seen, the points for the IMC/DDMC hybrid methoddoes not reproduce the exact solution (averaging above 0). However,the average is below 1%, and does not show a strongly increasingtrend with time. The distribution of the points never exceeds 4 per-cent. For this paper, this error lies below the numerical uncertaintyin our Monte Carlo statistics (Section 4.2) and the errors from ourmethods are always less than these statistical errors.
The goal of this paper is to study the effect of a coherent magneticfield on the particle flux at the tally plane. Although it alters theflux across this entire boundary, the primary affect is at the positionof the coherent magnetic field box (in the space 108 pc < y,z <
120 pc). With a finite set of particles, statistical uncertainties oftendominate the errors in a Monte Carlo approach. To test this, we firstran the particles through a purely isotropic magnetic field, usingthe diffusion algorithm alone. We then defined the particle flux F for a magnetic field model i as being the ratio between the numberof particles that hit this region N hit, i to the number that hit thisregion in the isotropic run N hit, back , relative to the total numbers ofparticles N tot, i and N tot, back run for both: F N hit, i N tot, i N hit, back N tot, back N tot, back N tot, i N hit, i N hit, back , (13)Statistical errors in Monte Carlo methods cause different runs to giveslightly different numbers for this quantity. However, a set of runswith the same magnetic field configuration should form roughly aGaussian distribution for the computed value of F , with 95 % of thevalues within σ of the mean. Because for every box configurationwe only wanted to do one run, we had to relate the σ value (cid:15) F ofthis distribution to an uncertainty measure in a single run based onthe number of particle counts.For counting statistics, the uncertainty in a counted number ofparticles is equal to the square root of the number of particles forthe statistic. We can therefore compute upper and lower uncertaintybounds for the computed flux above the background using countsfor the number of particles that exit the box region and counts thatexit the box region for the background run: F up N tot, back N tot, i N hit, i √︁ N hit, i N hit, back − √︁ N hit, back , (14) F low N tot, back N tot, i N hit, i − √︁ N hit, i N hit, back √︁ N hit, back . (15)The difference between these two quantities ∆ F gives another mea-sure for the σ uncertainty interval around the computed value forthe flux above the background: F ± σ F ± ∆ F F ± F up − F low . Inorder to see whether ∆ F computed for one run was approximatelyequal to what the value of (cid:15) F would be for the distribution, we setup a test configuration used in the full set of runs; a box extent of24 pc, offset of 0 pc, and B g B t . ; and did 100 separate particleruns for this same configuration. We did this for both Method 1 andMethod 2.For both configurations, we computed the mean flux µ and stan-dard deviation σ in the flux for the set of runs. We found that usingthese quantities, the data accurately fit a Gaussian distribution. Wealso confirmed that the average computed value of ∆ F for the setvery closely matched the σ value (cid:15) F of the Gaussian distribution,and consequently the average computed value of ∆ F F for the setwas close to the value of σµ from the Gaussian distribution; in fact,for the tests, they were only about 10 % apart from each other. Forour production runs we wanted to be within 10 % of the "correct"value of F, with 95 % certainty; so, for every box configuration, weran enough particles such that the σ value (cid:15) F F ∆ FF was lessthan 0.1.Figure 4 shows histograms of the computed flux F above thebackground for the configuration tested for Method 1 and Method2 along with a Gaussian fitted to the mean and standard deviationsfor each. Table 2 shows the results of the data from the Gaussiandistribution. MNRAS000
120 pc). With a finite set of particles, statistical uncertainties oftendominate the errors in a Monte Carlo approach. To test this, we firstran the particles through a purely isotropic magnetic field, usingthe diffusion algorithm alone. We then defined the particle flux F for a magnetic field model i as being the ratio between the numberof particles that hit this region N hit, i to the number that hit thisregion in the isotropic run N hit, back , relative to the total numbers ofparticles N tot, i and N tot, back run for both: F N hit, i N tot, i N hit, back N tot, back N tot, back N tot, i N hit, i N hit, back , (13)Statistical errors in Monte Carlo methods cause different runs to giveslightly different numbers for this quantity. However, a set of runswith the same magnetic field configuration should form roughly aGaussian distribution for the computed value of F , with 95 % of thevalues within σ of the mean. Because for every box configurationwe only wanted to do one run, we had to relate the σ value (cid:15) F ofthis distribution to an uncertainty measure in a single run based onthe number of particle counts.For counting statistics, the uncertainty in a counted number ofparticles is equal to the square root of the number of particles forthe statistic. We can therefore compute upper and lower uncertaintybounds for the computed flux above the background using countsfor the number of particles that exit the box region and counts thatexit the box region for the background run: F up N tot, back N tot, i N hit, i √︁ N hit, i N hit, back − √︁ N hit, back , (14) F low N tot, back N tot, i N hit, i − √︁ N hit, i N hit, back √︁ N hit, back . (15)The difference between these two quantities ∆ F gives another mea-sure for the σ uncertainty interval around the computed value forthe flux above the background: F ± σ F ± ∆ F F ± F up − F low . Inorder to see whether ∆ F computed for one run was approximatelyequal to what the value of (cid:15) F would be for the distribution, we setup a test configuration used in the full set of runs; a box extent of24 pc, offset of 0 pc, and B g B t . ; and did 100 separate particleruns for this same configuration. We did this for both Method 1 andMethod 2.For both configurations, we computed the mean flux µ and stan-dard deviation σ in the flux for the set of runs. We found that usingthese quantities, the data accurately fit a Gaussian distribution. Wealso confirmed that the average computed value of ∆ F for the setvery closely matched the σ value (cid:15) F of the Gaussian distribution,and consequently the average computed value of ∆ F F for the setwas close to the value of σµ from the Gaussian distribution; in fact,for the tests, they were only about 10 % apart from each other. Forour production runs we wanted to be within 10 % of the "correct"value of F, with 95 % certainty; so, for every box configuration, weran enough particles such that the σ value (cid:15) F F ∆ FF was lessthan 0.1.Figure 4 shows histograms of the computed flux F above thebackground for the configuration tested for Method 1 and Method2 along with a Gaussian fitted to the mean and standard deviationsfor each. Table 2 shows the results of the data from the Gaussiandistribution. MNRAS000 , 1–14 (2020) onte Carlo CR Transport Figure 3.
Shows, for the isotropic magnetic field verification, the rms distance errors from the analytic solution calculated at each time of the 100 runs done.The plot on the left uses just the DDMC algorithm, while the plot on the right uses the hybrid DDMC/transport algorithm. The line shows the average of thesepoints at each time. The average error never reaches 1 percent and tends to level out around 80000 years. The distribution of points at every time never spansmore than 4 percent.
Figure 4.
Distributions for the uncertainty test described in the text, showing histogram of flux uncertainties obtained for 100 runs of the two methods. Thevertical red line in each plot marks the average computed mean of the data and the other tick markers show the 1, 2, and 3 σ intervals around the mean. Note:The obtained fluxes are different for the cases due to the different methods used, which is discussed more in Section 5. Table 2.
Describes the data found by the uncertainty test described in the text.The second column shows the mean value of ∆ F F for the set of 100 runs foreach configuration (computed for each run individually). The next columnshows the ratio between the mean and standard deviation for the Gaussiandistribution of all the flux values computed. The fourth column shows thepercentage difference between these quantities, and the final column showspercentages within range of the mean for the Gaussian distribution of theflux values.M (cid:104) ∆ F F (cid:105) σµ (cid:104) ∆ F F (cid:105)− σµσµ , , σ With this code, we can calculate the range of effects varying thecoherent magnetic field properties outlined in Section 2.1. Figure5 shows two contour plots of the particle flux a coherent magneticfield produces above the background at the tally plane. As can beseen, the coherent magnetic field produces a noticeable flux abovethe background, shown by the yellow and green patches in the upperright corners of the plots. These patches span 4x4 grid zones in ourcase, and the flux F previously discussed contains the relative sumof all of the extra particles in this region.In these plots, the likelihood of the particles to follow the di-rection of the coherent field can clearly be observed. One noticeablefeature of the plots is that the effect on the tally plane is not alwayssimply a higher flux (darker patch) on the contour plot. Rather, themagnetic field configurations with a stronger effect tend to forma “ring" of increased flux due to the likelihood of the particles to MNRAS , 1–14 (2020)
Fitz Axen et al. follow the coherent field when they encounter it in the outer zonesof the coherent field region. Another feature of the stronger mag-netic field configurations is the decreased flux surrounding the boxand extending up into the right corner of the plot. This shadowingeffect was also observed in Harding et al. (2016), and is due to alack of particles scattering past the coherent field without enteringit and getting “trapped". Because our coherent field box contains aturbulent component, there is some expectation that the shadowingeffect would be muted, but it is still present in many cases.The computed flux results of our code are shown in Figure 6.The plot on the left shows the results using “Method 1" while the ploton the right shows the results using “Method 2". The figure shows thecomputed flux above the background F for every box extent, tallyplane offset, and coherent magnetic field ratio tested. Each colorshows a different coherent magnetic field box extent, with each colorhaving three different lines representing the three different coherentmagnetic field ratios. The dotted lines are put in for comparison andwere computed using the methods of Harding et al. (2016); i.e., noturbulent magnetic field in the coherent field box. In general, thetopmost line of one color is the highest coherent magnetic field andthe coherent field decreases moving down in height; this holds upexcept for the lower points close to the background of F , wherestatistical uncertainty causes more variation.These results demonstrate the large effect of a coherent mag-netic field component on the particle flux above the background.Indeed, one can see from Figure 6 that for all runs done in the studyof coherent magnetic field amplitude, coherent field box length, andcoherent field box offset from the tally plane, there was a noticeableflux observed above the background level; and above the observedanisotropy level. Only a few of the runs with the lowest flux valuesdid not have a σ range above F , which signifies a 95 % chanceof a noticeable observed flux. These points are shown by crossesinstead of dots in Figure 6. Both of the hybrid transport methods, Method 1 and Method 2,produce the same, generally predictable trends for the variationof the tested box configuration variables. Increasing the coherentmagnetic field ratio, increasing the box length, and decreasing itsdistance from the tally plane all increase the observed flux at theplane. However, the two methods do produce quantitatively-differentresults from each other, depending on the coherent magnetic fieldratio used.In addition to our main configuration parameter space, we didone study where we picked a single box offset of 0 pc and extentvalue of 6 pc and varied the value of the coherent magnetic fieldratio for 10 different values, ranging from 0.1 to 1.0. This is shown inFigure 7. It can be seen that, for higher magnetic field ratios ( B g B t (cid:38) . ), Method 2 predicts a higher particle flux above the backgroundthan Method 1 while, for lower magnetic field ratios ( B g B t (cid:46) . ),Method 1 predicts a higher particle flux above the background thanMethod 2. This means that over the span of magnetic field ratiovalues we tested for the main study, Method 2 has a much broaderrange in predicted flux values (because higher coherent magneticfield ratios predict a higher particle flux above the background ingeneral).These differences can be seen in the dependence of the resultson the spatial distribution of the coherent magnetic field. Figure 8shows the particle flux at tally plane as a function of the distanceaway from the center of the box region (108 pc < y,z <
120 pc).The top two plots were computed using Method 1, and the bottom two plots were computed using Method 2. The left column uses acoherent magnetic field ratio of B g B t . while the right columnuses a coherent magnetic field ratio of B g B t . . As can be seen,for the higher magnetic field ratio, Method 2 produces a higherparticle flux around the box region, and a greater deficit outside.However, for the lower field value, Method 1 still has an observableflux around the box while it is much harder to see for Method 2.The differences between these two methods provide an indi-cation of the uncertainty in the predicted particle anisotropy at thetally plane. For our tests, Method 1 predicts an anisotropy for allbox geometries with B g B t . and B g B t . , and for most of thegeometries with B g B t . . Method 2, on the other hand, only pre-dicts an anisotropy for all box configurations with B g B t . ; all ofthe configurations with B g B t . and many with B g B t . are notpredicted to be high enough to be seen within our 10% simulationerrors. Finally, we tested one magnetic field configuration using higherenergy particles, for both Method 1 and Method 2. We tested 100TeV, 1 PeV, 10 PeV, and 100 PeV particles. The magnetic fieldconfiguration we chose to test was a box extent of 24 pc, offset of 0pc, and coherent magnetic field ratio of B g B t . . We found that thetrend is that higher energy particles produce a smaller flux abovethe background in the box. This is shown in the Figure 9. The reasonfor this likely has to do with the higher energy particles having alonger path length, which means they are more likely to "escape"from the coherent magnetic field. We now attempt to approximate the effect that a coherent magneticfield might have on the observed cosmic ray flux at Earth. Thoughwe have demonstrated theoretically that a coherent magnetic fieldcan create anisotropies in the cosmic ray flux across a tally plane,this doesn’t correlate exactly to the anisotropies observed in thecosmic ray arrival direction at Earth. The observed flux at a pointrequires tallying the angular distribution of the flux as well. Forstatistical results, section 5 tallied the entire flux. In this section, wealso study the angular distribution.To understand the angular distribution of the material on ourtally plan, we ran one configuration where we also tallied the velocity u of the particles as they hit the plane. We ran particles using theconfiguration tested for Method 1 with a box extent=24 pc, offset=3pc, and B g B t . (though any of the runs with an offset greater than0 pc could have been used). For every particle that hit the plane, wecalculated cos θ u ˆ u · ˆ n , where ˆ n , , and θ u is the angle between itsvelocity and the plane surface normal. We then tallied the numberof particles between cos θ u and cos θ u ∆ cos θ u for values of cos θ u between 0 and 1. We used 15 bins in cos θ u , so ∆ cos θ u .We found that the velocity distribution of the particles beingemitted from the plane follows the relationship F ˆ u ∝ ˆ u · ˆ n cos θ u . (16)The results of this test are shown in Figure 10. Each point representsa total number of particles emitted in one bin. In addition to showing This is a Lambertian distribution in angle, which is expected for particlesfollowing isotropic scattering motion. MNRAS , 1–14 (2020) onte Carlo CR Transport Figure 5.
Shows two examples of the particle flux above the background F for our simulation runs. Both were done using the same coherent magnetic fieldconfiguration of B g B t . , box extent=24 pc, and box offset=0 pc, but the left plot was computed using Method 1 while the right plot was computed usingMethod 2. The coherent magnetic field region can be seen by the yellow and green patches in the upper right corner of the plots. Figure 6.
Shows the results of our code, with the particle flux above the background defined in Section 4.2 as a function of the coherent field offset fromthe tally plane. Each color represents a different box extent, which was tested for three different values of the magnetic field. The dotted lines shows valuescomputed using the methods of Harding et al. (2016). Every run is done to 10 % uncertainty. The points show the computed values, with the dots being pointswhose uncertainty range is above the background while the stars have an uncertainty range below the background. the calculation done for the entire tally plane, we also show thesame calculation for two regions on the plane; one that is the boxregion used to compute the flux, and another region that is in theopposite corner of the plane. All lines are normalized to the numberof particles that were considered for that region. These details areto emphasize the fact that although these regions have differenttotal numbers of particles and particle density, they all follow thesame distribution in momentum. The main difference between thedifferent regions is the offset of the points from the best fit line;the opposite corner region has a higher offset than the box regionor entire plane because there are less particles being considered,leading to higher statistical error.We can use the angular distribution of the particles at the tally surface inferred from this study to calculate the observed flux inthe observer frame (as a function of viewing angle). To do so, wedetermine an observer location ("position of the Earth") at coor-dinates x e , y e , z e . We calculate each particle’s position in sphericalcoordinates [ φ,θ ] from its Cartesian coordinates x p , y p , z p on thetally plane using tan φ y p − y e x p − x e and tan θ z p − z e x p − x e . Eachparticle is then in a spherical grid zone Ω z , where each zone Ω z bounds a two dimensional region in [ φ , θ ] space. By summing onlythese angle-dependent contributions from the tally surface, we cal-culate the observed sky distribution for part of the sky that can bemapped from the tally plane. For each spherical grid zone Ω z , wetallied the proportional flux of particles that would arrive at the MNRAS , 1–14 (2020) Fitz Axen et al.
Figure 7.
Shows the results of varying the magnetic field ratio for one boxgeometry configuration (extent=6 pc, offset=0 pc).
Earth point from that zone as R Ω z ∝ N Ω z i cos θ e , (17)where the sum is over all of the N Ω z particles that hit the plane inthe zone Ω z and θ e is the angle between the surface normal andthe vector connecting the particle hit position to the Earth point.For the plots we present, the Earth point is at position x e , y e , z e , , pc, and the angular resolution in both dimensions are ∆ θ π and ∆ φ π (50 angular bins each), with θ ∈ ,π and φ ∈ ,π .To demonstrate the cosmic ray flux in a way that is commonly doneobservationally, and reduce the effect of the higher noise floor dueto the geometric effect of the angular bins, for each solid anglebin we chose to consider the difference between the flux and thebackground run flux rather than the ratio: F Ω z R Ω z − R back Ω z ∝ N Ω z i cos θ e − N back, Ω z i cos θ e , (18)where here all calculations were done using the same number of par-ticles originally started for both the main and the background runs( N tot, i N tot, back ), since the number of particles will not normalizeout in these plots.Figure 11 illustrates this process for one box configuration us-ing Method 2. The top two panels show R Ω z and R back Ω z computedusing Equation 17. The higher particle counts in the bottom rightcorners of these panels appear because of the geometric location ofthe earth sky point compared to the tally plane and the location ofthe cosmic ray source. The placement of the earth sky point nearthe top left corner of the tally plane causes more zones on the planeto be included in the lower right angular bins, and fewer in thetop left angular bins. Additionally, there are more particle countsin the tally plane zones near the cosmic ray source. In the bottompanel, computed using Equation 18, these effects are removed andthe location of the coherent magnetic field can clearly be seen bythe higher particle flux.Finally, we calculate for each plot a generic noise level σ Ω ,where σ Ω | min F Ω z | . This simply means that the noise level for allof the angular bins for one plot is the magnitude of the most negativevalue computed, or the flux in the angular bin where the isotropicrun had the highest flux above the main run. Though the uncertaintyis technically a function of angular zone, this method is a simple way to calculate the uncertainty, and is an overestimate for most ofthe grid. Plotting F Ω z σ Ω allows us to neglect the proportionalityconstant in these expressions, as it should be the same for F Ω z and σ Ω . The results are presented in Figures 12 and 13, which showcontour plots of F Ω z σ Ω at the earth point for Methods 1 and 2.Each row in the figure shows the effect of varying one ofthe magnetic field configuration parameters; the box extent, coher-ent magnetic field ratio, and earth sky offset. In these plots, thecolored region in the upper left corner shows the location of the co-herent magnetic field box. The effect of varying the three magneticfield configuration parameters produced what might be the expectedchanges on the background flux at the Earth sky. Increasing the boxextent, decreasing its distance from the sky, and increasing the co-herent magnetic field ratio produce a higher particle flux observedon the sky. This is in contrast to Harding et al. (2016), where allcoherent magnetic field values produce the same result as there is noturbulent field in the box to drive the particles out.We have variedthe observer position and the magnitude of the anisotropy varies,but not significantly.The anisotropies presented here are much higher than the ob-served anisotropy in the cosmic ray flux, implying that the strengthof the coherent magnetic fields can be much lower than what weassumed. Even so, we are assuming the magnetic field energy inthe coherent fields is much less than that of the small-scale struc-tures. The strength of our coherent magnetic fields was not chosento match the observed anisotropies, but to demonstrate the role suchmagnetic fields can play on the flux observed at the earth. Statisticallimitations of our Monte Carlo method required higher coherentmagnetic field strengths than those needed to explain the observa-tions. Finally, we note that limitations with our grid setup, such asparticles not being allowed to scatter back across the tally plane,may also have a minor effect on our results. In this study , we have extended the work of Harding et al. (2016) toinclude more realistic magnetic field structures which included bothan isotropic and an anisotropic component. Additionally, we variedour magnetic field structures in a different way than was previouslystudied, varying the length of the field structure and moving it awayfrom the tally plane. To accomplish this hybrid solution, we proposetwo methods which "combine" the standard transport and diffusionalgorithms usually used. One of these methods involved the particlesswitching off between the two methods, picking the method of travelbased on a probabilistic approach. The other method involved theaddition of two vectors intended to represent the two componentsof the magnetic field separately.We have shown that the anisotropies in the particle flux atthe tally plane are much easier to create than might be expected, nomatter which method is used. In particular, with a strong enough co-herent magnetic field component, there are noticeable anisotropiesabove the background for all coherent field length scales and offsetvalues tested. Because the length scale of the coherent field is somuch larger than the mean free path of particle motion, even a weakcoherent magnetic field component can give the particle a strongenough directional push over time to move it into the box regionat the tally plane. Additionally, even magnetic field structures re-moved from the plane can cause a noticeable particle flux abovethe background. This is because the coherent magnetic field altersthe particle transport, effectively creating what appears to be a newcosmic ray source.
MNRAS000
MNRAS000 , 1–14 (2020) onte Carlo CR Transport Figure 8.
Shows the number of particles at the Earth sky as a function of box offset and radial distance from the center of the box region. The top plots werecomputed using Method 1 and the bottom plots were computed using Method 2. The left column uses a magnetic field ratio of B g B t . while the right columnuses a magnetic field ratio of B g B t . . Both have a box extent of 24 pc. Figure 9.
Shows the results of varying the energy for one box geometry andmagnetic field ratio (extent=6 pc, offset=0 pc, B g B t . ) Finally, we took the results at the tally plane and used themto make a construction of what an observer at Earth might see.Though this method had many limitations, we showed that thestronger, closer magnetic field configurations produce a noticeableanisotropy. Many other works have found a similar effect to the re-
Figure 10.
Shows the number of particles emitted in certain differentdirections as a function of cos( θ ), run for one box configuration. We showthe results for the entire x=150 pc plane (red), as well as the box regionalone (blue) and one other region at a different location in the grid (green).Each line is normalized to the number of particles being considered in thatregion. The black lines are a linear fit to the red points.MNRAS , 1–14 (2020) Fitz Axen et al.
Figure 11.
Plots of R Ω z (top left), R back Ω z (top right), and F Ω z (bottom) using Method 1 and B g B t =1.0, box extent=24 pc and box offset = 3 pc. sults presented here and in Harding et al. (2016). In several studies(Giacinti & Sigl (2012),Ahlers & Mertsch (2015), Ahlers (2014),López-Barquero et al. (2016)), it was shown that anisotropies in theflux can arise from local turbulent magnetic field structures. Oursimulations were only sensitive to anisotropies that are much largerthan those observed. The fact that modest coherent magnetic fieldscan produce strong anisotropies show that the features needed toexplain the observed fields can be quite small and the small ampli-tudes of the observed cosmic ray anisotropies place limits on thenature of these coherent magnetic fields.We stress that although we have made improvements over pre-vious studies modeling cosmic ray transport, we have not done the“ultimate" treatment of the fields. This involves resolving the tur-bulent magnetic field structure and using this in combination withthe coherent field, directly solving the Lorentz equation of motionfor every particle step. Although such full treatments are beyondcomputational power for large grids, it is possible to use small-scalecalculations to improve on the recipes (Methods 1 and 2) used here.We defer this more detailed study for a later paper. Another way to reduce the uncertainties in this study wouldbe to obtain better measurements of the magnetic field structuresin the solar neighborhood. With these structures, the cosmic rayanisotropies may be better understood. DATA AVAILABILITY
The data underlying this article are available in the article.
ACKNOWLEDGEMENTS
This work was supported by the US Department of Energy throughthe Los Alamos National Laboratory. Los Alamos National Labo-ratory is operated by Triad National Security, LLC, for the NationalNuclear Security Administration of U.S. Department of Energy(Contract No. 89233218CNA000001) This research was supportedin part by the National Science Foundation under Grant No. NSFPHY-1748958 and by NASA ATP grant 80NSSC20K0507.
MNRAS000
MNRAS000 , 1–14 (2020) onte Carlo CR Transport Figure 12.
Shows plots of F Ω z σ Ω for different box configurations using Method 1. The different rows show the effects on the Earth sky image of varying thebox extent, coherent magnetic field ratio, and offset. The top row uses box extent values of 6 pc and 24 pc while keeping B g B t . and offset=3 pc. The middlerow uses values of B g B t . , 0.5, and 1.0 while keeping the box extent=24 pc and box offset=3 pc. The bottom row uses box offset values of 15 pc, 9 pc, and 3pc while keeping the box extent=24 pc and B g B t =1.0. REFERENCES
Aartsen M. G., et al., 2013, ApJ, 765, 55Abbasi R., et al., 2010, ApJ, 718, L194Abbasi R., et al., 2011, ApJ, 740, 16Abdo A. A., et al., 2008, Physical Review Letters, 101, 221101Abdo A. A., et al., 2009, ApJ, 698, 2121Aglietta M., et al., 2009, ApJ, 692, L130Ahlers M., 2014, Physical Review Letters, 112, 021101Ahlers M., Mertsch P., 2015, ApJ, 815, L2Amenomori M., et al., 2005, ApJ, 626, L29Amenomori M., et al., 2006, Science, 314, 439Amenomori M., et al., 2010, ApJ, 711, 119Biermann P. L., Becker Tjus J., Seo E.-S., Mandelartz M., 2013, ApJ, 768,124Blasi P., Amato E., 2012, J. Cosmology Astropart. Phys., 1, 011Cui S., 2011, International Cosmic Ray Conference, 1, 6Desiati P., Lazarian A., 2013, ApJ, 762, 44Drury L. O. ., Aharonian F. A., 2008, Astroparticle Physics, 29, 420Erlykin A. D., Wolfendale A. W., 2006, Astroparticle Physics, 25, 183Evans T. M., Urbatsch T. J., Lichtenstein H., Morel J. E., 2003, Journal of Computational Physics, 189, 539Fryer C. L., Liu S., Rockefeller G., Hungerford A., Belanger G., 2007, ApJ,659, 389Giacinti G., Sigl G., 2012, Physical Review Letters, 109, 071101Guillian G., et al., 2007, Phys. Rev. D, 75, 062003Harding J. P., 2013, preprint, ( arXiv:1307.6537 )Harding J. P., Fryer C. L., Mendel S., 2016, ApJ, 822, 102Kim W.-T., Kim C.-G., Ostriker E. C., 2020, ApJ, 898, 35Kotera K., Perez-Garcia M. A., Silk J., 2013, Physics Letters B, 725, 196Krasheninnikov S. I., 2013, Journal of Plasma Physics, 79, 1011Kritsuk A. G., Ustyugov S. D., Norman M. L., 2017, New Journal of Physics,19, 065003Kulsrud R. M., Zweibel E. G., 2008, Reports on Progress in Physics, 71,046901Lazarian A., Desiati P., 2010, ApJ, 722, 188López-Barquero V., Farber R., Xu S., Desiati P., Lazarian A., 2016, ApJ,830, 19Malkov M. A., Diamond P. H., O’C. Drury L., Sagdeev R. Z., 2010, ApJ,721, 750Munakata K., Mizoguchi Y., Kato C., Yasue S., Mori S., Takita M., Kóta J.,2010, ApJ, 712, 1100MNRAS , 1–14 (2020) Fitz Axen et al.
Figure 13.
Shows plots of F Ω z σ Ω for different box configurations using Method 2. The different rows show the effects on the Earth sky image of varying thebox extent, coherent magnetic field ratio, and offset. The top row uses box extent values of 6 pc and 24 pc while keeping B g B t . and offset=0 pc. The middlerow uses values of B g B t . , 0.5, and 1.0 while keeping the box extent=24 pc and box offset=0 pc. The bottom row uses box offset values of 2 pc, 1 pc, and 0pc while keeping the box extent=24 pc and B g B t . .Pohl M., Eichler D., 2013, ApJ, 766, 4Rycroft C., Bazant M., 2005, Technical report, Introduction to RandomWalks and Diffusion. Department of Mathematics, Massachusetts In-stitue of Technology.Salvati M., Sacco B., 2008, A&A, 485, 527Saveliev A., Jedamzik K., Sigl G., 2012, Phys. Rev. D, 86, 103010Schlickeiser R., 2002, Cosmic Ray AstrophysicsSchlickeiser R., Miller J. A., 1998, ApJ, 492, 352Schwadron N. A., et al., 2014, Science, 343, 988Strong A. W., Moskalenko I. V., Ptuskin V. S., 2007, Annual Review ofNuclear and Particle Science, 57, 285Sveshnikova L. G., Strelnikova O. N., Ptuskin V. S., 2013, AstroparticlePhysics, 50, 33This paper has been typeset from a TEX/L A TEX file prepared by the author. MNRAS000