Magnetohydrodynamic-Particle-in-Cell Method for Coupling Cosmic Rays with a Thermal Plasma: Application to Non-relativistic Shocks
Xue-Ning Bai, Damiano Caprioli, Lorenzo Sironi, Anatoly Spitkovsky
DDraft version March 31, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
MAGNETOHYDRODYNAMIC-PARTICLE-IN-CELL METHOD FOR COUPLING COSMIC RAYS WITH ATHERMAL PLASMA: APPLICATION TO NON-RELATIVISTIC SHOCKS
Xue-Ning Bai , Damiano Caprioli , Lorenzo Sironi , Anatoly Spitkovsky Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, 60 Garden St., MS-51, Cambridge, MA 02138and Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544
Draft version March 31, 2018
ABSTRACTWe formulate a magnetohydrodynamic-particle-in-cell (MHD-PIC) method for describing the inter-action between collisionless cosmic ray (CR) particles and a thermal plasma. The thermal plasma istreated as a fluid, obeying equations of ideal MHD, while CRs are treated as relativistic Lagrangianparticles subject to the Lorentz force. Backreaction from CRs to the gas is included in the form ofmomentum and energy feedback. In addition, we include the electromagnetic feedback due to CR-induced Hall effect that becomes important when the electron-ion drift velocity of the backgroundplasma induced by CRs approaches the Alfv´en velocity. Our method is applicable on scales muchlarger than the ion inertial length, bypassing the microscopic scales that must be resolved in con-ventional PIC methods, while retaining the full kinetic nature of the CRs. We have implementedand tested this method in the Athena MHD code, where the overall scheme is second-order accurateand fully conservative. As a first application, we describe a numerical experiment to study particleacceleration in non-relativistic shocks. Using a simplified prescription for ion injection, we reproducethe shock structure and the CR energy spectra obtained with more self-consistent hybrid-PIC sim-ulations, but at substantially reduced computational cost. We also show that the CR-induced Halleffect reduces the growth rate of the Bell instability and affects the gas dynamics in the vicinity ofthe shock front. As a step forward, we are able to capture the transition of particle acceleration fromnon-relativistic to relativistic regimes, with momentum spectrum f ( p ) ∝ p − connecting smoothlythrough the transition, as expected from the theory of Fermi acceleration. Subject headings: methods: numerical — magnetohydrodynamics — plasmas — instabilities — accel-eration of particles — shock waves INTRODUCTION
Cosmic rays (CRs) are charged particles that possessmuch higher energy than particles in the surroundingmedium. They are a dynamically important constituentin the Galaxy, where their energy density is in roughequipartition with gas pressure and magnetic pressure.They also provide key observational diagnostics for as-trophysical phenomena via non-thermal radiation. Theorigin and transport of CRs are among the forefront top-ics of astrophysics, and involve collisionless kinetic effectswith complex non-linear interactions with magnetic fieldsand background thermal plasma.Particle-in-cell (PIC) methods (e.g., Birdsall & Lang-don 2005), which iteratively move particles and updatethe electromagnetic fields, have been commonly em-ployed to study the plasma physics of CRs from firstprinciples. Fully kinetic PIC codes treat both electronsand ions as particles (e.g., Fonseca et al. 2002, Bowerset al. 2008, Spitkovsky 2005). While being the mostself-consistent method that rigorously solves for every-thing from plasma oscillations to the propagation of lightwaves, the main disadvantage is that it has to resolvethe microscopic plasma scale: the plasma (electron) skin [email protected] NASA Hubble Fellow NASA Einstein Fellow depth c/ω pe ∼ . × cm ( n/ cm − ) − / , where ω pe isthe electron plasma frequency. Moreover, the ion scale iswell separated from the electron scale. Even with arti-ficially reduced ion-to-electron mass ratio, properly cap-turing the physics at both scales still demands tremen-dous computational cost.An alternative approach, called the hybrid-PIC, treatsall ions as kinetic particles, while electrons are includedas a neutralizing, massless and conducting fluid with aprescribed equation of state (e.g., Lipatov 2002; Winskeet al. 2003; Gargat´e et al. 2007; Kunz et al. 2014a). Bycompromising the electron-scale physics, the physics atthe ion-scale can be well captured with much reducedcomputational cost, and this method has been success-fully applied to study a wide range of laboratory, spacephysics and astrophysical problems. The hybrid-PIC ap-proach again is limited by the requirement of resolv-ing the characteristic ion scale, the ion inertial length c/ω pi ∼ . × cm ( n/ cm − ) − / , where ω pi is theion plasma frequency and n is the ion number density.This scale is typically much smaller than the CR gyroradius. The latter scale is particularly relevant for reso-nant scattering of CRs, yet it becomes progressively moredifficult to accommodate using the hybrid-PIC methodas the CR gyroradius becomes much larger than the ioninertial length with increasing CR energy.Another approach is to treat CRs as non-thermal test a r X i v : . [ a s t r o - ph . H E ] J un X.-N. Bai et al.particles, embedded in a thermal plasma (gas) describedby magnetohydrodynamics (e.g. Lehe et al. 2009; Beres-nyak et al. 2011). Being test particles, CRs only pas-sively respond to the electromagnetic field carried by thegas. On the other hand, CRs are well known to influencethe dynamics of the background gas and magnetic fields,which in turn affect their propagation (e.g., Kulsrud &Pearce 1969; Skilling 1975; Parker 1992).In this paper, we consider another approach, whichwe call the MHD-PIC method. We treat only CRs asparticles using the conventional PIC technique, whilethe rest of the thermal plasma (ions and electrons) istreated as gaseous fluid described by magnetohydrody-namics (MHD). Different from the test particle approach,the CR feedback to the gas is included following a rig-orous formulation that enables the mutual interactionsbetween the thermal gas and CRs. With CRs treatedas particles, kinetic effects of CRs can be fully captured.With fluid treatment of the thermal plasma, we avoid re-solving microscopic scales, allowing us to accommodatemuch larger and even macroscopic scales in our simula-tions with only modest computational cost .We note that similar methods as we describe in thispaper have been proposed and implemented to study theCR streaming instabilities. Zachary (1987) derived theformulation from a plasma point of view (published inZachary & Cohen 1986), with CR feedback reflected inthe momentum and energy equations as external Lorentzforce and heating/cooling due to CR charge and current.Lucek & Bell (2000) used the same formulation and dis-covered in their simulations what was later termed theBell instability (Bell 2004), which is a non-resonant in-stability driven by the CR current. Reville & Bell (2012)implemented the same method and studied in detail thefilamentation instability of Bell (2005). Reville & Bell(2013) and Bell et al. (2013) adopted the same MHDformulation, but developed a Vlasov-Fokker-Planck ap-proach for CRs with approximate closure relations (in-stead of using particles) to study CR diffusion, escapeand magnetic field amplification in the upstream of col-lisionless shocks.In this work, we re-derive the formulation rigorouslyin the MHD framework, and further show that there isan additional CR feedback term which modifies the in-duction and energy equations. For sufficiently strong CRcurrents, this term, which we call the CR-induced Hallterm, may appreciably affect the growth rate of the Bellinstability and become important in the upstream fluidas it approaches the shock front. We have implementedthe MHD-PIC formalism described in this paper intothe state-of-the-art, fully parallelized MHD code Athena(Stone et al. 2008). This opens up a variety of CR-relatedproblems that can be addressed using non-linear com-puter simulations, especially the long-term evolution ofnon-relativistic shocks, to be introduced in Section 1.1. While the thermal plasma can be collisionless or weakly col-lisional, we employ the MHD framework which implicitly assumesefficient inter-particle “collisions” via plasma instabilities at mi-croscopic scales so that the plasma is well isotropized at relativelylarge scales that we consider.
Applications
One of our main applications is to study the structure,evolution and particle acceleration in non-relativistic col-lisionless shocks, which have long been considered as themost promising sites for producing CRs by means of thefirst-order Fermi acceleration mechanism (Fermi 1949;Bell 1978; Blandford & Ostriker 1978). In particular,diffusive shock acceleration (DSA) of electrons and nu-clei at supernova remnant (SNR) blast waves is expectedto produce power-law energy distributions of acceleratedparticles extending over many energy decades (e.g. Mor-lino & Caprioli 2012; Ackermann et al. 2013), which isthought to be responsible for most of the Galactic cos-mic rays (CRs), possibly up to the “knee” part of theCR spectrum ( ∼ eV).Theoretical investigation of CR acceleration in non-relativistic collisionless shocks has a long history. Mostearlier works either analytically or numerically solve theCR transport equations (e.g., Berezhko & Volk 1997;Malkov 1997; Kang et al. 2002; Amato & Blasi 2005;Zirakashvili & Aharonian 2010; Caprioli 2012) or adopta Monte-Carlo approach (e.g., Ellison et al. 1990, 1995;Niemiec & Ostrowski 2004; Vladimirov et al. 2006).While these works generally yield results consistent withobservations, they rely on quasi-linear theory for mag-netic field amplification with simple prescriptions of CRdiffusion coefficients and injection efficiency, reflectingthe poorly known microphysics of particle injection andturbulence. In particular, the turbulence that is respon-sible for scattering CRs should be largely generated bythe accelerated CRs themselves, which in turn affects theshock structure and the injection process. The complexinterplay among these physical processes calls for self-consistent numerical simulations.PIC approach is at the frontier of this research. Simu-lations have been carried out using full PIC codes (e.g.,Amano & Hoshino 2007; Riquelme & Spitkovsky 2011;Narayan et al. 2012; Guo et al. 2014), which mainly ad-dress electron acceleration, as well as hybrid-PIC codes(e.g., Quest 1988; Giacalone et al. 1992; Giacalone 2004;Gargat´e & Spitkovsky 2012; Guo & Giacalone 2013),which mainly address ion acceleration. With growingcomputing power, simulations with relatively large do-mains and relatively long duration are becoming moreaffordable (e.g., Caprioli & Spitkovsky 2013), allowingthe physical processes such as particle injection, diffusionand acceleration efficiency to be studied in greater detail(Caprioli & Spitkovsky 2014a,b,c). However, despite thesuccesses of these state-of-the-art simulations, they stillcapture only the initial stages of shock evolution and par-ticle acceleration. Larger spatial scales and longer timeevolution are needed as particles accelerate to higher en-ergies to accommodate the interactions of these particleswith turbulence. In addition, current hybrid-PIC codesdeal only with non-relativistic particles. Continued ac-celeration would in reality lead to transition into the rela-tivistic regime, giving a different power-law slope in par-ticle energy spectrum according to the Fermi accelerationtheory (from f ( E ) ∝ E − / to E − in the test particlelimit). This transition can not be captured in hybrid-PIC simulations, and remains to be confirmed. More-HD-PIC Method for CR-Gas Interaction 3over, since the total energy in CRs ( (cid:82) Ef ( E ) dE ) woulddiverge as E → ∞ , we expect the fraction of CR particles( (cid:82) f ( E ) dE ) to decrease with time as particles are accel-erated to higher energies. This reduction is importantfor understanding the long-term evolution of collision-less shocks, and again is too computationally demandingto be addressed by using the hybrid-PIC approach.Without the need to resolve plasma scales, our MHD-PIC approach can substantially outperform hybrid-PICcodes, allowing us to conduct large simulations follow-ing the long-term evolution of non-relativistic collision-less shocks with self-consistent particle acceleration. Onthe other hand, since the gas and CRs are treated sepa-rately, the main disadvantage of the MHD-PIC approach(as well as similar approaches, e.g., Reville & Bell 2013)is that injection of CR particles has to be handled withad hoc prescriptions. However, this may be remediedby properly designing controlled experiments to calibratethe injection recipes with hybrid-PIC simulations, as wedefer to future work. In this paper, as the first applica-tion of our code to study particle acceleration in shocks,we adopt a simple and easily controllable injection pre-scription to demonstrate the capability of our code tohandle highly non-linear regime of CR-gas interaction.Even with such a simple prescription, we are able to re-produce most of the essential features of shock structureand particle acceleration as in hybrid-PIC simulations.Most interestingly, we have performed one simulationwith unprecedented long duration which followed particleacceleration from non-relativistic to relativistic regimes.In addition, our formalism is well suited to study theinteraction between the CRs and the gas in the ISM, po-tentially providing microphysical input for studying CRpropagation and transport in the Galaxy (Strong et al.2007). In particular, the non-linear wave-particle inter-action of propagating CRs (Kulsrud & Pearce 1969) isexpected to play an important role in self-confinementof CRs in the Galaxy. In addition, given the dynamicalimportance of (especially low-energy) CRs in the ISM,the back-reaction from CRs may also affect the cascadeof MHD turbulence (e.g., Lazarian & Beresnyak 2006),and in turn their own propagation and transport in theGalaxy.This paper is organized as follows. We derive the for-mulation of the CR-gas system in Section 2. Implemen-tation of this formulation to the Athena MHD code isbriefly described in Section 3, with details and code testsprovided in Appendices A and B. As a special applica-tion of our formulation, we revisit the linear dispersionrelation of the Bell instability in Section 4 and highlightthe relevance of the CR-induced Hall term. In Section 5,we describe the setup of shock simulations with a sim-ple prescription of particle injection. Simulation resultspresented and analyzed in Sections 6 and 7, highlightingthe code capabilities and transition of particle acceler-ation from non-relativistic to relativistic regimes. Wesummarize and conclude in Section 8. FORMULATION
We consider the interaction between a thermal plasma(gas) and cosmic ray (CR) particles. The gas is fullyionized, non-relativistic, and is treated as a magnetized fluid. It consists of ions, which contain all the inertia,and massless electrons, which mainly conduct current.CR particles are relativistic and collisionless, so that theyare only subject to the Lorentz force. We assume thatCRs are particles with only one sign of charge (positive).While CRs do not interact with the gas directly, theyprovide indirect feedback via interactions with electro-magnetic fields. In this section, we formulate the equa-tions for CR-gas interaction. We rigorously derive theform of CR feedback to the gas in the MHD frameworkin Section 2.1, where CRs are treated as fluid, and inSection 2.2, we discuss the motion of CRs as individualparticles.
MHD Equations with CR Feedback
The effect of collisionless CRs on the thermal gas canbe well represented by the CR charge density n CR andcurrent density J CR (zeroth and first moments of mo-mentum distribution, treated as a fluid). They are re-lated by J CR = n CR u CR , (1)where u CR can be interpreted as the mean drift velocityof CR particles.CRs are subject to the Lorentz force, and the forcedensity they experience is F CR = n CR E + 1 c J CR × B , (2)where E and B are the electric and magnetic fields.The total current density in the system comprises ofcontributions from the gas (ions and electrons) and theCRs J tot = c π ∇ × B = J g + J CR = n i v i + n e v e + n CR u CR , (3)where n i and v i ( n e and v e ) are the charge density andmean velocity of the ions (electrons). We are interested inscales much larger than the Debye length, where chargeneutrality applies, n i + n e + n CR = 0 . (4)Note that being charge densities, we have n i > n e <
0. While we generally consider CR ions ( n CR > n CR can be either positiveor negative.In general, we expect | n CR | (cid:28) n i ≈ | n e | , and this is theregime where our formulation is applicable: essentiallyall electrons come from the gas, thus they are largelythermal and can be treated as a massless fluid, obeyingthe force-free condition n e (cid:18) E + v e c × B (cid:19) − ∇ · P e = 0 , (5)where P e is the electron pressure tensor. X.-N. Bai et al.Combining the equations above, we arrive at the gen-eralized Ohm’s law in the presence of CRs: E = − v g c × B + 1 | n e | c J tot × B − | n e | ∇ · P e − n CR | n e | ( u CR − v g ) c × B , (6)where v g ≡ v i is the bulk velocity of the gas (ions). Inthe above, the first term is the normal inductive termin MHD, the second corresponds to the usual Hall term,and the third is the gradient of the electron pressure.The last term is a novel term due to the presence ofCRs, which we call the CR-induced Hall term (hereafter,CR-Hall term). It describes the relative drift betweenCRs and the bulk of the plasma.Both the Hall term and the electron pressure gra-dient term typically become important only on scalessmaller than the ion inertial length c/ω pi , where ω pi =(4 πn i /ρ ) / is the ion plasma frequency, and ρ is the gasmass density. To see this, we compare the Hall termwith the induction term, and the Hall term becomesimportant when J tot /n e is comparable or larger thanthe Alfv´en velocity v A = B/ √ πρ . Replacing J tot by cB/L H , we see that the Hall term dominates on scales L H (cid:46) cB | n e | v A ≈ cω pi , (7)where we have used | n e | ≈ n i (cid:29) | n CR | . Similarly, as-suming that the electron pressure is comparable to thegas pressure (which again generally requires | n e | ≈ n i ),one can further replace ∇ · P e by ρc s /L P , where c s isthe sound speed. It is straightforward to show that theelectron pressure gradient term dominates the inductiveterm when L P (cid:46) (cid:18) c s v A (cid:19)(cid:18) cω pi (cid:19) . (8)For typical plasmas where thermal and magnetic energiesare in rough equipartition, we find L P ∼ L H ∼ ( c/ω pi ).Under the MHD framework, we are interested in scalesmuch larger than c/ω pi , thus we ignore the conventionalHall and electron pressure gradient terms. We note thatthese two terms are explicitly included in the formulationfor hybrid-PIC codes, which aim at resolving the physicson the scale of ion inertial length. In these codes, someprescription of the electron pressure tensor (i.e., electronequation of state) is needed as a closure relation to ac-count for the unresolved electron microphysics. At muchlarger scales, however, MHD assumption is expected tohold (as a closure relation), and we arrive at the idealMHD Ohm’s law with the corrections due to the CR-Hall term E = − v g c × B − n CR | n e | ( u CR − v g ) c × B . (9)Note that unlike the conventional Hall term, the CR-Hallterm does not provide a typical lengthscale above whichit can be safely neglected, and should be always kept inthe formulation. More insight into the importance of CR-Hall term canbe gained by rearranging the above equation in two ways.First, we define R ≡ n CR | n e | = n CR n i + n CR , (10)the ratio of the CR charge density to the total chargedensity (of the electrons). Note that our formulationrequires that R (cid:28)
1. Equation (9) is then equivalent to E = − [(1 − R ) v g + R u CR ] × B /c . (11)This equation has a clear physical meaning, namely, theadvection velocity of the magnetic field is the charge-density-weighted mean velocity of the composite gas-CRsystem. We can also see that the CR-Hall term becomesdynamically important when the following quantity,Λ ≡ R ( u CR − v g ) v A , (12)approaches unity. While we expect R (cid:28)
1, CRs mayachieve highly super-Alfv´enic drift velocities relative tothe gas so that Λ (cid:38)
1. Physically, from Equations (3)and (4), this requires the electron-ion drift velocity in thethermal gas to be on the order of v A or larger. Note thatthe super-Alfv´enic streaming between background elec-trons and ions may lead to plasma streaming instabilitiesat the electron scale, e.g., Weibel-like filamentation wasfound at about 10 electron skin depths in full-PIC simu-lations by Riquelme & Spitkovsky (2009). These effectsare not captured in either MHD or hybrid-PIC methods.While our formulation remains formally valid for Λ > E ≡ − v g c × B , (13)as the electric field in ideal MHD without CRs. Substi-tuting Equation (9) into Equation (2), we can re-expressthe overall Lorentz force felt by CRs as F CR = (1 − R ) (cid:18) n CR E + 1 c J CR × B (cid:19) . (14)This is an explicit manifestation of the effect of electro-magnetic feedback: the overall Lorentz force experiencedby CRs as a fluid is reduced by a fraction R comparedwith the case without the CR-Hall term. Plugging Equa-tion (14) back to Equation (9), the electric field can berewritten in a particularly simple form E = E − F CR n g , (15)where n g ≡ n i is the bulk charge density of the (ion)fluid.So far we have focused on the electromagnetic feed-back from the CRs, which determines the evolution ofthe magnetic field via the induction equation, ∂ B ∂t = − c ∇ × E . (16)HD-PIC Method for CR-Gas Interaction 5Below we discuss the rest of the MHD equations for thegas.CRs streaming in the gas also exchange momentumwith the gas. Due to charge-neutrality, the gas (ionsand electrons) possesses the opposite charge density toCRs. Correspondingly, the gas experiences exactly theopposite electric force density as CRs in any local vol-ume of the system, namely, − n CR E . In addition, asin Equation (3), the total current density in the sys-tem has contributions from both CRs, J CR , and thegas, J g . The magnetic force experienced by the gasis J g × B /c , while the ideal MHD equations adopt( ∇ × B ) × B / π = J tot × B /c . Therefore, we needto compensate for the difference given by − J CR × B /c ,which is exactly the opposite of the magnetic force expe-rienced by CRs. Combining these two facts, we see thatthe momentum feedback from CRs on the gas is sim-ply to “kick” the gas with the opposite of the Lorentzforce F CR experienced by CRs. This is an explicit man-ifestation of Newton’s third law which guarantees totalmomentum conservation.The energy feedback from CRs can be understood asthe work done by the external force (CR momentum feed-back) per unit time, namely, v g · ( − F CR ). Using previousequations (2), (11), (13), and (14), we can easily derive v g · F CR = (1 − R ) v g · ( n CR E + J CR × B /c )= (1 − R ) J CR · E = J CR · E = u CR · F CR . (17)We note that u CR · F CR can be understood as the rate ofenergy gain by CRs due to the work done by the Lorentzforce. The above equality means that the same amountof energy is lost by the gas, which is an explicit manifes-tation of energy conservation of the composite gas andCR system. Moreover, the equality with J CR · E indi-cates that the energy exchange with the gas is achievedvia Joule heating (or cooling).Based on the discussions above, the remaining of theideal MHD equations for the gas can be written as follows(in conservative form) ∂ρ∂t + ∇ · ( ρ v g ) = 0 , (18) ∂ρ v g ∂t + ∇ · ( ρ v Tg v g − B T B π + P ∗ g I )= − (1 − R )( n CR E + J CR × B /c ) = − F CR , (19) ∂E∂t + ∇ · (cid:20) ( E + P ∗ ) v g − ( B · v g ) B π + c π ( E − E ) × B (cid:21) = − (1 − R ) J CR · E = − u CR · F CR , (20)where P ∗ g ≡ P g + B / π , I is the identity tensor, and thetotal energy density is defined as E = P g γ − ρv g + B π . (21)In the above, ρ , P g are gas density and pressure, and γ is the adiabatic index. There is also a source termin the energy flux in the energy equation (20), whichaccounts for the correction to the ideal MHD Poyntingflux ( c E × B / π ) due to the CR-Hall term.Note that the CR feedback source terms are generallyproportional to n CR ∝ R . The requirement of R (cid:28) c/ω pi , andwhen the CR charge and current density satisfy R (cid:28) (cid:38) Particle Treatment of CRs
We adopted Gaussian units in our formulation in Sec-tion 2.1, which contains factors of speed of light, c , in theexpressions. This factor has no physical significance innon-relativistic MHD and can be absorbed by redefiningthe electric field and charge/current density as c E → E and J /c → J . On the other hand, our CR particles canbe fully relativistic. Therefore, we artificially specify thespeed of light C for the CR particles. For consistencywith MHD, the value of C must be chosen such that it ismuch larger than any characteristic MHD velocities. Westill keep the irrelevant factor c so as to retain consistencyand to contrast with the factor C . For individual particle j , we use u j to denote its ve-locity. We further define the vector component of itsfour-velocity to be v j ≡ γ j u j , (22)where the Lorentz factor γ j = (cid:113) C + v j C = C (cid:113) C − u j . (23)We define CR particle kinetic energy E k as the differencebetween its total energy and rest mass energy: E k,j ≡ ( γ j − C = v j γ j + 1 , (24)which reduces to standard expressions of particle energyin both non-relativistic and relativistic regimes. One may continue to use c to specify particle speed of light. Weprefer to introduce a different symbol for conceptual clarity: in thenon-relativistic MHD framework, C is a freely-variable numericalparameter rather than a physical “constant.” X.-N. Bai et al.The particle equation of motion reads d v j dt = (cid:18) qmc (cid:19) j (cid:18) c E + u j × B (cid:19) . (25)where ( q/mc ) j represents particle charge-to-mass ratio.For the electric field, instead of using its primitive expres-sion (9), we use the alternative expression from Equation(15). It has the clear advantage that the correction fromthe CR-Hall term is proportional to the momentum feed-back ( F CR ), and hence no extra computation is needed.Note that here, and essentially in all other occasions, thefactor c can be absorbed into the ( q/mc ) factor. Thisfactor incorporates the particle intrinsic properties andis the only physical parameter to distinguish different CRparticle species. IMPLEMENTATION
We have implemented the MHD formulation describedin § LINEAR ANALYSIS OF THE BELL INSTABILITY
Before proceeding to applications, we first consider thenon-resonant CR streaming instability of Bell (2004),and point out the relevance of the CR-Hall term thatwas previously neglected. This instability occurs whenthe CR drift velocity exceeds the local Alfv´en velocity,where right-handed circularly-polarized modes becomeunstable as the Lorentz force exceeds magnetic tension.This instability has unstable wavelengths much smallerthan the Larmor radius of the CR particles, and its freeenergy comes from the relative streaming of gas and CRs.The Bell instability is considered to play an impor-tant role in the upstream of high Mach number non-relativistic collisionless shocks, such as SNR shocks. Ac-celerated CRs drift with respect to the upstream fluidat velocity U s , which is of the order the shock velocityor higher and is much larger than v A . Since the typi-cal growth time is much shorter than the advection timein the CR-induced shock precursor, the instability typ-ically enters its strongly non-linear stage; the preshockmagnetic field may be effectively amplified, and the self-generated magnetic turbulence feeds back on CRs, en-hancing their diffusion. Hybrid-PIC simulations have re-cently shown that such a non-linear interplay betweenCRs and magnetic fields shapes the large-scale structureof the shock, and promotes rapid acceleration of energeticparticles (Caprioli & Spitkovsky 2014a,b,c).We note that the CR-Hall term was not included ei-ther in the original MHD derivation of the instability(Bell 2004), or in the kinetic approach by Amato & Blasi(2009). Here we re-derive the dispersion relation thatproperly includes the CR-Hall term. Suppose the back-ground gas is uniform with density ρ and zero veloc-ity v g = 0 (upstream reference frame). Let B be thebackground magnetic field, and J CR = n CR U s be theexternal current density provided by the CR particles.Both of them are along the ˆ x direction. Now we con-sider incompressible perturbations in the gas, assumingall physical quantities to evolve as exp i( kx − ωt ) (notethat the wave vector k is also along the ˆ x direction) anduse u , b to denote first-order perturbations of gas veloc-ity and magnetic field. For incompressible modes, bothquantities are in the transverse direction. From the mo-mentum equation (19) (or more explicitly see Equationof (1) of Bell 2004), we have − i ωρ ( u × k ) = +i( k · B )( b × k ) − (1 − R )[( J CR × b ) − n CR ( u × B )] × k , (26)where for conciseness we have ignored factors contain-ing c and 4 π , which will be irrelevant to the results.Also note the factor (1 − R ) that comes from EquationHD-PIC Method for CR-Gas Interaction 7(14). The perturbation equation for the induction equa-tion (16) reads − ω b = (1 − R )( k · B ) u − R ( k · U s ) b . (27)The above equation gives the relation between b and u .Combining the two equations, and noting that k , B and J CR are parallel, we find after some algebrai[ ω − ωkRU s − (1 − R ) k v A ]( b × ˆ x )+(1 − R ) kJ CR v A B (cid:18) − ωkU s (cid:19) b = 0 . (28)Note that all terms containing R are due to the CR-Hall term, and the most significant modification lies inthe term ωkRU s , which leads to a non-negligible shift inwave frequency. Our Λ-parameter defined in Equation(12) corresponds to Λ = RU s /v A .We provide the full derivation of the dispersion relationin Appendix C. In the following, we consider the limit R (cid:28) ω ± (cid:15)k v A ) = f (cid:18) k ± f k (cid:19) v A − ( k v A ) (cid:18) f − (cid:15) (cid:19) , (29)where the two ± signs must be taken to be the same, andwe have defined˜ ω ≡ ω − (Λ / v A k , k ≡ J CR B , (30) f ≡ / , (cid:15) ≡ v A U s . (31)Note that for highly super-Alfv´enic shocks, (cid:15) (cid:28)
1. Whenthe CR-Hall term is not included or for Λ (cid:28)
1, we have f = 1, and the dispersion relation (29) reduces to Equa-tion (7) of Bell (2004) (regime II, since we treat CRs asa fixed external current). In the general case, the systembecomes unstable when (cid:15) < /f . The most unstablemode corresponds to k m = ± k /f , with fastest growthrate of the order k v A / √ f .Now we discuss the role of the CR-Hall term, set-ting (cid:15) = 0, valid for highly super-Alfv´enic shocks. Ourlinear analysis of the Bell instability leads to qualita-tive changes in two aspects. First, the most unstablewavenumber k m is reduced by a factor of about f . With-out the CR-Hall term, we have k m = k ∝ J CR ∝ Λ,which increases with J CR without bound. Includingthis term, and considering the limit Λ (cid:29) f ∼ ( RU s / v A ) ) we have k m = k /f ∼ J CR v A B R U s ∼ (cid:18) cω pi (cid:19) − . (32)Matching the two limits of Λ (cid:28) (cid:29)
1, we see thatas J CR (hence Λ) increases from zero, k m ∝ J CR , andthe most unstable wavelength decreases until it reachesthe scale of ion inertial length, at about Λ ∼
2. Fur-ther increasing J CR would increase the most unstablewavelength. Therefore, the most unstable mode neverfalls below the scale of c/ω pi . This result may have im-portant consequences for the resonant scattering of CRs in the vicinity of the shock front. Second, the fastestgrowth rate ∼ k v A / √ f , is reduced by a factor of about √ f . Without the CR-Hall term, the fastest growth rateincreases with J CR without bound since k ∝ J CR . In-cluding this term, and again in the limit Λ (cid:29)
1, we have ω max ∼ J CR v A B RU s ∼ Ω c . (33)Therefore, the maximum growth rate saturates at Ω c ≡ qB /m i c , the ion cyclotron frequency of the backgroundplasma.We note that using full kinetic PIC simulations,Riquelme & Spitkovsky (2009) showed that when Λ (cid:38) c and migrates to larger scales at Λ (cid:29)
1, itis not surprising that the Bell mode is overtaken by thefilamentation mode, which may not be captured in ourMHD-PIC formulation.We will show in our shock simulations that Λ can ap-proach order unity in the vicinity of the shock front.While this already corresponds to the non-linear stage ofthe Bell instability, our derivations in the linear regimestill provide a useful quantitative assessment on the im-portance of the CR-Hall effect on the gas dynamics inthis region. PARTICLE ACCELERATION IN COLLISIONLESSSHOCKS: SIMULATION SETUP
As introduced in Section 1.1, we conduct a proof-of-concept study of the evolution and particle accelerationin non-relativistic collisionless shocks to demonstrate ourcode performance, where the non-linear development ofthe Bell instability in the shock upstream plays a vitalrole (e.g., Bell et al. 2013; Caprioli & Spitkovsky 2014b).The basic setup of the shock problem involves collid-ing a highly super-Alfv´enic flow into a conducting andstatic wall (left boundary in the Figures). The shockforms instantaneously, propagating away from the wall.After the collision, the downstream gas has zero longi-tudinal velocity, hence the simulation takes place in thedownstream frame. Let ˆ x denote the direction of theshock, and the conducting wall is located at x = 0. Thecolliding gas (upstream) moves from the right ( x > ρ , pressure p and velocity − v ˆ x . We consider parallel shocks withmagnetic field along the ˆ x direction, whose strength is B . For convenience, we define the shock Alfv´enic Machnumber to be M A ≡ v /v A , where v A = B / √ ρ isthe background Alfv´en velocity . The angle between thebackground field and shock normal (ˆ x ) is denoted by θ which, together with M A , are the main shock parameters(the sonic Mach number is comparable to M A ).For the purposes of this paper, we only consider highMach number, parallel shocks with M A = 30. Forhighly super-sonic and super-Alfv´enic shock, we have More strictly, shock Mach number is defined in the shockframe, our definition simply follows the convention of Caprioli &Spitkovsky (2013).
X.-N. Bai et al. ρ v (cid:29) P + B /
2, which leads to shock jump condi-tions (expressed in the downstream frame) under idealMHD ρ = rρ , P = rr − ρ v ,B x = B x , v sh = v / ( r − . (34)where r ≡ ( γ + 1) / ( γ −
1) is the shock compression ra-tio, v sh is the velocity of the shock in the simulation(downstream) frame. For adiabatic index γ = 5 /
3, wehave r = 4 for strong shock as usual. We use subscripts“0” and “1” to denote upstream and downstream quan-tities measured in the simulation (downstream) frame,hence v = 0 by definition. The shock converts thekinetic energy in the upstream gas into thermal en-ergy. The gas temperature can be defined as p/ρ , hence T = P /ρ , while T = v / ( r −
1) is independent of T for strong shocks. We also define the shock kineticenergy E sh ≡ v /
2, which characterizes the energy ofdownstream thermal particles (per unit mass) and willbe used as a scale to normalize the energy of CR parti-cles.The natural units for length, time and velocity arethe ion inertial length c/ω pi = ( q/mc ) − ρ − / , ion cy-clotron frequency Ω c = ( q/mc ) B , and the Alfv´en ve-locity v A = B ρ − / . They are all set to 1 in code units,with ρ = B = ( q/mc ) = 1. We also set T = 1 (or P = 1), which is irrelevant as long as T (cid:28) v . We aremainly interested in highly super-Alfv´enic shocks so that v (cid:29) v A (we choose v = 30). This means that in thedownstream frame, the gyro-radius of a thermal particlearound background field is much larger than c/ω pi , by afactor of ∼ M A . Initially, we focus on particle acceler-ation in the non-relativistic regime, and set the particlespeed of light to C = 10 (cid:29) v (so that particles are notable to achieve relativistic velocities within the durationof simulations). In addition, we also run one simulationwith much reduced speed of light C = 10 √ v = 424 . E k, trans ∼ C / E sh , which is optimalfor studying particle acceleration trasitioning to the rel-ativistic regime.Being a higher-order Godunov code, the MHD integra-tor in Athena works very well for shock capturing. Withuniform and homogeneous upstream medium, the shockis captured with only 2-3 grid cells using the CTU inte-grator with third-order reconstruction. The microphysicsunderlying the sharp transition in the MHD shock isthe efficient particle thermalization process, mediatedby electromagnetic turbulence as a result of plasma in-stabilities on microscopic scales ( (cid:46) c/ω pi ). While suchplasma micro-instabilities are not captured in the MHDframework, we expect the scattering of CR particles to bedominated by turbulence on larger scales in highly super-Alfv´enic turbulence, which is captured in the MHD-PICapproach when streaming instabilities such as the Bellinstability are fully developed.Here we describe our simple prescription for particleinjection. In every time step, we keep track of the loca- TABLE 1List of main shock simulation runs
Run Domain size resolution η C L x × L y ( c/ω pi ) ( c/ω pi per cell)R1 (1 . × ) × − R2 (1 . × ) × × − R4 (1 . × ) × × − R2-hr (1 . × ) × × − R2-REL (3 . × ) × × − . M A = 30. tion of the shock front x s by computing the transverselyaveraged profile of v x along ˆ x . This 1D profile is fur-ther smoothed with a Gaussian kernel whose width is 4grid cells. The shock front x s is identified as the loca-tion where v x is reduced by about 40% with respect tothe upstream value. With x s identified, the amount ofmass traversed by the shock can also be obtained. Weinject a spatially uniform, isotropic distribution of CRparticles at the shock front in the co-moving frame withthe shock, whose energy is set to be 10 E sh in this frame.The amount (mass fraction) of CR particles we inject is asmall and fixed fraction η of the amount of mass swept bythe shock. Note that E k = 10 E sh is approximately thethreshold energy for a particle to be considered as non-thermal and to participate in the DSA process (Caprioli& Spitkovsky 2014a; Caprioli et al. 2015). In brief, η isour only parameter for particle injection, which encap-sulates our ignorance of the physics of particle injection.Physically, the injected particles originate from the gasthat is newly processed by the shock. While we do notcapture the injection physics with the simple prescrip-tion, we do compensate for the mass, momentum andenergy of injected particles by subtracting them off fromthe gas at each timestep over certain width at the down-stream side of the shock , so that the total mass, mo-mentum and energy in the full system are conserved. Fortypical injection fraction η = 2 × − (see below), be-cause the energy of injected particles only amounts to avery small fraction of the shock energy ( ∼ j CR ). Since the growth rate of the Bell We further define the location of the tail end of the shock x t based on the 1D smoothed pressure profile (where P achieves2 P / x t and x s . HD-PIC Method for CR-Gas Interaction 9
Fig. 1.—
Snapshots of gas density at four different times (as indicated) from our fiducial run R2 which show the evolution of theshock. Only a small fraction of the simulation box is shown which is adjusted to cover the vicinity of the expected shock location at x = v sh t = v t/ instability is proportional to j CR , this transient CR flowwould lead to growth of the Bell instability by a finiteamplitude, as a result of the unrealistic initial condition.To remedy this artifact, we record the injection time t i for every CR particle. All particles with t i < − c (which is more than sufficient for the Bell instability tobe fully developed near the shock front) are removed attime t = 2 t i . Following this procedure, we expect ourresults to be largely unaffected by initial conditions at t > − c .We perform 2D shock simulations with fiducial domainsize L x = 1 . × ( c/ω pi ) and L y = 3000( c/ω pi ). Weuse conducting boundary condition at x = 0, and in-flow boundary condition at x = L x , where all quanti-ties are fixed at the initial value. With M A = 30, hence v sh = 10, we run the simulations up to time t = 3000Ω − c so that at the end of the simulation, the gas initially atthe right boundary just reaches the shock. We choosefiducial numerical resolution of 12 c/ω pi per cell, whichwill be shown to well capture the most unstable modesof the Bell instability. This is also much smaller thanthe Larmor radius of the lowest energy CR particles or-biting background field (for E k = 10 E sh particles, theLarmor radius R L ∼ c/ω pi ). We also perform a reso-lution study with up to four times higher and two timeslower resolution (see Section 6.4), while we mainly dis- cuss the run with twice resolution. Note that in hybrid-PIC simulations, the typical numerical resolution is 2cells per c/ω pi in order to properly capture the essentialmicrophysics (e.g., Gargat´e & Spitkovsky 2012; Capri-oli & Spitkovsky 2013). With the MHD-PIC approach,the gain in efficiency is tremendous. For the simula-tion with reduced speed of light C = 10 √ v = 424 . L x =3 . × ( c/ω pi ), L y = 4 . × ( c/ω pi ), and run thesimulation to t = 11520Ω − c .We choose η = 2 × − as the standard CR massfraction in our simulations. This is motivated by thenumbers obtained from hybrid-PIC simulations (Capri-oli & Spitkovsky 2014a). In these simulations, it wasfound that in a parallel shock, about ξ ≈
15% of theshock kinetic energy is converted to non-thermal parti-cles. With η = 2 × − , the energy contained in ourinitially injected particles is about 3% of the shock en-ergy. This is sufficiently small so as not to affect theshock structure, leaving enough room for these particlesto further gain energy from the Fermi acceleration pro- In addition, the particle timestep in our MHD-PIC simula-tions is mainly constrained by resolving the Larmor time and grid-crossing time, and with Υ = 0 . t > . − c is more than 6 times larger thanused in hybrid-PIC. λ m ∼ π/ ( ηM A ) ∼ c/ω pi ), which is well resolvedwith ∼
13 cells. Larger λ m is expected further upstreamof the shock due to the reduction of net CR current den-sity. As a parameter study, we also perform simulationswith η = 10 − and η = 4 × − . For all simulations,the number of particles we inject is equivalent to 4 par-ticles per cell at background density ρ , or 16 particlesper cell on average in the downstream. Note that de-spite the reduced spatial resolution, the number densityof CR particles in physical units is comparable to thatin typical hybrid-PIC simulations (where spatial resolu-tion is much higher but only a very small fraction of theparticles enters Fermi acceleration), and we have verifiedthat simulation results converge with respect to numberof injected particles. In the simulations, we allow a max-imum of 5 particles steps per MHD step to speed up thecalculation.In sum, we mainly discuss 5 simulation runs: threeruns with C = 10 (cid:29) v using standard resolutionwith η = 10 − , × − and 4 × − , one high res-olution run, and one run with reduced speed of light C = 10 √ v = 424 .
3. These runs are summarized in Ta-ble 1, and each run is assigned a run name as indicated.We consider run R2 as the fiducial simulation and otherruns as variations. All runs except for run R2-REL arerelatively cheap computationally, with our fiducial runR2 taking about 2700 CPU hours on a HP Beowulf clus-ter with 3.47 GHz Intel Westmere cores, yet the simula-tion box size and duration are already much larger thanin the state-of-the-art hybrid-PIC simulations. PARTICLE ACCELERATION IN COLLISIONLESSSHOCKS: NON-RELATIVISTIC REGIME
We focus on runs with C (cid:29) v in this section, where allparticles remain non-relativistic at all times. We beginby discussing simulation results from our fiducial run R2in Sections 6.1 and 6.2, followed by parameter study inSection 6.3. Shock Evolution and Structure
In our simulations, the Bell instability is first excitedin the shock upstream due to the streaming CRs, andexhibits clear signature of circularly polarized modes intransverse magnetic and velocity fields, as expected fromthe linear eigenvector shown in Appendix B.2 (EquationB2). Growth into the non-linear stage leads to strongdensity variation and magnetic field amplification. InFigure 1, we show snapshots of gas density at four dif-ferent times from our fiducial run R2. Note that only asmall fraction (15%) of the simulation box in the vicinityof the shock front is shown. While the Bell instability isessentially incompressible in the linear regime, the gas isgradually evacuated to produce cavities towards the non-linear regime as a result of the filamentation instability(Bell 2005; Reville & Bell 2012). The density structurefound here is qualitatively very similar to that found inhybrid-PIC simulations of Caprioli & Spitkovsky (2013),justifying the validity of our MHD-PIC approach. Wenote that over relatively long time evolution, the systemdevelops larger and larger structures, which can only be accommodated with large transverse simulation domainsize adopted in our simulations (as opposed to typicaldomain size of (cid:46) c/ω pi in hybrid-PIC simulations).In Figure 2, we further show snapshots of the gas den-sity and magnetic field strength at time t = 2400Ω − c from all our simulations with non-relativistic particles(all runs except R2-REL). We focus on our fiducial run(second row) in this section. Accompanied by the Bell in-stability with cavitation and filamentation, we find thatupstream magnetic fields are amplified by a factor ofup to ∼ − t =2300Ω − c to t = 2500Ω − c (with proper positional shift toaccount for shock motion). Focusing on run R2 here, wesee that prior to the shock front, the mean field strengthhas increased by a factor of ∼
2. The downstream meanfield strength is about 6 −
10 stronger than the initialfield B = 1, which is again consistent with hybrid-PICsimulations of Caprioli & Spitkovsky (2014a). We alsonote that while the downstream profile of mean fieldstrength is highly time variable, which is also reflectedin its large spatial variations, the profile shown in theFigure is rather typical.We next show, in the upper panel of Figure 4, thetransversely averaged profile of longitudinal CR currentdensity j CR ,x . It is measured in the co-moving frameof the gas (based on transversely averaged gas velocityin ˆ x ), and is the source of free energy in the Bell in-stability. The measured value of j CR ,x slowly decreasestowards the far upstream of the shock to the level of < .
01 in code units, a factor of 8 smaller than the ini-tial injected CR current density (4 / M A η ∼ .
08. Thisis because most of the particles are scattered back down-stream within a short distance (diffusion length) fromthe shock front, and the far-upstream CR current is car-ried by escaping particles accelerated to much higher en-ergies whose number density is much smaller (see thetop panels of Figure 5 discussed in the next subsectionfor more information). These features are all consis-tent with hybrid-PIC simulations (Caprioli & Spitkovsky2014b). We can infer the most unstable wavelengthfor the Bell instability away from the shock front to be λ m ∼ πB /j CR ,x ∼ . × ( c/ω pi ), which is very wellresolved in our simulations. This scale is also consistentwith the typical scale of density and magnetic fluctua-tions in the shock upstream present in Figure 2. To-ward the shock front, the value of j CR ,x increases to itspeak value of near 0 .
1, which can be considered as thesum of injected CR current and the current from highenergy CRs. It drops quickly to about zero in the down-HD-PIC Method for CR-Gas Interaction 11
Fig. 2.—
Simulation snapshots of gas density (left) and total magnetic field strength (right) at t = 2400Ω − c from four of our simulationruns R1, R2, R2-hr and R4 with increasing injection efficiency η as labeled in each panel. The third panels from top correspond to thehigh resolution run. The shock is approximately located at x = 24000. ω pi ] B t o t η =10 −3 η =2 × −3 η =2 × −3 , hi−res η =4 × −3 Fig. 3.—
Transversely averaged profiles of magnetic field strengthin the vicinity of the shock front around time t = 2400Ω − c for fourof our simulation runs R1, R2, R2-hr and R4. Fig. 4.—
Results from our fiducial simulation R2 at t = 2400Ω − c .Top: transversely averaged x -component of the CR current densityprofile along ˆ x in the co-moving frame. Bottom: snapshot of Λdefined in Equation (12) at t = 2400Ω − c , which measures theimportance of the CR-Hall term. Vertical dashed lines in the twopanels indicate the location of the shock front. stream, since no relative motion between gas and CRs isexpected.In Sections 2 and 4, we have highlighted the poten-tial importance of the CR-Hall term and its relevance tothe Bell instability in collisionless shocks. In the bot-tom panel of Figure 4, we further show the map of Λdefined in Equation (12) at t = 2400Ω − c , which mea-sures the importance of the CR-Hall term. The spatialdistribution of Λ is very non-uniform. Comparing withthe corresponding panels in Figure 2, we see that Λ is thelargest in the density cavities, and achieves order unitynear the shock front. We also see that at up to ∼ ion inertial lengths ahead of the shock, Λ already hasnon-negligible deviations from zero. This result justifiesthe necessity of including the CR-Hall term to capturethe role of streaming CRs on the evolution of magneticfields. On the other hand, since Λ ∝ j CR ,x by definition,we have Λ (cid:28) Particle Acceleration
Since we inject particles with sufficiently high energy E k = 10 E sh into the shock upstream, most of the in-jected CR particles directly enter the DSA process. Aslong as the scattering process is isotropic in the co-moving frame, the expected spectrum of accelerated par-ticles does not depend on the details of the scattering,but only on the compression ratio r , with the momen-tum spectrum of accelerated particles to be f ( p ) ∝ p − q ,where q = 3 r/ ( r − (cid:39) r (cid:39)
4. The energy spectrumis related to the momentum spectrum via f ( E ) = 4 πp f ( p ) dpdE . (35)For non-relativistic particles considered here, the energyspectrum is expected to have the form f ( E ) ∝ E − / .In the bottom panel of Figure 5, we show the timeevolution of the downstream particle energy spectrum2 X.-N. Bai et al. Fig. 5.—
Energy spectrum (in dimensionless form Ef ( E )) ofthe CR particles from our fiducial run R2. The top panel showsthe energy spectrum as a function of x in units of E sh = v / t = 2400Ω − c . The bottom panel shows the time evolution ofthe downstream particle energy spectrum, as marked with differentcolors indicated by the color table. The spectrum is extracted byaveraging through a layer at a fixed distance behind the shockfront. As an example, vertical black dashed lines in the top panelsmark the location of the shock front, and the vertical white dashedlines mark the layer where the downstream spectrum is extracted.The black dashed line shows the thermal energy spectrum withtemperature T = 0 . T . for our fiducial run R2. At early times before turbulenceis fully developed ( (cid:46) − c ), the particle spectrum re-flects the initial injected energy distribution, which peaksat E = 10 E sh . The spectrum then significantly broad-ens, extending substantially to the high energy side witha very small fraction scattered to low energy side. A high-energy power-law tail is gradually built up, with spectralslope consistent with − /
2. The high-energy tail extendsto higher and higher energies with time, with a roughlyexponential energy cutoff. We note that the normaliza-tion of the downstream spectrum appears to vary slightlyin time and space, which reflects the stochasticity of theshock. For example, by performing additional simula-tions with different random seeds, we quote a relativeuncertainty of ∼
20% in the normalization.We also show in black dashed line the expectedenergy spectrum for a thermal particle distribution(Maxwellian) in the downstream Ef ( E ) = 4 × (2 / √ π )( E/E th ) / e − E/E th , (36)where we choose E th = 0 . T (note that T = v / E sh to about 10 E sh (Caprioli & Spitkovsky 2013). Thepresence of the supra-thermal particles is related to theparticle injection process that bridges the gap betweenthermal and non-thermal particle populations (Caprioliet al. 2015). While this piece of physics is currently miss-ing in our simple injection prescription, our current pre-scription may serve as a first approximation to mimic theinjection process, with injected E = 10 E sh particles pre-sumably fed by the supra thermal particle population.We plan to develop and employ more realistic injectionprescriptions in the near future.In the upper panel of Figure 5, we show the spatialdistribution of the particle energy spectrum at relativelylate stage of the shock with t = 2400Ω − c . We see thathigh-energy CRs with energies ∼ E sh or higher pen-etrate into the shock upstream and provide the source ofCR current to drive the Bell instability. The lower pro-trusion at E k ∼ E sh into the shock upstream is mainlydue to the artificial injection procedure adopted here,and we see that the initially injected CR particles aremostly confined in regions close to the shock front. Wemeasure the particle energy spectrum at a fixed distanceof about 2400 c/ω pi behind the shock, as indicated in theFigure. As discussed in Caprioli & Spitkovsky (2014a),this allows the shape of the particle energy spectrum tobe fully settled.We have also examined the evolution of the maximumparticle energy. While we can always identify the parti-cle with the maximum energy, the time evolution of thisenergy can be rather noisy due to random scattering anddue to rapid changes when such particle leaves the do-main. Here, we instead consider the following definitionbased on the full energy spectrum¯ E max = (cid:82) E n +1 f ( E ) dE (cid:82) E n f ( E ) dE , (37)where n is an integer number of our choice. Notethat if the energy distribution function takes the form f ( E ) ∝ E − m exp( − E/E cut ), then this integral yields¯ E max ≈ ( n + 1 − m ) E cut . Empirically, we choose n = 6,which corresponds to ¯ E max being ∼ − . E cut if the cut-off is exponential. We also find that the number obtainedthis way is approximately half the absolutely maximumparticle energy.The time evolution of ¯ E max is directly related to therate of particle diffusion across the shock, and henceclosely probes the properties of the MHD turbulence self-generated by streaming CRs. This is shown in Figure6. Excluding the time before t = 960Ω − c (which is ex-pected to be affected by initial conditions), we see that¯ E max increases approximately linearly with time. Thisis consistent with efficient turbulent diffusion/scatteringof particles with diffusion coefficient D p ∝ v p R L , where v p and R L are particle velocity and gyro-radius (Capri-oli & Spitkovsky 2014c). Being a proof-of-concept studywith simplified injection prescription, performing morequantitative analysis on the properties of magnetic tur-bulence and particle diffusion is beyond the scope of thiswork, but we expect the detailed analysis performed inCaprioli & Spitkovsky (2014b,c) to hold in general.Based on the energy spectrum, we can measure the ef-HD-PIC Method for CR-Gas Interaction 13 Ω c−1 ) w e i gh t ed E m a x / E s h η =10 −3 η =2 × −3 η =2 × −3 , hi−res η =4 × −3 Fig. 6.—
Time evolution of the weighted maximum particle en-ergy ¯ E max from four of our simulation runs R1, R2, R2-hr and R4with increasing CR injection efficiency η . Also shown are linearfits to the data with t > − c . ficiency of particle acceleration ξ by calculating the frac-tional energy residing in the CR particles compared withthe shock kinetic energy ( mv sh / − c of our simulations,we find ξ ≈
13% from our fiducial run R2 (with relativeuncertainty of ∼ Parameter Dependence
The outcome of the simulations mainly depends on theprescribed CR injection fraction η , an artificial param-eter of the simulations. Snapshots of the shock struc-ture for all simulations with non-relativistic particles at t = 2400Ω − c are shown in Figure 2. For relatively smallinjection fraction η = 10 − (run R1), the CR currentstreaming through the upstream, and hence the growthrate of the CR-driven instabilities, are relatively small.Correspondingly, it takes longer for cavitation and fila-mentation to develop, with smaller density contrast be-fore reaching the shock front. While turbulence is gener-ated in the upstream, the shock profile remains sharp atall times in our simulation.When increasing η , cavitation becomes stronger, lead-ing to higher density contrast and stronger magnetic fluc-tuations in the shock upstream. As a result, the shockfront is more disturbed, leading to thicker shock transi-tion layer (which is more evident in the temperature pro-file not shown in the plots). For run R4 with η = 4 × − ,violent turbulence already dominates the shock near up-stream, and the shock front substantially smoothed butis still identifiable. We have also performed simulationswith η = 5 × − , and find that the shock front is somuch disrupted that it becomes difficult to identify the − − − E/E sh E f ( E ) E − Maxwellian (cid:100) =10 − (cid:100) =2 × − (cid:100) =2 × − , hi − res (cid:100) =4 × − Fig. 7.—
CR particle energy spectrum at t = 3000Ω − c in theshock downstream from our four simulation runs R1, R2, R2-hrand R4. The black dashed line shows the thermal energy spectrumwith temperature T = 0 . T . location of the shock. In realistic shocks, such strongmodification of shock structure is expected to lead to areduction of the injection efficiency, and we conclude thatimposing a constant η (cid:38) × − leads to an unphysicalover-injection.Different levels of magnetic field amplification in theshock upstream can be clearly identified in Figure 3,where the mean field strength near the shock front onthe upstream side increases with increasing η , as a resultof stronger upstream CR current. Correspondingly, thedownstream field strength also increases with η withinthe duration of our simulations.Particle acceleration is very efficient in all simulations.Figure 7 shows the cross-comparison of downstream par-ticle energy spectrum measured at t = 3000Ω − c , mea-sured from our four simulations with increasing η . For η (cid:46) × − , the particle energy spectra follow the f ( E ) ∝ E − / profile before the cutoff, while with rela-tively large η = 4 × − , the spectra slope is less con-sistent, especially at higher energies. Clearly, the totalenergy contained in the CR particles increases with theinjection parameter η , and we quote 6% and 20% as theacceleration efficiency ξ from runs with η = 10 − and4 × − , respectively. The former is smaller than ob-tained from hybrid-PIC simulations, indicating under-injection, while the latter is higher, indicating over-injection. The deviation from the expected spectralshape in run R4 may be the result of over-injection.In all cases, the maximum particle energy increasesapproximately linearly with time, as seen from Figure6. The rate of increase is shallower with smaller η , sug-gesting less efficient particle diffusion. This is in linewith the fact that there is less magnetic field amplifi-cation/fluctuation with smaller η , which leads to largerdiffusion coefficient. Increasing η to 4 × − , the rate atwhich ¯ E max grows follows more closely with our fiducialrun, and also shows larger fluctuations. This also re-flects the saturation of magnetic field amplification andparticle acceleration efficiency, again suggestive of over-injection.4 X.-N. Bai et al.Overall, our parameter study suggests that the sharp-ness/smoothness of the shock depends on the prescribedparticle injection efficiency η , which also determines theefficiency of particle acceleration ξ . Our fiducial choice of η = 2 × − appears to yield results most consistent withhybrid-PIC simulations, at least within the duration ofour runs. Over-injection of particles would over-smooththe shock. Since this is not observed in self-consistentsimulations, particle injection must be suppressed whenthe shock transition layer is over-smoothed, so that in re-ality the shock structure remains reasonably sharp, andthe particle acceleration efficiency should be kept to bewithin a certain level (cid:46) Convergence
Our high-resolution run R2-hr is in most aspects verysimilar to our fiducial run R2. Comparing their shockstructures in Figure 2, we see that run R2-hr developscavities at similar locations relative to the shock front asthe run R2. The density contrast in these cavities andthe level of magnetic field amplification are also similarbetween the two runs. While the shock transition regionin the high resolution run appears to be thinner than inthe low resolution run at this particular snapshot , this ismostly due to the very dynamic nature of the shock withhighly chaotic turbulence. Over longer timescales, thetwo runs show no significant difference in shock appear-ance. More definitive evidence of convergence is revealedin Figures 6 and 7, where we see that the particle energyspectrum (hence acceleration efficiency), together withthe evolution of maximum particle energy, agree quanti-tatively between the low and high resolution runs. Thisin turn implies that the particle diffusion coefficients arequantitatively similar between the runs.We have also computed the spectral energy distribu-tion F ( k x ) of transverse magnetic fluctuations followingthe approach of Caprioli & Spitkovsky (2014b) (see theirEquation (3)), and show the results in Figure 8. Here, F ( k x ) is dimensionless and measures the energy densityof magnetic fluctuations per logarithmic wavenumber,normalized to background field energy density. To betterdemonstrate convergence, we performed two additionalruns with half the resolution of R2 and twice the res-olution of R2-hr, and hence in total we have four runswith resolution spanning from 24 to 3 c/ω pi per cell. Wecompute F ( k x ) at fixed distances at the upstream anddownstream sides of the shock for each run (see the cap-tion of Figure 8 for details), averaged over time intervalbetween 2160 and 2400 Ω − c .We see that in the upstream, the Bell instability iswell resolved at all resolutions (modulo uncertainties inthe overall normalization as discussed previously). Inall cases, the spectrum peaks at k x ∼ . c/ω pi ) − ,corresponding to the fastest growing Bell mode aheadof the shock . In the downstream, the lowest resolu-tion simulation have somewhat less power and may beunder-resolved (but can also be due to larger noise and This corresponds to longer wavelength than the naive estimatein Section 5, since most injected particle are scattered back shortlywithout reaching large distance into the upstream.
Fig. 8.—
Dimensionless spectral energy distribution F ( k x ) oftransverse magnetic fluctuations from our resolution study of shocksimulations, averaged over a time period between 2160 and 2400Ω − c . Upper lines correspond to downstream spectrum, computedover a distance of 2400 c/ω pi in x starting from 2000 c/ω pi behindthe shock. Lower lines correspond to upstream spectrum, com-puted over a distance of 2400 c/ω pi in x starting from 4000 c/ω pi ahead of the shock. Different colors correspond to different resolu-tions, as indicated by the legend, where black solid lines correspondto our fiducial resolution run R2. The rightmost point of each linecorresponds to grid scale of the respective runs. variability since there are fewer resolution elements com-pared to higher resolution runs). For our fiducial run R2,the power spectrum matches well with higher resolutionruns at small k x . Increasing resolution, the spectrumextends to higher k x . We note that for typical grid-basedMHD code, ∼
16 cells are required to properly resolve asingle wave mode without strong numerical dissipation(e.g., see Figures 7 and 12 of Stone et al. 2008 and dis-cussions therein). This is exactly the scale where thepower spectrum of our fiducial run R2 starts to deviatefrom higher-resolution runs at k x ∼ . c/ω p ) − . Wealso see from the highest resolution run that interest-ingly, the downstream spectrum peaks at a scale aboutfour times smaller than the scale of the upstream spec-tral peak. This may be simply understood as a result ofshock compression.Overall, since CR acceleration relies on the turbulentpower at larger scales, which are well resolved in all runs,this leads to consistent calculation of the particle spec-trum and of the maximum energy achieved by acceler-ated particles (see Figures 6 and 7). Our convergencestudy thus gives us confidence that our fiducial resolutionof 12( c/ω pi ) per cell is sufficient to capture the essentialphysics. PARTICLE ACCELERATION IN COLLISIONLESSSHOCKS: TRANSITION TO RELATIVISTIC REGIME
In this section, we discuss our last run R2-REL, wherewe use a reduced particle speed of light C = 10 √ v While the spectra do not line up exactly at lowest k x , thisis most likely due to small number statistics of long-wavelengthmodes. HD-PIC Method for CR-Gas Interaction 15to study particle acceleration with transition from non-relativistic to relativistic regime. Note that the initialparticle energy of E k = 10 E sh corresponds to typicalparticle velocity of v ∼ . v ∼ . C . The correspond-ing γ ≈ . ≈
1, hence particles can still be consideredas non-relativistic. Transition to the relativistic regimeoccurs approximately when E k, trans ∼ C / E sh ,where γ ≈ .
5. This way, we expect the particle en-ergy spectrum to span nearly one decade before becom-ing relativistic. We run the simulations sufficiently longto achieve another decade of energy span in the relativis-tic regime . Shock Structure and Evolution
In Figure 9, we show the time evolution of the gas den-sity in the vicinity of the shock front. While the over-all sequence shares similarities with our fiducial run R2shown in Figure 1, there are several points worth furtherdiscussion.First, we clearly see that the waves and filamenta-tion/cavity structure in the shock upstream progress tolarger sizes over time. This is the main reason that wehave enlarged our box size to 4800( c/ω pi ) to accommo-date the growing structures. Near the end of the run(bottom panel), the size of the largest structure is al-ready a significant fraction of the transverse box size.Our additional test simulations with smaller transversebox size tend to saturate into a state with a single domi-nant transverse mode in the shock upstream, which mayaffect the scattering of energetic CR particles.Second, we notice that the shock upstream in this runappears to be less perturbed and has a sharper shocktransition layer compared to our fiducial run R2 (e.g., attime t = 2400 − − c ). This is mainly due to theuse of a smaller speed of light than in run R2. Note thatthe turbulence in the shock upstream is mainly excitedby the escaping CRs, whose contribution is mainly dueto the net CR current density. However, with a reducedspeed of light, the CR velocity saturates at a smallervalue of C : higher-energy CR particles become relativis-tic and do not contribute as much current as they doas non-relativistic particles. This effect becomes morepronounced with time, since the particles that penetrateinto the shock upstream have higher energies.Third, we notice that at later times the shock ve-locity slightly deviates from the standard jump condi-tion of v sh = v /
3, dropping by about 10% from time t = 4800Ω − c to 10800Ω − c . As a result, the shock frontlags from the expected position (Figure 9), and the shockcompression ratio density is r ≈ .
3, slightly larger thanexpected from standard jump condition (Equation 34).This is a direct consequence of channeling of a sizablefraction of the upstream ram pressure into acceleratedparticles that do not feel the standard entropy jump atthe shock, as discussed in § The time we terminate the simulation corresponds to the max-imum time when the flow that reaches the shock front is unaffectedby the boundary condition at the rightmost end of the box. Since ittakes time for the accelerated particles to reach the right boundary,this time is larger than 3 L x / v = 9720Ω − c . tion of the energy density, a further increase of the com-pression ratio is expected because of the contribution ofrelativistic CR fluid (whose adiabatic index is 4/3), butthis effect is still negligible in the simulations presentedhere. We also find the development of a CR-induced pre-cursor, in which the pre-shock fluid is slowed down andheated up because of the pressure in diffusing CRs andbecause of turbulent heating (see Section 6.1 of Caprioli& Spitkovsky 2014a). Note that both the shock precur-sor and the increase of shock compression ratio are alsoobserved in our previous non-relativistic runs. Particle Acceleration
The CR spectrum at the end of the simulation togetherwith its time evolution are shown in Figure 10. Since weexpect the particle momentum spectrum f ( p ) ∝ p − tobe a universal consequence of Fermi acceleration regard-less of being in the non-relativistic or relativistic regimes,we show the time evolution of the particle momentumspectrum in the bottom panel of the Figure. Given ourdefinition of particle kinetic energy E k from Equation(24), and using Equation (35), particle momentum spec-trum is related to energy spectrum by4 πp f ( p ) = Ef ( E ) C p E ( E + C ) . (38)We see that the momentum spectrum, plotted in dimen-sionless form 4 πp f ( p ) /v , exhibits approximately flatshape both in the non-relativistic regime with p/v < .
2, and in the relativistic regime toward higher mo-menta p/v (cid:38)
20, consistent with standard theory ofFermi acceleration. The relativistic part of the spec-trum has slightly smaller normalization than the non-relativistic part, while they connect smoothly throughthe transition. We also show the energy spectrum in themiddle panel of the Figure. Note that for a momentumspectrum of f ( p ) ∝ p − , the energy spectrum shouldnot clearly transition to the f ( E ) ∝ E − scaling untilLorentz factor γ (cid:38)
10 ( ∼ E sh ). Since our CR energyspectrum already drops off at E (cid:38) E sh , we do notyet cover the fully relativistic portion of the spectrum.In the top panel of Figure 10, we see that CRs escap-ing into the far upstream near the end of the run almostcompletely consist of relativistic particles whose energy iswell above 200 E sh . We also find that more contributionscome from non-relativistic particles toward earlier times.These observations are consistent with the reduction of j CR ,x with time discussed earlier and is mainly responsi-ble for the weaker turbulence in the shock upstream.From the simulation, we find that the acceleration ef-ficiency ξ , slowly increases at early time, and saturatesto about 20% in the later half of the run. The relativelylarge acceleration efficiency is in line with the reductionof shock velocity discussed earlier. We also note thatthe value of ξ in this run is about 15% at t ∼ − c ,comparable to 13% from our fiducial run R2. While we also expect the spectral slope to slightly deviatefrom standard value due to the change in shock compression ratio,the deviation is at the level of 0 .
1, which is not distinguishable dueto noise in the particle energy and momentum spectra.
Fig. 9.—
Snapshots of gas density at four different times (as indicated) from our run R2-REL with reduced speed of light C = 10 √ v at four different times as labeled in each panel (unit is Ω − c ). Only a small fraction of the simulation box is shown which is adjusted tocover the vicinity of the expected shock location at x = v sh t = v t/ In Figure 11, we show the evolution of weighted max-imum particle energy. We first note that acceleration isslower than our fiducial non-relativistic run R2 (indicatedin black dash-dotted line). This is most likely due to thereduction in scattering efficiency as particle transitionsinto relativistic regime (e.g., smaller gyro-frequency, andgyro-radius increases with energy more quickly), result-ing in slower acceleration. Since the threshold energyof E = 200 E sh is achieved at very early stage of thesimulation with t < − c , the bulk of the growthin E max occurs in the relativistic regime (note that ourdefinition of E max is a factor of several larger than thecutoff energy). In addition, we find that E max increaseswith time faster at early times ( t (cid:46) − c ), but be-comes slower later. We therefore provide two linear fitsfor t < − c and t > − c respectively, wherethe slopes differ by a factor of ∼
2. The difference mostlikely related to the reduction of CR current as escap-ing particles gradually become dominated by relativisticparticles, which in turn affects the strength of the turbu-lence. While more detailed analysis of particle diffusionis beyond the scope of this work, we have performed fur-ther experiments with larger and smaller particle speedof light, and confirm the effects discussed above. SUMMARY AND DISCUSSION
In this paper, we have rigorously derived an MHD for-mulation which describes the interaction between colli-sionless CRs and a thermal plamsa. Backreaction from CRs to the thermal plasma is incorporated in the formof momentum, energy and electromagnetic feedback, asEquations (19), (20) and (15). While momentum andenergy feedback have been well understood as a resultof momentum and energy conservation in the compositesystem of CRs and gas, we point out the previously ne-glected electromagnetic feedback, which we term as theCR-induced Hall (CR-Hall) term. It becomes importantwhen the relative drift velocity between electrons andions induced by CR streaming approaches the Alfv´en ve-locity, or their ratio Λ defined in (12) approaches or-der unity. Our formulation is applicable (i) on scalesmuch larger than the ion inertial length c/ω pi , which isthe scale that MHD is considered applicable to describethermal plasma, and (ii) when CRs constitute a negligi-ble fraction of the total gas mass (which is almost alwaysthe case in real systems). While our formulation remainsvalid when Λ (cid:38) C can be specified for CR particles so that they canbe either non-relativistic or fully relativistic, as long as C is much larger than any MHD velocity. All ingredi-ents of our MHD-PIC code are well tested with carefullydesigned problems and have achieved excellent perfor-mance. Note that for conventional full PIC or hybrid-HD-PIC Method for CR-Gas Interaction 17 Fig. 10.—
Energy and momentum spectrum of the injectedCR particles from our run R2-REL with reduced speed of light C = 10 √ v . The top panel shows the energy spectrum as a func-tion of x at time t = 11088Ω − c (about the end of the run) in unitsof E sh = v sh /
2. The middle (bottom) panel shows the time evo-lution of the downstream particle energy (momentum) spectrumas marked with different colors indicated by the color table. Theenergy and momentum spectra are shown in dimensionless form Ef ( E ) and 4 πp f ( p ) /v ), respectively. The spectrum is extractedby averaging through a layer at a fixed distance behind the shockfront. As an example, vertical black dashed lines in the top pan-els mark the location of the shock front, and the vertical whitedashed lines mark the layer where the downstream spectrum is ex-tracted. The dash-dotted lines in the middle and bottom panelsmark the transition from non-relativistic to relativistic regimes. Inthe middle panel, the expected energy spectrum from a f ( p ) ∝ p − momentum spectrum is also shown. PIC codes, the requirements to accommodate the largeLarmor radius and diffusion length of the CR particlesand to resolve microphysical plasma scales (on the orderof the ion or electron inertial length) become increasinglyincompatible with higher CR energy. By circumventingthe tiny plasma scales, our MHD-PIC approach enablesus to study the physics of CR-gas interaction on muchlarger, potentially macroscopic, scales with dramaticallyreduced computational cost, while it still retains the fullkinetic nature of the CR particles.Our MHD-PIC code is well suited to investigate a num-ber of astrophysical problems, as discussed in Section 1.In this paper, we have mainly focused on particle accel-eration in non-relativistic shocks, where the Bell insta-bility plays a dominant role in exciting upstream turbu-lence and scattering CR particles. We first point out thatwhen the CR-Hall term is properly taken into account, t ( Ω c-1 )0 2000 4000 6000 8000 10000 12000 w e i gh t ed E m a x / E s h Fig. 11.—
Time evolution of the weighted maximum particle en-ergy ¯ E max from our relativistic run R2-REL. Also shown in dashedlines are linear fits to data at t < − c (red) and t > − c (green). The black dash-dotted line is a reproduction of the fittingresult to run R2 in Figure 6: acceleration is slower in the relativisticcase. the growth rate of the most unstable mode of the Bell in-stability is reduced with increasing Λ, and it is eventuallyovertaken by the filamentation instability upon Λ (cid:38) M A = 30, using an artificialparticle injection prescription with fixed injection effi-ciency η . We observe the development of the Bell insta-bility followed by filamentation and cavitation as the up-stream flow approaches the shock front, leading to strongmagnetic field amplifications and vorticity generation.We find efficient acceleration of the CR particles, whichquickly develop a power-law tail in the energy spectrumwith the expected slope of f ( E ) ∝ E − / when all par-ticles are non-relativistic. The maximum particle energyis found to increase linearly with time, consistent withefficient particle scattering in the Bohm limit. All thesefeatures are consistent with findings from recent hybridsimulations (Caprioli & Spitkovsky 2014a,b,c), justifyingthe ability of our MHD-PIC approach at capturing thenon-linear evolution of CR-driven instabilities and Fermiacceleration in non-relativistic shocks. We also show thatthe CR-Hall term can become important in the vicinityof the shock front in the upstream, especially in the cav-ities, and for high Mach number shocks. On the otherhand, the CR-Hall effect becomes negligible in the shockfar upstream.By choosing a reduced particle speed of light C , we fur-ther perform a shock simulation with substantially largersimulation box and longer duration, and follow particleacceleration from non-relativistic regime to relativisticregimes. We confirm the expected momentum spectrum8 X.-N. Bai et al.of f ( p ) ∝ p − in both regimes, with a small drop in nor-malization in the relativistic part of the spectrum. Wefind that maximum particle energy increases at a slowerrate in the relativistic regime (yet still linear in time),which is likely due to less efficient particle diffusion. Inaddition, escaping CR particles are less effective currentcarriers as they become relativistic, leading to weakerupstream turbulence and slower particle acceleration.A major limitation of our initial study of particle accel-eration in non-relativistic shocks is that particle injectionis handled artificially. We find that with relatively smallinjection efficiency η = 10 − , the shock is less perturbedand particle acceleration efficiency is low ( ξ ∼ η = 4 × − accompaniedby rapid acceleration of particles to ξ (cid:38) η = 2 × − appears to be optimal andgives acceleration efficiency ξ ∼ − η to gradually decreasewith time as particles are accelerated to higher energies,because the energy content of the power-law spectrumwould diverge with increasing maximum particle energy.Therefore, our conclusion of optimal injection efficiency η ≈ × − mainly applies to the initial stage of particleacceleration. Nonetheless, for this effect to be prominent,the CR energy must span several decades in the energyspectrum, which is not easily achieved in current simu-lations.In the future, we will aim to re-design the particle in-jection prescriptions in a way that mimics the particleinjection process in hybrid-PIC simulations. Our cur-rent simple prescription fully relies on the shock detec-tion and tracking algorithm, and its details may have subtle long-term feedback to the shock structure in thevicinity of the shock front. Such dependence will be re-lieved in the new prescriptions. Once properly calibratedwith hybrid-PIC simulations, we will be able to system-atically explore the parameter space of particle acceler-ation in non-relativistic shocks, including Alfv´enic Machnumber and magnetic obliquity, and follow the long-termevolution of the shocks towards larger scales. In partic-ular, the Alfv´enic Mach number in supernova remnantshocks is typically several hundred or higher, which isout of reach of current hybrid-PIC simulations, but willbecome feasible with our MHD-PIC approach.We thank the anonymous referee for a very useful re-port that led to several improvements to the paper, par-ticularly on code implementation, convergence checks,and overall clarity. X.-N.B is supported by NASAthrough Hubble Fellowship grant HST-HF2-51301.001-A awarded by the Space Telescope Science Institute,which is operated by the Association of Universitiesfor Research in Astronomy, Inc., for NASA, under con-tract NAS 5-26555. D.C. is supported by NASA grantNNX14AQ34G. L.S. is supported by NASA through Ein-stein Fellowship grant number PF1-120090 awarded bythe Chandra X-ray Center, which is operated by theSmithsonian Astrophysical Observatory for NASA un-der contract NAS8-03060. This work was partially sup-ported by Simons Foundation (grants 267233 and 291817to AS), and was facilitated by Max Planck/PrincetonCenter for Plasma Physics. Computation for this workwas performed on computational resources supported byPICSciE and the OIT’s High Performance ComputingCenter at Princeton University. This work also usedresources on Stampede at Texas Advanced ComputingCenter through XSEDE grant TG-AST140001. APPENDIX
A: IMPLEMENTATION OF THE MHD-PIC SCHEME: DETAILS “Mass” of CR particles
In our implementation, each CR particle is assigned a “mass” m CR , which in reality represents the CR mass densitycarried by the particle (as a swarm of real particles). The charge and current density carried by particle j in the codeis connected to its “mass” again via the q/mc factor n CR ,j /c = ( q/mc ) j m CR ,j , J CR ,j /c = ( q/mc ) j m CR ,j u j . (A1)The total CR charge and current density n CR and J CR at all grid points can be obtained by particles deposits usingthe TSC scheme.We also define the charge-to-mass ratio ( q/mc ) g for the thermal ions in the gas, which is treated to be independentfrom that for the CR particles. The charge density of the thermal ions n i (or n g ) is then obtained from the gas density ρ via n g /c = ( q/mc ) g · ρ . (A2)The charge density ratio factor R defined in (10) can then be readily obtained based on (A1) and (A2). The Particle Integrator
We briefly describe our implementation of the Boris integrator and how it is coupled with the MHD integrator. Let x ( n ) j and v ( n ) j be the location and velocity of the j th particle at step n . The integration scheme to step n + 1 followsHD-PIC Method for CR-Gas Interaction 19the leap-frog pattern in the form of “drift-kick-rotate-kick-drift”: x ( n +1 / j = x ( n ) j + 12 u ( n ) j ∆ t , v j = v ( n ) j + h c E ( n + ) j , v j = v j + 2 v j + v j × b b × b , v ( n +1) j = v j + h c E ( n + ) j , x ( n +1) j = x ( n +1 / j + 12 u ( n +1) j ∆ t , (A3)where h ≡ ( q/m j c )∆ t , b ≡ h B ( n +1 / j /γ j , and γ j is the Lorentz factor for v j and v j (which are the same). Note therelation between u j and v j defined in (22).Electromagnetic fields are extracted from the half time step in the MHD integrator (with electric field modified bythe CR-Hall term), which we denote as E ( n +1 / and B ( n +1 / . In addition, while cell-centered field satisfies E · B = 0,this is not guaranteed for the interpolated fields, which may cause spurious particle acceleration. Therefore, we alwaysclean the parallel component of the interpolated electric field E → E − E · B /B before it is applied to push the particles.The use of mid-point electromagnetic field (which is fixed during the integration) and the symmetric integrationscheme guarantees second-order accuracy in time. The Boris integrator is time-reversible, and it preserves the geometricproperties of particle gyration in the absence of electric field. Although being straightforward, we do not implementthe Vay’s integrator (Vay 2008). While it better conserves particle energy in the ultra-relativistic regime and whenthe electric field strength approaches the magnetic field strength, it is computationally more expensive. Moreover,in non-relativistic MHD, we always have | c E| (cid:28) | C B | (by construction via the choice of C ). In this regime, the Vaypusher behaves essentially the same as the Boris pusher. The MHD-PIC Integration Scheme
We now describe the way to combine the particle integrator with the MHD integrator in Athena that achieves second-order accuracy in the composite system. For our purpose, the second-order CTU MHD integrator in Athena can beconsidered as a simple predictor-corrector scheme. Our MHD-PIC scheme is built on top of the existing predictor-corrector framework. To simplify the notation, we use f and g to represent particle and gas quantities respectively,and their evolution satisfies the differential equations dfdt = F ( f, g ) , dgdt = G ( f, g ) . (A4)The former describes particle equation of motion under the electromagnetic field provided by the gas, and the latterrepresents MHD equations with CR feedback. Using this notation, our integration scheme can be expressed as follows g ( n +1 / = g ( n ) + G ( f ( n ) , g ( n ) ) ∆ t , (A5a) f ( n +1) = f ( n ) + F ( f ( n +1 / , g ( n +1 / )∆ t , (A5b) g ( n +1) = g ( n ) + G ( f ( n ) + f ( n +1) , g ( n +1 / )∆ t . (A5c)In the above, (A5b) represents the particle integrator described in the previous subsection. Below, we mainly focuson (A5a) and (A5c), which represent the predictor and corrector steps for the gas evolution with CR feedback.The basic procedure to treat backreaction from CRs to the gas involves adding source terms to the gas equationsby depositing relevant physical quantities (momentum, energy, etc.) from the location of individual particles toneighboring grid cell centers. In the predictor step (A5a), the CR feedback is computed using the first equality ofEquations (19) and (20): we deposit the charge and current densities of individual particles from their initial positions x ( n ) j to the grid to obtain n ( n )CR and J ( n )CR , as well as R ( n ) based on the gas density ρ ( n ) . Together with MHD quantities E ( n )0 and B ( n ) , the momentum and energy feedback are directly calculated at cell centers.In the corrector step, we take advantage of the fact that all particles have evolved for a full time step, and computethe momentum and energy difference for each particle j over one time step: d M j = m CR ,j ( v ( n +1) j − v ( n ) j ) .dE k,j = m CR ,j ( E ( n +1) k,j − E ( n ) k,j ) . (A6)0 X.-N. Bai et al.The opposite of these quantities are deposited to the gas, which account for momentum and energy feedback from CRs.The deposit is exerted from particle location at half step x ( n +1 / j . This guarantees the consistency of our algorithmsince the feedback is exerted at the same location as where the particle experiences the Lorentz force. With Athenabeing a Godunov code, exact conservation of total momentum and total energy of the composite gas-CR system isachieved. Moreover, this procedure also ensures our MHD-PIC algorithm to be second-order accurate.Regarding the electromagnetic feedback due to the CR-Hall term, Equation (15) indicates that the correction to theelectric field is directly proportional to the momentum feedback. Therefore, no additional particle deposits are neededand we simply correct the electric field and the energy flux due to the CR-Hall term based on the momentum feedbackin both the predictor and corrector steps. The corrections are implemented in accordance with the intrinsic algorithmthat the CTU integrator handles the induction and energy equations. Timestepping and Code Performance
The timestepping for the MHD part of the code obeys the standard Courant-Friedrichs-Lewy condition. Additionalconstraints on the timestepping from particles include the following.First, particle Larmor time 1 / Ω L = γ ( qB/mc ) − must be resolved. In practice, we define a number Υ so that wedemand Ω L,i ∆ t ≤ Υ for all particles. Our tests in Section B.1 suggest that Υ (cid:46) . .
3, and we have also performed test simulations with Υ = 0 . . .Since CR particles can become relativistic, with velocities much larger than typical fluid velocity, the code can beslowed down substantially. To reconcile this issue, we have implemented sub-cycling for the CR particles: we evolvemultiple steps of particles per MHD step. In other words, the single-step particle evolution illustrated in (A5b) is nowreplaced by multiple steps of particle updates. When evolving the particles, we always use the electromagnetic fieldat half MHD step (i.e., g ( n +1 / ), hence the MHD-PIC scheme remains second-order accurate over a full MHD step.Particles deposit feedback at each sub-step, and the feedback accumulated from all sub-steps will finally be added tothe gas in the corrector step as usual.We have tested the code performance on a single CPU as well as in parallel. For particles, most of the computationtime is spent on the Boris integrator, particle feedback and various interpolations. The serial test takes place on aIntel Xeon E5 2.3GHz processor. Without sub-cycling, it takes about 1 µ s to integrate one particle in 3D, and 0 . µ s tointegrate one particle in 2D. Combined with isothermal MHD using the CTU integrator with third-order reconstructionand the HLLD solver, we find that the code spends roughly equal amount of time on the MHD integrator and on theparticles when using approximately 7 (for 3D) or 6 (for 2D) particles per cell. When run in parallel, because of extrasteps for feedback exchange, the particle numbers quoted above are reduced to 6 (for 3D) and 5 (for 2D). The codeperformance on the particles can be enhanced by using sub-cycling. With the same test on the single processor, wefind that the averaged time to integrate one particle for one particle sub-step is reduced by 30% when using 10 particlesub-steps per MHD step. B: CODE TESTS
The two main ingredients of our MHD-PIC code are the Boris integrator and the particle feedback. We outline heretwo problems that test the two ingredients separately, and show that the code performs well in these respects. Sincein Athena the units for magnetic field are such that magnetic permeability µ = 1, factors of 4 π will be left out here aswell. Gyration Test
We test the Boris integrator with the particle gyration problem. We set up a uniform background gas that movesat constant velocity v in the ˆ y direction. The velocity can be “relativistic” based on the artificial speed of light C ,with Lorentz factor γ = C / (cid:112) C − v . To facilitate analytical calculation, we prescribe other parameters in a framethat is co-moving with the gas. In the co-moving frame, the electric field is zero, and there is a uniform magneticfield B along the ˆ x direction. We include one particle whose parallel and perpendicular four-velocities in the co-moving frame are given by v (cid:107) and v ⊥ , with Lorentz factor γ = (cid:113) C + v (cid:107) + v ⊥ / C . Then we have the Larmor radius R L = v ⊥ ( qB /mc ) − , gyro-frequency Ω L = ( qB /γ mc ), parallel velocity u = v (cid:107) /γ . We compute the particle orbitin the lab frame using the code and compare the results with analytical orbits in the co-moving frame.In Figure 12, we show the test results for non-relativistic particles. We initialize the problem with B = q/mc = v ⊥ = 1, v (cid:107) = 0, C = 10, and hence R L = 1 and Ω L = 1. Considering the fact that the time step in our MHD-PIC codeis generally not a constant, the time step we adopt for this test is ∆ t = 0 . . α , where α is a random number.Two tests are performed, one with zero drift velocity v = 0 (blue dashed), and the other with drift velocity v = 1 (red With the particle module turned on in Athena, the innermosttwo (otherwise just one) ghost zones are updated during the pre- dictor step. The TSC scheme requires that particles can move tothe innermost but not the next ghost zone in the predictor step.
HD-PIC Method for CR-Gas Interaction 21 E ’ z Fig. 12.—
Gyration test for non-relativistic particle orbits with v ⊥ = 0 . C = 1 using variable time step with Ω L ∆ t = 0 . ± .
1. Left: timeevolution of particle kinetic energy in the co-moving frame (indicated by the prime). Right: evolution of particle position in the co-movingframe. Blue dashed lines correspond to zero background flow velocity ( v = 0), red dash-dotted lines correspond to mildly relativistic driftof the background gas ( v = 1). Black solid lines represent analytical results. Note that we have chosen to a large time step to exaggeratethe truncation errors. E ’ z Fig. 13.—
Same as Figure 12, but for relativistic particles with v ⊥ = 10 C . dash-dotted). Shown on the left is the evolution of particle kinetic energy in the co-moving frame, while on the right isthe particle trajectory in the co-moving frame. We see that for v = 0, the integrator conserves energy exactly despitethe variable time step, with relatively small phase error. With drift v = 1, energy is no longer conserved. There issystematic oscillation in E (cid:48) is at 0 .
1% level, which is due to non-zero background electric field. The energy E (cid:48) alsoundergoes a random-walk type drift over time, which is due to the variable time step. Note that both oscillation anddrift in E are inherent to the Boris pusher under these circumstances.In Figure 13, we show the results from the same test but for relativistic particles. We increase particle velocity to v ⊥ = 100 ( γ ≈ R L = 100 and Ω L ≈ .
1, and wemodify the time step to ∆ t = 5 ± cos α so that Ω L ∆ t remains the same as before. We see again that without E × B drift, the integrator captures particle orbits perfectly well. In the presence of drift, the fractional error in the energyis smaller than in the previous case, while the phase error remains similar.Note that we have used relatively large time step (Ω L ∆ t = 0 .
5) to exaggerate the truncation errors. The timestep in real applications is typically much smaller. Moreover, random variation in the timestep also enhances energydrift in the tests, while in real simulations, timestep variation is smooth. Therefore, our particle gyration tests havedemonstrated a reliable performance of the Boris integrator for astrophysical applications.
Linear Growth of the Bell Instability
We test the CR feedback and the CR-Hall term using the linear dispersion relation of the Bell instability, as discussedin Section 4. The test is conducted by evolving the exact eigen-vectors of the most unstable mode.Two sets of tests are performed. In the first test, we suppress the CR-Hall term by using sufficiently small R suchthat Λ ≡ RU s /v A (cid:28)
1. In this case, the most unstable wavenumber is k , and the dispersion relation (29) gives ω = k v A ( (cid:15) + i (cid:112) − (cid:15) ) = k v A exp (cid:20) i (cid:18) π − θ (cid:19)(cid:21) , (B1) The level of energy oscillation increases rapidly with increas-ing gas flow velocity v once v exceeds particle velocity. Never- theless, for all potential applications of our code, we expect CRparticle velocity to be (much) larger than v . ω / k θ ω / k θ ω / k θ Re( ω /k)Im( ω /k) Fig. 14.—
Linear dispersion relation test of the Bell instability in 1D, 2D and 3D with Λ (cid:28) θ , definedas sin θ = (cid:15) , where (cid:15) ≡ v A /u , the ratio of Alfv´en velocity to the CR drift velocity. Symbols denote measured values from the code test,curves represent analytical expectations. See Section B.2 for details. −1 −1 ω / k Λ −1 −1 ω / k Λ −1 −1 ω / k Λ Re( ω /k)Im( ω /k) Fig. 15.—
Linear dispersion relation test of the Bell instability in 1D, 2D and 3D with (cid:15) (cid:28)
R/(cid:15) . Symbols denote measured values from the code test, curves represent analytical expectations. See Section B.2 for details. where we have defined sin θ ≡ (cid:15) . The corresponding eigenvector of this mode reads b y = b ⊥ cos ( k x − (cid:15)k v A t ) ,b z = b ⊥ sin ( k x − (cid:15)k v A t ) ,u y = v A B b ⊥ sin ( k x − θ − (cid:15)k v A t ) ,u z = − v A B b ⊥ cos ( k x − θ − (cid:15)k v A t ) , (B2)where b and u denote transverse perturbations of the magnetic field and fluid velocity, and the wave amplitude growsas b ⊥ ∝ exp (cos θk v A t ).To setup the problem, we consider background gas with uniform density ρ and zero velocity. Let B be thebackground magnetic field, J CR = n CR u be the external uniform current density provided by the CR particles, bothof which are along the ˆ x direction. Note that the CR current is provided by a spatially uniform distribution of CRswith number density n CR that drift at velocity u . For our code test purpose, we treat the CR particles as havinginfinite inertia so that the their effects are simply to provide constant external current and charge.In our code test, we set up the exact eigenvector of the most unstable mode (B2), and test the dispersion relationHD-PIC Method for CR-Gas Interaction 23(B1) as a function of (cid:15) (or θ ). To elaborate, we consider the following parameters: B = ρ = 1 (hence v A = 1).We fit one most unstable wavelength λ into the simulation box, with λ determined by the box size set by the user.This fixes the CR current density: J CR = 2 B k = 4 π/λ . We initiate with a cold population of particles moving atconstant velocity u = v A /(cid:15) along the direction of B . There is one particle per grid cell located cell centers. Themain parameter is (cid:15) ranging from 0 to 1. Correspondingly, the growth rate of the instability is simply √ − (cid:15) k . Inorder for the Lorentz force to be negligible on particle motion over a few e -folding time, we set q/mc ≡ − k for theCR particles. This further determines the value of m CR , which is simply given by m CR = J CR / [( q/mc ) u ] = 2 × (cid:15) .We finally set the gas isothermal sound speed to be c s = 1, and speed of light C = 100 u , both of which are irrelevantfor this problem but needed in the code.The problem generator deals with all 1D, 2D and 3D cases, where in multi-dimensions, the wave vector of the modeis not grid-aligned. The box size in these cases are 1, √ × √ /
2, 3 × . × .
5, resolved by 32, 64 ×
32 and 96 × × λ = 1 with approximately 32 cells perwavelength. We vary (cid:15) from 0 .
01 to 1 .
0, and show the measured real and imaginary parts of ω/k in Figure 14 asa function of θ where sin θ = (cid:15) . Theoretically, the real and imaginary parts should be exact sine and cosine curves,and we see that in all three dimensions the code very well reproduces the analytical phase velocity and growth rate.Although not obvious in the plot, deviations from analytical results in 2D and 3D cases are comparable and largerthan the 1D case, but are well within 2%.In the second test, we relax the constraint on Λ while keeping R (cid:28)
1, which enables the CR-Hall term withoutviolating our formulation. We also keep (cid:15) (cid:28)
1. In this case, the most unstable wavenumber is k m = k /f , where f = 1 + (Λ / , and the dispersion relation reads ω = ( k m v A )(Λ / (cid:112) f ) . (B3)The corresponding eigen-vectors are b y = b ⊥ cos φ ,b z = b ⊥ sin φ ,u y = v A B b ⊥ [(Λ /
2) cos φ + (cid:112) f sin φ ] ,u z = v A B b ⊥ [(Λ /
2) sin φ − (cid:112) f cos φ ] , (B4)where φ ≡ k m [ x − (Λ / v A t ], and the wave amplitude grows as b ⊥ ∝ exp( k m v A √ f t ).We use the same numerical setup as the previous test, except for resetting the eigenvectors. By varying Λ from 0.2to 20, we show the measured real and imaginary parts of ω/k in Figure 15. We see that for all values of Λ, the phasevelocity and growth rate are accurately reproduced from the regime Λ (cid:28) (cid:29)
1. In particular, the measuredgrowth rates deviate from analytical results by less than 0 .
5% for all cases at our numerical resolution. Although ourcode is may not be physically reliable if Λ (cid:38)
C: LINEAR DISPERSION RELATION OF THE BELL INSTABILITY
In this Appendix, we provide more detailed derivation of the linear dispersion relation of the Bell instability takinginto account of the CR-Hall term, following Equation (28) in Section 4. We first define˜ ω ≡ ω − (Λ / v A k , ˜ U s ≡ U s (1 − R/ , (C1)and two factors g ≡ (1 − R )(1 − R/ , f ≡ − R + (Λ / . (C2)Now Equation (28) can be rewritten with the new variablesi(˜ ω − f k v A )( b × ˆ x ) + g kJ CR v A B (cid:18) − ˜ ωk ˜ U s (cid:19) b = 0 . (C3)Non-trivial solution requires b y = ± i b z , corresponding to left/right circularly polarized modes, and it is straightforwardto obtain the dispersion relation ˜ ω − f k v A = ± g kJ CR v A B (cid:18) − ˜ ωk ˜ U s (cid:19) . (C4)To make it more instructive, we further define k ≡ J CR B , (cid:15) ≡ v A ˜ U s , (C5)4 X.-N. Bai et al.and the above dispersion relation can be rewritten as(˜ ω ± (cid:15)k v A ) = f (cid:18) k ± gf k (cid:19) v A − ( k v A ) (cid:18) g f − (cid:15) (cid:19) , (C6)where the two ± signs must be taken to be the same. It reduces to Equation (29) when R (cid:28)
1, or g = 1.= 1.