INTRODUCTION

Cosmic rays (CRs) are relativistic charged particles that can exchange energy and momentum with the surrounding medium, and thus can potentially play a significant role in the kinematic and thermal evolution of astrophysical fluids, as well as being interesting in their own right as probes of high-energy astroparticle physics. CRs are expected to be in equipartition with the turbulent and magnetic energy density in the ISM (Boulares & Cox 1990 ). They have been shown to be influential at galactic and extragalactic scales in a number of recent numerical studies, with implications for galaxy formation (Chan et al. 2019 ;Su et al. 2019 ;Buck et al. 2020 ;Hopkins et al. 2020c ), galactic outflows/winds (Ruszkowski, Ya n g & Re ynolds 2017 ;F arber et al. 2018 ;Hopkins et al. 2020a ), cosmic inflows and virial shocks (Ji et al. 2021b ), and the phase structure of the circumgalactic medium (Salem, Bryan & Corlies 2016 ;Butsky et al. 2020 ;Ji et al. 2020 ;Wang, Ruszkowski & Ya n g 2020 ). In simulations of CRs on interstellar medium (ISM) or star formation or galaxy scales, as well as classical models of Galactic CR transport that are compared to Solar system CR experiments (Strong & Moskalenko 2001 ;J óhannesson et al. 2016 ;Evoli et al. 2017 ), a fluid or Fokker-Planck type approximation for CRs is necessarily adopted, usually with some simplified assumptions for the ef fecti ve 'dif fusion ⋆ E-mail: suoqing@shao.ac.cn (SJ); phopkins@caltech.edu (PH) coefficient' or 'streaming speed' (or some combination thereof). The studies abo v e, and others comparing different assumptions for these CR 'transport coefficients' (Butsky & Quinn 2018 ;Hopkins et al. 2020b ) have shown that the effects of CRs and their interactions with the background plasma, let alone Solar system observables, are highly dependent on these coefficients. Therefore, it is critical to obtain a better understanding of CR transport physics and CR 'feedback' effects (coupling to the gas via magnetic fields) at a microscopic level.

In recent years, hybrid magnetohydrodynamic (MHD)-particlein-cell (PIC) simulations have been developed to investigate the interactions between CRs and background plasma, where CRs are evolved kinetically with the PIC treatment, while the ionized gas is treated as MHD fluid (Bai et al. 2015 ;Mignone et al. 2018 ). This is possible because, for many astrophysical applications including all those abo v e, the background plasma (1) is non-relativistic (bulk speeds u ≪ c ), ( 2) is sufficiently well-ionized that ideal MHD is a good approximation, and (3) has gyro radii, which are vastly smaller than the CRs (e.g. for typical ISM protons and electrons, the gyro radii are factors ∼10 5 -10 6 times smaller than the ∼ GeV CRs carrying most of the CR energy density). Compared with previous pure PIC simulations where both CRs and ions (and sometimes electrons as well) are simulated as particles (e.g. Spitko vsk y 2005 ; Gargat é & Spitko vsk y 2011 ; K unz, Stone & Bai 2014 ), the hybrid MHD-PIC method can simulate CR kinetic effects and feedback to gas by resolving the CR gyro-radius ( ∼au for GeV CRs), without needing to MNRAS 516, 5143-5147 (2022) e xplicitly resolv e the e xtremely tin y background plasma ion inertial length of ∼ 10 -6 au or even smaller background electron skin depth of ∼ 10 -8 au . This method has been pro v en capable of capturing the CR streaming instability (Kulsrud & Pearce 1969 ;Skilling 1971 ;Bai et al. 2019 ) and accounting for the ion-neutral damping effects (Kulsrud & Cesarsky 1971 ;Zweibel & Shull 1982 ;Bambic, Bai & Ostrik er 2021 ;Plotnik ov, Ostrik er & Bai 2021 ). Therefore, the MHD-PIC simulations have a unique advantage in following the evolution of CRs and their feedback/coupling to gas in many astrophysical contexts on much greater scales and over much longer durations than pure PIC simulations.

Even with this approximation, ho we ver, one dif ficulty in simulating CRs is that CRs are relativistic with a velocity close to the speed of light c , far greater than other characteristic MHD velocities (e.g. the sound speed c s , the Alfv én speed v A ), which severely limits the simulation time-step (for any standard Courant-type condition) and thus makes simulations computationally prohibitiv e. F or instance, the speed-of-light c is involved in computing the Lorentz force acting on CRs as follows:

where q is the CR charge, v the CR velocity, B the magnetic field ( E electric field), and c the true speed of light possibly with c ≫ c s , v A , etc. F or conv enience, an artificial-speed-of-light (ASOL) ˜ c is usually used to substitute the true speed of light c in equation ( 1), and the Lorentz factor is accordingly rewritten as Bai et al. 2015 ;Mignone et al. 2018 ;v an Marle, Casse & Marco with 2018 ;Bai 2022 ). With ASOL, as long as large enough scale separations are retained, 'in-code' simulated properities can be properly rescaled to and interpreted as the physical regimes in reality. In other contexts, such as radiation-hydrodynamics (RHD), application of an ASOL is commonly used to allow larger time-steps for the same reason (photons moving 'across the grid' at c require very small numerical time-steps). But in the RHD case, since the modelLing of macroscopic systems allows little flexibility to properly rescale simulated properties, considerable effort has gone into carefully formulating a version of the RHD equations that ensures results are not systematically biased if the ASOL ˜ c ≪ c. In fact, most current RHD formulations (and some applications in other fields; see e.g. Hopkins, Squire & Butsky 2021 ) ensure that steady-state solutions are independent of ˜ c and that the 'true' ( ˜ c = c) values of the photon energy and flux/momentum density can al w ays be reco v ered (see e.g. Kimm & Cen 2013 ;Rosdahl et al. 2013 ;Skinner & Ostriker 2013 , and references therein).

Ho we ver, when dealing with multiscale problems where physical processes with different characteristic lengths/times/energies are involved, simply substituting c with ˜ c in MHD-PIC simulations, hereinafter refereed to as a reduced-speed-of-light (RSOL), might cause physical inconsistencies between simulated and 'real' properties. For instance, since the expressions of the CR energy density, momentum density, and mass/number density contain factors of ˜ c 2 , ˜ c 1 , and ˜ c 0 respectively, simply replacing c with ˜ c for any arbitrary choice of ˜ c = c makes it impossible to simultaneously match the CR energy density, mass/number density, and momentum density to their correct physical values (those with ˜ c = c). In other words, the naive substitution of c in RSOL leads to changes in multiple dimensionless numbers, e.g. the CR-to-gas energy/momentum/mass ratios, ρcr c 2 (Buckingham 1914 ), since characteristic dimensionless numbers in the simulations are changed after the native substitution of ˜ c , simulations with and without RSOL end up with solving totally different systems. The 'back-reaction' forces from CRs to gas thus can be incorrectly calculated due to the CR mass/energy/momentum density mismatch, and the mismatch in RSOL cannot be corrected by 'post-simulation' rescaling once the 'back-reaction' forces have already taken effects on the fly in simulations. Moreo v er, the CR gyroradius also varies with the choices of ˜ c , which leads to inconsistency in related length-scales, e.g. scales of CR scattering, in simulations with RSOL adopted. There is no guarantee, even, with naive substitution of c , that the RSOL universally ensures that the time-step can be increased by a factor ∼ c/ ˜ c (which is the goal of introducing it in the first place), nor that steady-state solutions are independent of ˜ c . These challenges moti v ate us to de velop an accurate formulation of the MHD-PIC equations with a RSOL, which resolves all of the issues abo v e. The paper is organized as follows. Section 2 presents the deri v ation of the ne w RSOL formulation. We compare to previous applications of a RSOL in MHD-PIC simulations in Section 3 , discuss the numerical implementation in Section 4 , and finally conclude in Section 5 .

DERIVATION

Begin by considering the general Vlasov equation for a population of CRs:

where ∇ x , p represent the gradients in position and momentum space, respectively, F is the force, f is the distribution function, the term on the right-hand side represents collisional processes (which we will neglect for now, and return to later), t is time, and this is defined in the lab (simulation) frame. We expect Lorentz forces (equation 1 ) to dominate in the regimes of interest, but note that our deri v ation is robust to the form of F . Following the standard practice for RHD and newly developed methods that evolve the Fokker-Planck equation for a population of CRs (Skinner & Ostriker 2013 ;Hopkins et al. 2021 ), we multiply this equation by 1/ c but then ad hoc replace the value 1/ c only in front of the time deri v ati ve ∂ t f with a different value ˜ c , to introduce the RSOL.

. the time variation of f is systematically slowed by a factor of ˜ c /c, equivalent to a rescaling of time 'as seen by' the CRs. This is what allows us, fundamentally, to take larger time-steps by a factor c/ ˜ c , as the time variation is slower. But this also ensures that one still reco v ers exactly the correct steady-state solutions ( ∂ t f → 0) for the distribution function f , independent of the choice of ˜ c . Now we can turn this into the equation for a 'single' CR group. More accurately for a MHD-PIC method that Monte Carlo samples the dynamics of a CR population, we assume that f is a sum of δfunctions, each representing a group or 'packet' of CRs (which are represented in-code by some 'super-particles' or other tracer field):

where s denotes the species (with individual CR mass m s ), etc., j and N j is the label and total number of individual CR particles of each species s j , W , δ, and f s are the spatial kernel weighting

functions, δ-function, and species function, respectively, and represents MNRAS 516, 5143-5147 (2022) an average of individual CR particles o v er a whole CR 'packet'. Inserting this into equation ( 3 ) and inte grating o v er d 3 x d 3 p d s... times x , p , s ,..., we obtain the evolution equation for each 'packet.'

Ignoring collisions, we have d t N j = 0, d t s j = 0, and d t m s j = 0, where d t ψ j = d ψ j / d t = ψ j represents the Lagrangian time deri v ati ve with each packet of some quantity ψ, i.e. CR number, and individual charge, species, and mass are conserved along the trajectories of individual CRs (since we hav e ne glected spallation and other catastrophic processes).

For x , we obtain: ˜ c -1 ∂ t x j = c -1 v j , or ∂ t x j = ( ˜ c /c) v j . This makes it clear how certain terms should be interpreted. Here, and in all the equations abo v e, v represents the true CR speed, which is an intrinsic property of the individual CR particles. However, ∂ t x j is the effective CR advection speed of CR packets across the grid (which we can defined as v eff j ), which is reduced by our introduction of ˜ c < c.

, so it functions like a RSOL as desired. The momentum equation is (using various properties of δ-functions and cross-products to simplify the algebra):

where E = -( u /c) × B to leading order here (with u the fluid velocity, and ρ the fluid density), and in the last equality we have inserted the assumption that F comes primarily from Lorentz forces, but our expression d t p j = ( ˜ c /c) F j is general and allows for the introduction of any other forces (e.g. gravity).

We can make this more clear by defining the following for notation purposes. Take quantities such as the Lorentz factors β ≡| β| and γ , and likewise the energy E and momentum p of individual CRs to have their true values, so they are intrinsic properties of each CR particle. Then we have:

E j ≡ γ j m s j c 2 (9)

By definition, then the total CR number, mass, momentum, and energy represented by a given 'packet' are: N j , M j ≡ N j m s j , E j ≡ N j E j , and P j ≡ N j p j . We can then rewrite the momentum equation for the CR 'packet' as:

So this is exactly the 'true' acceleration, reduced by ˜ c /c. We now have a momentum equation which can e volve v eff j , the ef fecti v e adv ection speed of a CR 'packet' across the grid, entirely written in terms of background MHD quantities (which retain exactly their usual meaning) and well-defined properties of the 'packet.'

We can derive the 'back-reaction' force on the gas by going back to the original Vlasov equation, calculate the force or change in momentum for ˜ c = c, integrate over the distribution function, and apply the equal-and-opposite change to gas, as in Hopkins et al. ( 2021 ). Or, because we are assuming tight coupling between nonrelativistic ions and magnetic fields, we could also calculate a current effect and use Ampere's law, to get to the force on the gas, as in Zweibel ( 2017 ) and Thomas & Pfrommer ( 2019 ), and obtain the identical answer. The important thing is that just like the radiation case, the'normal' value of c is what appears here, which ensures that the force on the gas/magnetic fields has its correct behaviour (identical to ˜ c = c) when the CR distribution function is in steady state. This gives:

where n j ≡ N j / V j is the local CR number density with V j some ef fecti ve dif ferential volume assigned to the 'packet'.

With these insights, we can also incorporate collisional effects. Returning to equation ( 3 ), ignoring the terms in F , we simply have d t ψ j = ( ˜ c /c) d t ψ j | coll for some ψ and rele v ant collisional effect (e.g. catastrophic losses). So, for example, if catastrophic losses remo v e CRs at a 'true' rate ṅ j = -σv j n n n j (where n n is the local nucleon number density), the total number/mass/momentum/energy of a 'packet' is reduced accordingly as d t ψ j → -( ˜ c /c) σv j n n ψ j for ψ j = ( N j , M j , E j , P j ).

Upon examination, a few points are clear. First, as intended, the manner in which ˜ c enters is equi v alent to a rescaling of time as seen by the CRs, so it guarantees that we can al w ays increase the time-step by c/ ˜ c . This means that the gyro-frequency is reduced in the lab frame (necessarily) by the same factor (as eff g = g ( ˜ c /c)), so the actual effective time it requires to complete a gyro orbit is longer. Ho we ver, this means that for a given true v j or equi v alently CR energy or γ j , the gyro radius on the grid is entirely independent of ˜ c . We can explicitly confirm either if we think about the RSOL as a rescaling in time, or simply examine the 'ef fecti ve' transport speeds across the grid, that the Courant criteria for CR 'advection' (gyro orbits) become e.g. t < C/ eff g = ( c/ ˜ c ) t( ˜ c = c); likewise for advection 'through a cell' (which restricts both CR and MHD time-steps),

Just like with RHD, the 'conserved momentum' here upon MHD-CR exchange is ρ u + ( c/ ˜ c ) P (where P ≡ j P j ), not ρ u + P (which is only true if ˜ c = c). From the form of the Vlasov equation and force on gas written in terms of the distribution function f , it is also immediately obvious that so long as the distribution function is in steady-state , i.e. ∂ t f → 0, then we are guaranteed to reco v er the correct steady-state solutions for the distribution function, i.e. the correct steady-state distribution of CR velocities, phase and pitch angles, momenta/Lorentz factors, etc, as well as the correct force on the gas. We are also guaranteed that processes such as scattering occur on the correct length scales: e.g. scattering rates (or e.g. collisional/catastrophic loss rates) will be slowed by a factor ˜ c /c, but so will transport speeds, so by the time a CR has travelled a given distance ℓ down a magnetic field line, it will have an identical integrated scattering/catastrophic loss probability, independent of ˜ c .

Downloaded

from https://academic.oup.com/mnras/article/516/4/5143/6694104 by California Institute of Technology user on 20 August 2023 MNRAS 516, 5143-5147 (2022)

COMPARISON TO NAIVE RSOL SUBSTITUTION

While 'pure PIC' methods are quite mature, MHD-PIC methods for applications such as CR dynamics have only been widely developed relatively recently. As previously discussed, the ASOL method can be interpreted as changing the unit system while retaining the same dimensionless numbers of a problem. While dealing with multiscale problems where prohibitively small time-step might become an issue, simply treating the RSOL ˜ c , as if the 'true' SOL c were reduced for the CRs specifically, does enable much larger time-steps, but leads to several conceptual difficulties and inconsistencies in an effort to match the physical properties in real systems, including:

(1) it is impossible to simultaneously reproduce the true CR energy and/or momentum and/or number densities, regardless of whatever limit the system is in, since for example in their formulation e j ≡ d E j / d 3 x = n j γ j m s j ˜ c 2 , so either n j or e j must be incorrect (and likewise for momentum density); (2) the transition between relativistic and ultrarelativistic behaviour occurs at the incorrect CR energy/momentum; (3) this leads to a momentum equation (in our notation) d t ( γ j v eff j ) → ( q j / m s j c) ( c E + v eff j × B ), which is different from equation ( 11) by ( c/ ˜ c ) 2 in the electric-field term and ( c/ ˜ c ) 1 in the magnetic-field term, breaking the physical correspondence between the two unless one redefines the electric field E to be different from that assumed in the MHD derivation; (4) likewise if one rescales this momentum equation to keep ( q/m s c) fixed, then this leads not to a rescaling of with ˜ c , but gives the same (requiring the same time-step as simulations with ˜ c = c, for gyro orbits), instead with a smaller gyro radius ˜ c ; (5) there is no guarantee that when the distribution function is in steady-state (even for the simplest homogeneous gas and CR configuration), the distribution function and/or the net back-reaction force on the gas will actually be correct, compared to their value with ˜ c = c. As a result, in these approaches, one obtains certain results that vary systematically with ˜ c and must be extrapolated (e.g. by running a large suite of simulations with different ˜ c ) to estimate their correct ˜ c → c values. In contrast, our approach allows for formal convergence of steady-state properties with ˜ c ≪ c. We finally note that our proposed RSOL formulation actually derives the 'correction terms', which are needed to add to the electric field, CR acceleration equation, and 'back-reaction' force equation, and generically guarantees that the CR distribution function will give the correct steady-state answer. Ho we ver, simply replacing c with RSOL ˜ c , by definition, will not solve the correct CR distribution function in steady state. The previous method is equi v alent to a simultaneous rescaling of both the velocities and spatial units of the original problem, so by decreasing c , it ends up with solving a different problem from what is intended to solve originally, specifically a problem with a different range of spatial scales and a different ratio of v A / c , c s / c and u / c . Whereas, our proposed method, equi v alent to a rescaling of time ('as seen by' the CRs) only, can obtain the right answer for the actually desired dimensionless numbers of the problem simultaneously (and all physical properties as well, e.g. CR mass, momentum, and energy densities), so long as the CR distribution function is in a steady state. Therefore, whenever an RSOL is needed for enabling larger timesteps, our proposed RSOL method should al w ays be implemented in order to achieve the best possible consistency within simulated systems.

NUMERICAL IMPLEMENTATION AND TIME-STEP CONSTRAINTS

The form of the RSOL proposed here can be immediately and trivially implemented in any MHD-PIC code with minimal effort, as it amounts to a straightforward renormalization of the timederi v ati ve terms, which would be computed an yways. We hav e e xplored e xisting implementations in the code GIZMO in Ji, Squire & Hopkins ( 2021a ), but note it immediately applies to any other MHD-PIC implementation, for example both the 'fullf ' and ' δf ' methods in Bai et al. ( 2015Bai et al. ( , 2019 ) ). In these methods, the MHD (non-relativistic ion/neutral/electron) dynamics are solved in some fluid limit with whatev er solv er is most useful, additional physics (e.g. radiation transport or dust as in Ji et al. 2021a ) can be added as well, and then the CRs are integrated on this background with a 'super-particle' approach. In this, one simply identifies each super-particle with a 'packet' abo v e, and immediately obtain their evolution equations. The only difference between our approach and one with no RSOL is the insertion of appropriate factors of ˜ c /c in the code. With this formulation, we see also that when initializing the system, one should simply initialize the correct N j such that e.g. N ≡

N j = M total cr / m s , where m s is the true individual CR mass. Or equi v alently N = E total cr / E j = E total cr / γ j m s j c 2 . This ensures that the correct total CR energy and momentum densities ǫ and P are also initialized -t h e s e are all automatically consistent with one another. Inconsistency only arises if we incorrectly identify γ j m s j v eff j as the individual CR momentum, but this is not correct because the relation between v eff j and the momentum as it appears in conservation and dynamics equations is changed by the introduction of the RSOL. As detailed abo v e, quantities like p , v , γ , β, E , etc. all retain their usual intrinsic physical meaning (that they would have with ˜ c = c), but the actual advection velocity ∂ t x j of a 'superparticle' is v eff j , and the CR momentum equation is modified with the appropriate ˜ c /c factors. The result of this is that all time-step constraints for the CRs also rescale by c/ ˜ c , as all the time-deri v ati ve terms for CRs rescale with this. Obviously, one requires t

), where r g j is the gyro radius and C is a Courant factor, or more generally t

, and if there are some catastrophic losses, similarly t j < C N j / | Ṅ j | ∼ ( c/ ˜ c ) / ( σv j n n ). And for the Courant condition of CRs propagating through the gas with some cell size x , we have t

For the case with non-relativistic CRs, even though it does not gain significant benefit from the RSOL regarding computational time saving, this RSOL implementation is still recommended: it naturally returns to the ASOL formulation under such regimes, and can be safely and conv eniently e xtrapolated to relativistic regimes where an RSOL might be needed.

Finally, just like any other RSOL implementation, convergence to correct solutions requires choosing a sufficiently large value of ˜ c compared to other (e.g. MHD) signal speeds in the problem. Although formally steady-state solutions are independent of ˜ c for all ˜ c with this method, if one chose ˜ c too small compared to e.g. some rms turbulent velocity of the MHD fluid, ˜ c ≪| δu | , then the plasma would artificially 'outpace' the CRs and the system will not actually be able to converge to a steady-state distribution function f in the first place. This is also consistent with the requirement of sufficient scale separation in order to ensure, e.g. a much smaller gyro-period than the eddy turno v er/Alfv én crossing time-scales.

CONCLUDING REMARKS

In this paper, we propose an accurate, well-defined formulation of an RSOL that can be applied to MHD-PIC simulations. The major advantages of the proposed formulation are that: (1) the time-step can be safely increased by c/ ˜ c as intended (regardless of if the time-step restriction from gyro orbits, adv ection o v er the grid, grid response to CRs, or collisional terms dominates). Meanwhile, (2) the key properties of CRs in numerical simulations, such as energy density, momentum density, mass/number density, and CR gyroradius, are independent of the choices of ˜ c , and (3) they can simultaneously reproduce their true, physical values (those one would have with ˜ c = c). And perhaps most importantly, (4) solutions for CR bulk properties, distribution function, and their interactions with/forces on the magnetic fields and background fluid are independent of ˜ c as well (even if ˜ c ≪ c) so long as the CR distribution function f approaches a local steady-state. These advantages naturally recommend the formulation here as a standard and straightforward implementation for MHD-PIC simulations whenever an RSOL is required.