Spectrally resolved cosmic ray hydrodynamics -- I. Spectral scheme
Philipp Girichidis Christoph Pfrommer Michal Hanasz Thorsten Naab



firstpage–lastpage
Probing muonic forces with neutron star binaries
    Toby Opferkuch
    December 30, 2023
================================================





Cosmic ray (CR) protons are an important component in many astrophysical systems. Processes like CR injection, cooling, adiabatic changes as well as active CR transport through the medium strongly modify the CR momentum distribution and have to be taken into account in hydrodynamical simulations. We present an efficient novel numerical scheme to accurately compute the evolution of the particle distribution function by solving the Fokker-Planck equation with a low number of spectral bins (10 - 20), which is required to include a full spectrum for every computational fluid element. The distribution function is represented by piecewise power laws and is not forced to be continuous, which enables an optimal representation of the spectrum. The Fokker-Planck equation is solved with a two-moment approach evolving the CR number and energy density. The low numerical diffusion of the scheme reduces the numerical errors by orders of magnitude in comparison to classical schemes with piecewise constant spectral representations. With this method not only the spectral evolution of CRs can be computed accurately in magnetohydrodynamic simulations but also their dynamical impact as well as CR ionisation. This allows for more accurate models for astrophysical plasmas, like the interstellar medium, and direct comparisons with observations.



cosmic rays – methods: numerical – MHD – plasmas – astroparticle physics










§ INTRODUCTION


CRs are an important energy component in many astrophysical systems from proto-planetary discs, to the interstellar medium in galaxies, and to galaxy clusters <cit.>. Of particular importance are CR protons in galaxies because their energy densities are comparable to magnetic, thermal and kinetic energy densities and because of their resulting dynamical and chemical impact on the gas. Early theoretical models  <cit.> highlight the acceleration of CRs at strong shocks via diffusive shock acceleration, which has been successfully modelled numerically <cit.>. For Galactic CRs the most abundant shocks are supernova remnants with evidence of hadronic particle acceleration <cit.>, see also <cit.> and <cit.>.

The coupling of CRs to the gas is mediated via the generation of Alfvén waves and the resulting scattering. The simplest way to describe CRs as a fluid is the one-moment approach in the scattering angle μ under the assumption that the particle distribution function is isotropic <cit.>. A further simplification for numerical models is a grey approach, in which the distribution function times the kinetic energy per particle is integrated over momentum space and the resulting CR energy density is evolved in time and space. Recently, <cit.> and <cit.> extended the traditional one-moment schemes in the scattering angle to two-moment descriptions that include a more self-consistent coupling of CRs to the plasma and captures streaming and diffusion relative to the gas rest frame.

Previous studies have modelled the dynamical impact of CR protons in the galactic ISM by including them as a relativistic fluid with an effective adiabatic index, see e.g. <cit.>. The first dynamical coupling was performed by <cit.> studying the Parker instability. A prominent application of CR hydrodynamics is their role in driving galactic outflows <cit.>, but see also applications in galaxy clusters <cit.>. Previous hydrodynamical simulations used global disc setups <cit.> with a focus on the large-scale dynamics as well as stratified boxes of a representative fractions of the ISM <cit.> with a focus on the chemical evolution, the detailed CR coupling and the relative importance between CRs and SNe as a driver.

All of the previous studies reveal that CRs can provide relevant pressures and accelerations of the gas. However, to what extent they provide a main contribution to the dynamical evolution depends on the system under consideration and the details of the CR parameters. In particular the combination of CR losses and their spatial transport might be important in determining the global impact. Both processes are strong functions of the spectral energy distribution of CRs, which is not resolved in current CR-MHD simulations but only integrated to yield a total CR energy <cit.>.

The dynamical work of CRs on the gas – and vice versa – results in adiabatic changes that are connected with compression and expansion of the gas. In addition low-energy CRs lose energy via Coulomb collisions with the thermal particles of the gas. Collisions of high-energy CRs at energies above GeV with the gas result in catastrophic hadronic losses via the production of neutral pions and their decay into γ-ray photons. Strong shocks further accelerate CRs. All together the CR energy and the spectral distribution is constantly changing.

In the interstellar medium the efficient coupling of CRs with the gas via magnetic fields results in a non-negligible effective CR pressure that thickens the galactic disc and launches outflows from the galaxies <cit.>. This direct dynamical impact is mainly due to CRs with momenta of a few GeV/c. Low-energy CRs suffer from strong Coulomb losses, which reduces their energy density and results in a negligible impact via their pressure. However, a perceptible increase of the cross section of MeV-to-GeV CRs with the thermal gas causes efficient ionisation of the gas <cit.>. As CRs can penetrate deeply into dense molecular clouds, they influence the formation of stars and the observational signatures.

Above a total energy of E_thr=1.22 GeV CRs are energetic enough to produce pions, which in turn decay into γ-ray photons, secondary electrons and neutrinos. The secondary electrons can emit radio synchrotron emission in ubiquitous magnetic fields and Compton upscatter ambient radiation fields into the X-ray to gamma-ray regime. A predictive modelling of the resulting non-thermal emission processes calls for self-consistent spectral modelling of the CR spectrum in time and space. While the CR spectrum, its composition, and the non-thermal radiative signatures are modelled in CR propagation codes <cit.>, such approaches assume the Galaxy to be static and adopt observationally inferred distributions of the gas density, magnetic fields, and CR sources that are not necessarily emerging from a self-consistent simulation of a dynamically evolving galaxy.

As CRs have a dynamical impact it is thus favourable to follow the spectrum together with the hydrodynamical evolution, i.e. to compute a full spectrum for every computational cell. This poses strong constraints on the numerical scheme. We would like to follow a large dynamical range in CR energy from below MeV to above TeV. The CR spectrum itself is very steep, i.e. covers a large dynamical range in amplitude. Nonetheless we can only represent the spectrum with a low number of spectral bins (∼10-30), which requires a relatively complex numerical scheme compared to standard methods with orders of magnitude larger spectral resolution.

In a series of papers we introduce a novel implementation of the spectral CR energy distribution, which allows a dynamical coupling of CRs with the gas as well as an accurate evolution of the CR spectrum for the relevant mechanism by only using a low number of spectral bins compared to classical spectral approaches. In Section <ref> we present the theoretical background and present analytical solution to idealised cases. In Section <ref> we outline the spectral discretisation of the particle distribution function. In Section <ref> we derive a numerical scheme for the time evolution of the CR spectrum and show one-zone tests of it in Section <ref>. We present a one-dimensional test of energy dependent spatial diffusion in Section <ref> and conclude in Section <ref>.



§ EVOLUTION OF CRS



Before presenting the discretisation scheme and the numerical algorithms, we review the time evolution of CRs including the individual loss processes as well as the combined solutions of freely cooling and the steady state spectrum for continuous injection, which we will use to compare and scrutinise our numerical simulations.



 §.§ Theoretical background


Cosmic rays are charged particles and therefore interact with the magnetic field. Quasi-linear theory in the frequent scattering limit leads to the Fokker-Planck equation for the phase space particle distribution function <cit.>. Throughout the paper we use the three-dimensional form of the distribution function, f=f^(3). Other studies <cit.> use the equivalent one-dimensional form f^(1), where f^(1) = 4π p^2 f. The Fokker-Planck equation then reads

    ∂ f/∂ t =    -u f_advection + D_xx f_diffusion + 1/3up∂ f/∂ p_adiabatic process
        + 1/p^2∂/∂ pp^2b_l f + D_pp∂ f/∂ p_other losses and Fermi II acceleration + j_sources,

where f = f(x,p,t) is the isotropic part of the distribution function and D_xx and D_pp are the spatial diffusion tensor and the momentum space diffusion coefficient, respectively. The losses are described as b_l(x,p,t)= p/ t. We note that formally only continuous losses can be cast into this form but not impulsive losses like, e.g. hadronic losses because they do not conserve the number of particles <cit.>. However, we will be taking moments of equation (<ref>), which implies integrating over the distribution function  that consists of a large population of individual particles. It is therefore justified to take the continuous limit of the collection of
these interactions. The term j(x,p,t) describes CR sources. Here, we neglect CR streaming for simplicity <cit.>. Because second-order Fermi acceleration is a slow process that only mildly shifts the CR spectrum toward larger momenta, we postpone a treatment of turbulent reacceleration to future work and set D_pp=0. For the remainder of the paper we therefore only discuss spatial diffusion (D_xx) and omit the subscript xx. Generally, the distribution function and all its components depend on the position in space, x and time, t. In order to simplify the equations, we omit these dependencies unless explicitly needed.

The number density n and energy density e are then given by the appropriate moments of f,

    n_CR   = ∫_0^∞ 4π p^2 f(p)   p
    
      e_CR   = ∫_0^∞ 4π p^2 f(p) T(p)  p.

Here, T(p) is the kinetic energy of the CRs,

    T(p) =   √(p^2c^2 + m_p^2c^4) - m_p c^2,

where m_p is the proton mass and c is the speed of light. The CR pressure is given by

    P_CR   = ∫_0^∞4π/3  c p^3 β(p) f(p)  p
       = ∫_0^∞4π/3  f(p) p^4c^2/√(m^2c^4+p^2c^2) p

where β(p) = p/√(p^2+(m_pc)^2).



 §.§ Adiabatic changes


We start with the simplest process, which is the adiabatic one,

    ∂ f/∂ t_ad = -1/3u ∂ f/∂ln p.

The divergence of the velocity field u is constant over a spectral time integration because we apply operator splitting for the hydrodynamics and the spectral evolution. We note that the adiabatic process is simply equivalent to an advection in logarithmic space with the advection speed -1/3 u. This means that the local slope in a comoving frame u does not change. Using

    ∂ f/∂ t = ∂ f/∂ p  p/ t

in equation (<ref>) we find 

    p/ t = - 1/3u p.

Separation of variables and integrating from p (t_0) to p(t), which corresponds to the temporal evolution from t_0 to t=t_0+Δ t, yields

    p(t) = p(t_0)  exp-∫_t_0^t_0+Δ t1/3u  t.




 §.§ Coulomb losses






The total energy loss per proton is <cit.>

    -[ T(p)/ t]_C   = ω_pl^2 e^2/β c[ ln( 2m_e c^2β p/ħω_pl m_pc) - β^2/2],

where ω_pl=√(4π e^2 n_e/m_e ) is the plasma frequency and n_e is the electron number density. Coulomb and hadronic losses are shown as a function of momentum for different electron densities in Fig. <ref>. We note that the weak scaling with n_e in the logarithmic term can be neglected in comparison to the linear scaling of the Coulomb loss term (∝ω_pl). To illustrate the scaling we can express the losses as a function of p instead of T and can use the simplified approximation

    b_C   ≡[ p/ t]_C
       ≈  -10^-18 erg cm^3 s^-1n_e/c[1+(p/GeV/c)^-1.9],

which is accurate to 17% over the range shown. The scaling with p^-1.9 for low momenta leads to a finite cooling time, at which the momentum reaches zero. The cooling times to p=0 are plotted as a function of initial momentum in Fig. <ref> for different electron densities. Typical hydrodynamical time steps (Δ t_hydro) are indicated by the grey area, which illustrates that for momenta below 1 GeV/c the hydrodynamical time step might be larger than the cooling time. For t_cool≲Δ t_hydro we can use the numerical cooling scheme. If cooling occurs on much shorter time-scales compared to the simulation time step we can directly evaluate the steady state solution without intensive numerical integration. For the tests in this paper we can set the time steps independently of any hydrodynamical simulation and therefore use the fully numerical solution.



 §.§ Hadronic losses



Hadronic losses occur as a result of inelastic reactions of CRs with the gas. The interactions mainly result in the production of pions if the CR energy exceeds the threshold momentum, p_thr≈0.78 GeV/c. The loss rate of kinetic energy is given by

    - T/ t   = c n_Nσ_ppK_pT(p)θ(p-p_thr).

Here, n_N is the target nucleon density of the ISM, σ_pp is the total pion cross section and K_p≈ 1/2 is the inelasticity of the reaction <cit.>. The CR momentum losses are shown in Fig. <ref>. The cooling times for hadronic losses are shown in Fig. <ref>. We note that the typical cooling times for hadronic losses are long compared to typical hydrodynamical time steps. We note that the losses scale with the kinetic energy, T(p), which asymptotes to a linear scaling in momentum, p, for relativistic CRs. As a consequence, the spectral changes due to hadronic losses in the relativistic case are

    ∂ f/∂ t_hadr   ∝1/p^2 ∂/∂ p(p^2 b_hadrf),

where

    b_hadr =  p/ T ( T/ t)_hadr.

If f is a power-law in p, we note that the hadronic losses have the same scaling with momentum as f,

    ∂ f/∂ t_hadr   ∝ f,

which leaves the slope of the particle distribution function unchanged.



 §.§ Fermi-I acceleration


Fermi-I or diffusive shock acceleration is an important gain process of energy, in which thermal particles can be accelerated to super-thermal energies by means of adiabatic compression and expansion together with spatial diffusion at shocks <cit.>. In principle both processes are included in the Fokker-Planck equation, so that Fermi-I acceleration can be accounted for analytically. However, a numerical solution requires resolutions of the kinetic physics, in particular the scattering of particles due to the non-resonant hybrid instability on small spatial scales around the shock front <cit.>, which is impossible to resolve in our applications of astrophysical fluid dynamics. We therefore need to treat the acceleration of CRs as a subgrid model and effectively describe it as an injection of CR energy in regions of strong shocks. Given the shock compression ratio r=ρ_post/ρ_pre, with the pre and post-shock densities ρ_pre and ρ_post, we expect injection with a spectral index

    _acc = 3r/r-1.

The effective Fermi-I acceleration is thus encoded in the sources j with

    j_acc = A_acc p^-_accexp(-p/p_acc),

where the acceleration efficiency and thus A_acc and the maximum momentum p_acc depends on the local shock conditions <cit.>, see also <cit.>.



 §.§ Free cooling





  §.§.§ General theoretical considerations for the spectral behaviour


We can investigate the expected spectral changes by assuming conservation of the total CR number density in the same volume,

    n(t_1)    = n(t_0),
    
    
    4π∫_p_1^∞p^2 f(p,t_1) p    = 4π∫_p_0^∞p'^2 f(p',t_0) p'.

Here, p_0 and p_1 are related via the change in momentum over time, which we address below. Assuming that f(∞) vanishes, i.e. the number density of CRs is finite, we can assume conservation of the number density in an infinitesimal momentum interval to find

    f(p_1,t_1)    = f(p_0,t_0) p_0^2/p_1^2  p_0/ p_1,

where the differential changes in p, p, depend on the momentum and thus on the time, such that we can write

    p_0/ p_1   =  p/ t(p_0)  p/ t(p_1)^-1.

The changes in momentum are described by the generalised loss term

    p/ t_X   = b_X.

where we assume a simple power-law scaling with momentum p

    p/ t_X   = b_X,0 p^h,   h≠-1   for simplicity.

We solve the differential equation by separation of variables and integration

    ∫_p_0^p_1  p'/b_X(p)   = (t_1-t_0)

and we can solve for the momentum p_1

    p_1 = p_0^1-h + (1-h)b_X,0(t_1-t_0)^1/(1-h).

For losses b_X,0 < 0, p_1<p_0 and depending on h the momentum will cool to p_1=0 within a finite time. Inserting the momentum as well as the loss rates into equation (<ref>) yields

    f(p_1,t_1)    = f(p_0,t_0) p_0^2/p_1^2 p_0^h/p_1^h.




  §.§.§ Application to approximate Coulomb cooling


Coulomb cooling for low momenta scales as p^-1.9, i.e. h=-1.9. Hence, we find

    f(p_1,t_1)    = f(p_0,t_0) p_0^0.1p_1^-0.1
       =f(p_0,t_0) p_0^0.1p_0^2.9 + 2.9 b_C,0(t_1-t_0)^-0.1/2.9.

As long as we do not cool the momentum to zero we can assume that p_0^2.9 > 2.9 b_C,0(t_1-t_0). For p_0^2.9≫ 2.9 b_C,0(t_1-t_0) we can neglect the second term in the sum and find again

    f(p_1,t_1)    = f(p_0,t_0),

a vanishing scaling with p, i.e. a flat slope.



 §.§ Steady state spectrum



Here we focus on the steady state solution resulting from cooling and continuous CR injection. The left-hand side of Equation (<ref>) thus vanishes and the resulting equation reads

    ∂ f/∂ t    = 0 = 1/p^2∂/∂ pp^2 b_l(p) f(p) + j(p),

which can be rewritten to yield

    f(p)    = -1/p^2 b_l(p)∫_p^∞p^2 j(p)  p.

The losses include hadronic and Coulomb interactions

    b_l(p)     =  p/ t_tot =  p/ t_hadr+ p/ t_Coul.

The injection j(p) is modelled as a power-law spectrum with

    j(p)    = Ap^-.

The steady state solution converges for >3 and the steady state solution reads

    f(p)    = A p^-+1/(-3) b_l(p)




§ SPECTRAL DISCRETISATION OF THE PARTICLE DISTRIBUTION FUNCTION








 §.§ Discretization in momentum


In principle we can evolve the Fokker-Planck equation using f directly, similar to previous approaches as in GALPROP <cit.>, PICARD <cit.>, DRAGON2 <cit.> or CREST <cit.>. However, a standard discretisation with piecewise constant values for f requires relatively high spectral resolution in order to obtain accurate results <cit.>. Reasonably low errors are achieved with approximately 50-100 bins per momentum decade, i.e. much more than a few hundred bins for a full spectrum ranging from p∼10^-1 to 100 GeV/c. This high number stems from the fact that we need to cover a large range in momentum and  – due to the steep CR spectra – a large dynamical range of f.

For our CR module we focus on CR protons and their dynamical impact on the hydrodynamical evolution. This requires to solve the spectrum in every hydrodynamical cell and evolving hundreds of momentum bins is not feasible. Instead we aim for a relatively low number of spectral bins of order 10-20. As a consequence, we need a numerical scheme, which can accurately treat the large dynamic range with only several cells and does not suffer from strong numerical diffusion. We therefore chose a logarithmic spacing for the spectral discretisation and describe the particle distribution function as piecewise power laws.

The spectral distribution in principle covers the entire momentum space. We restrict our computation to a finite range p_min<p<p_max. We discretise the spectral distribution with N_bins bins between p_min and p_max as illustrated in Fig. <ref>. The logarithmically spaced bins include two buffer or boundary bins at the lower and the upper end of the spectrum, where the lower buffer bin ranges from p_min to p_- and the upper one from p_+ to p_max. The bins between p_- and p_+ are equally spaced in log p. We define cell centred quantities with index i and corresponding left-hand cell-faced quantities with index i-1/2. We adopt a distribution function f as a piecewise power-law

    f(p) = f_i-1/2p/p_i-1/2^-_i,

with the slopes _i (blue lines in Fig. <ref>). We note that this functional approximation of f has two degrees of freedom per bin, namely the amplitude f_i-1/2 and the slope _i. We therefore also investigate two moments of f per bin, i.e. number and energy density, which are given by

    n_i    = ∫_p_i-1/2^p_i+1/2 4π p^2 f(p)   p
    
      e_i    = ∫_p_i-1/2^p_i+1/2 4π p^2 f(p) T(p)  p.




 §.§ Motivation for a two-moment approach




Several methods have been proposed to solve the Fokker-Planck equation using a piecewise power law representation <cit.>. If the numerical scheme is only based on one moment the two degrees of freedom per bin need a closure relation. A simple and intuitive assumption is to require a continuous function for f. This reduces the degrees of freedom for N bins to N+1, i.e. there is only one additional condition that needs to be set. <cit.> force the slopes in the first two bins to be equal. <cit.> assumes the proton spectrum to be of constant curvature, i.e. _i+1-_i = _i -_i-1. However, there is a fundamental problem with a continuous description of f in particular for a locally varying spectrum. Let us assume a steady state spectrum as shown in the blue curve in Fig <ref>. This spectrum can nicely be represented by a continuous function of piecewise powerlaws. If we now inject energy at the three lowest bins and force the spectrum to be continuous after the injection (dark red curve), we force changes of the local slope across the entire spectrum. The final continuous representation then alternates between a concave and a convex spectrum. Avoiding this alternating behaviour would require to reshape also the high-energy part of the spectrum, which should effectively be unchanged if only low-energy CRs are injected. Without the restriction of a continuous distribution function we can still model the spectrum with a physically useful description (black lines). Even if for most physical applications (injection, cooling, diffusion) there are possibilities to keep a continuous spectrum, the discontinuous representation allows for more freedom and a more stable numerical treatment. The discontinuous modelling, however, requires to constrain two degrees of freedom, which we simply chose to be the two moments of the particle distribution function.



 §.§ Moments of the distribution and their time evolution


Instead of evolving the momentum-integrated systems in time we need to solve for changes in each momentum bin separately. One special case is advection with the gas, which does not involve spectral changes, but simply advects the spectrum across all bins. We translate the time evolution of the Fokker-Planck equation to the evolution equation for the CR number and energy density in bin i. For clarity, we omit the subscript CR in the following two equations. The time evolution of n_i reads

    ∂ n_i/∂ t   = ∂ n_adv,i/∂ t +∂ n_diff,i/∂ t + ∂ n_ad,i/∂ t + ∂ n_l,i/∂ t + j_n,i
       = -un_i + D_n n_i
       +1/3up + b_l(p) 4π p^2 f_p_i-1/2^p_i+1/2 + j_n,i,

where the individual terms describe advection (subscript adv), diffusion (diff), adiabatic changes (ad), and losses (l). Sources are indicated by j.
The energy evolution is given by

    ∂ e_i/∂ t   = ∂ e_adv,i/∂ t +∂ e_diff,i/∂ t + ∂ e_ad,i/∂ t + ∂ e_l,i/∂ t + j_e,i
       = -ue_cr,i + D_e e_cr,i
           + 4π/3u(T(p)p^3f_p_i-1/2^p_i+1/2.
                   .-∫_p_i-1/2^p_i+1/2 fp^4c^2/√(m^2c^4+p^2c^2) p)
           + 4πTp^2b_l(p)f_p_i-1/2^p_i+1/2
           -4π∫_p_i-1/2^p_i+1/2p^3c^2b_l(p)f/√(m^2c^4+p^2c^2) p + j_e,i

In both equations (for n_i and e_i) we rewrote the diffusion term, such that it formally takes the form of a simple diffusion equation with modified diffusion tensors D_n and D_e,

    ∂ n_diff,i/∂ t   =∫_p_i-1/2^p_i+1/2 4π(D f)p^2  p
    
         = D_n n_i

with

    D_n =  n_i^-1∫_p_i-1/2^p_i+1/2 4πD  f p^2  p

and analogously for the energy equation,

    D_e =  e_i^-1∫_p_i-1/2^p_i+1/2 4πD  f p^2T(p) p.

We note that the the inverse of the gradient needs to be computed for every individual component. We discuss the effective diffusion coefficient in more detail in Section <ref>.



§ DERIVATION OF A NUMERICAL SCHEME FOR THE TIME EVOLUTION





 §.§ Computing changes in number and energy density






We use operator splitting for the individual parts of the time evolution of the Fokker-Planck equation, in particular the evolution in space and the spectral evolution. In this paper we focus on the spectral evolution and discuss the integration of the method into hydrodynamics in a subsequent paper. We describe the individual parts for the physical processes in terms of the discretised momentum bins with piecewise power laws following <cit.>. In general we convert the evolution in time into an evolution in momentum,

    p/ t = F(p, t, T(p), u),

where F is a function that depends on the individual physical processes. In particular we cover:


  * adiabatic gains and losses,

  * injection via diffusive shock acceleration,

  * hadronic losses and Coulomb losses, and

  * momentum-dependent diffusion.


For any given physical process and a given integration time step Δ t we can compute the change in momentum. Without loss of generality, we would like to illustrate this for a loss term, b_l,

    p/ t = b_l < 0,

which we can rewrite as

    ∫_p_loss^p_x p/b_l = ∫_t_0^t_0+Δ t t = Δ t

if the loss process does not explicitly depend on time. The initial momentum p_loss cools during Δ t to p_x. For the adiabatic process Eqs. (<ref>) and (<ref>) explicitly show the closed formulation. We chose p_x to be the momentum at the bin boundary, p_x = p_i+1/2, and compute the corresponding momentum p_loss, which is illustrated in Fig. <ref>.

At this point we would like to discuss the constraints on the maximum integration time step. As in other explicit numerical methods that account for only the immediate neighbour cells, the flux of number and energy density can extend to at most one bin. For the gain and loss momentum in our setup this means that p_i+1/2≤ p_loss < p_i+3/2 and p_i-1/2< p_gain≤ p_i+1/2. Depending on the loss and gain process, the maximum time step is implicitly given via equation (<ref>) by limiting the maximum fraction of the bin that p_gain or p_loss should occupy. We find satisfactory results for

    p_loss/p_i+1/2   ≤ 0.4(p_i+3/2/p_i+1/2)  and
    p_i+1/2/p_gain   ≤ 0.4(p_i+1/2/p_i-1/2).


 
The shift of momentum from p_loss to p_i+1/2 corresponds to a transport of particle and energy density across the spectral bin boundary i+1/2,

    Δ n_i+1/2   = 4π∫_p_s,gain^p_i+1/2  p^2 f(p)  p,
    Δ e'_i+1/2   = 4π∫_p_s,gain^p_i+1/2  p^2 f(p) T(p)  p,

 illustrated by the red area in Fig. <ref>. This is the essential part of the time evolution of the new method. We evolve the number and energy density in bin i+1 as

    n(t+Δ t)    = n(t) - Δ n_i+1/2 + Δ n_i+3/2,
    
    e'(t+Δ t)    = e(t) - Δ e'_i+1/2 + Δ e'_i+3/2.

The transport of Δ n from bin i+1 to bin i does not require any further correction. Contrary, for the energy term we need to take into account that the loss in momentum connects to a loss in energy, i.e. shifting the spectrum towards lower momenta, we need to additionally correct for that shift,

    Δ e_i+1/2 = Δ e'_i+1/2 T(p_i+1/2)/T(p_loss,1).

The transfer into and out of cell i+1 is illustrated in Fig. <ref>. We can compute the energy correction factors independently for Δ e at interface p_i+1/2 and p_i-1/2. In case of identical ratios of p_loss,1/p_i+1/2 and p_loss,2/p_1+3/2 the energy correction factors for Δ e'_i+1/2 and Δ e'_i+3/2 are similar. In the non-relativistic limit as well as in the relativistic limit, in which the particle energy has a simple scaling with the momentum (∝ p^2 and ∝ p, respectively), the same ratio for p_loss,1/p_i+1/2 and p_loss,2/p_1+3/2 results in the same correction factor for the transfer into and out of cell i+1, i.e.

    e(t+Δ t) = e'(t+Δ t) T(p_i+1/2)/T(p_loss,1).

In the transition region between the classical and the relativistic energy, the factors differ. If the correction factors at the left and right boundary of the bin are not largely different, we can simply apply Δ e'_i+1/2 and Δ e'_i+3/2 and then correct the entire modified bin i+1 by an arithmetic average of the correction factors,

    e(t+Δ t) = e'(t+Δ t) 1/2(T(p_i+1/2)/T(p_loss,1) + T(p_i+3/2)/T(p_loss,2)).


Depending on the process the ratios p_loss,1/p_i+1/2 and p_loss,2/p_i+3/2 behave differently. For adiabatic expansion (and compression) both ratios are the same as well as for the hadronic losses as long as all momenta are above the threshold momentum for hadronic losses. In these cases, the energy correction factor for the bin is trivial. The case of Coulomb losses at sub-relativistic energies is more difficult to compute. As the loss rate scales with approximately p^-1.9, the loss rate at the lower momentum boundary can be significantly larger than at the high momentum boundary, in particular, if the bins span half an order of magnitude. In this regime we therefore compute the correction term numerically using a subgrid interpolation. We subdivide the bin in 10 logarithmically spaced sub-bins with momenta p_j(t_0) and compute the momenta after Δ t, p_j'(t_0=Δ t), using Eqs. (<ref>) and (<ref>). We then average over the individual values ⟨ p_j'(t_0+Δ t)⟩_j.

In principle, we can compute a total shift for each bin for all terms in the Fokker-Planck equation. However, it is simpler to compute a separate shift momentum for the individual processes. The adiabatic process can also be used to analytically show that the method of discretisation yields correct results, which we demonstrate in appendix <ref>.



 §.§ Reconstruction of the particle distribution function



After computing the temporal changes we end up with a modified number and energy density in each bin. We now need to reconstruct the new amplitude and slope of the particle distribution function in every bin. We first compute the new slope by solving the ratio

    e_i/n_i   = ∫_p_i-1/2^p_i+1/24π p^2 f(p) T(p)  p/∫_p_i-1/2^p_i+1/24π p^2 f(p)  p
       = ∫_p_i-1/2^p_i+1/24π p^2 f_i-1/2 (p/p_i-1/2)^-_i  T(p)  p/∫_p_i-1/2^p_i+1/24π p^2 f_i-1/2 (p/p_i-1/2)^-_i  p

for the slope _i numerically using the Newton-Raphson method. The solution of _i is unique as long as number density and energy density scale differently as a function of the slope. In Fig. <ref> we plot the ratio e/n as a function of  for different energies. For a physically relevant range of slopes there is a one-to-one mapping allowing a unique reconstruction. We tabulate values at the beginning of the simulation to speed up the computation. As the integral for the number density has a simple closed form, we can trivially solve for the amplitude f_i-1/2 analytically and compute it using the new slope _i. We could also use the integral for the energy to find f_i-1/2, however, the integral in general cannot simply be solved for f_i-1/2 analytically.



 §.§ Spectral boundary conditions


In order to conserve CR energy, one would need to apply closed boundaries. In case of losses the energy would accumulate in the lower buffer bin. Adiabatic gains due to strong compression would be stored in the upper buffer bin. Although this seems like a reasonable process for one step it bears difficulties over a longer simulation time. Let us assume a long-term cooling period followed by a strong compression as a simple gedanken experiment. The cooling period will result in an accumulation of significant CR energy in the lowest bin with a large amplitude f_0. Even if the CRs cooled entirely to p=0 (see Fig. <ref>) over that period, they would still be buffered in the lower buffer bin with effectively non-negligible momentum. The strong compression will then push this high-amplitude bin to larger momentum bins resulting in an artificially large population of CRs in bins i that are affected by the compression, cf. equation (<ref>).

We therefore effectively use different boundary conditions for the individual processes. We note that hadronic losses only occur for CR momenta above p_thr, so for spectral configurations with p_min<p_thr the hadronic losses do not interfere with the lower boundary. At the high-p boundary we acknowledge that our numerical spectrum is only a small cut-out of the total CR spectrum, which extends as a power-law tail up to CR momenta of p∼10^20 eV/c. We therefore use inflow boundary conditions, which reflect the continued spectrum towards higher momenta. In practice we simply keep the slope at the highest bin constant during the hadronic cooling step, i.e. it is determined by a combination of the injection spectrum from the previous time step and by the slope from the next lower momentum bin in case of adiabatic compression. We note that energy dependent spatial diffusion would also change the slope inside every bin. However, we currently only account for energy dependent diffusion of each independently, see Section <ref>. For Coulomb losses, we apply outflow boundary conditions at small momenta, resulting in a loss of energy for CRs that cool below p_min, which mimics the
thermalisation process of these CRs. At the high-momentum boundary we follow the same reasoning as for the hadronic losses and compute an inflow of energy based on the continued spectrum outside of our spectral range. For the adiabatic process we assume that the slope in the inflowing buffer bin (lower buffer bin for compression and vice versa) does not change and we allow CRs to enter and leave the spectrum during compression and expansion.



§ ONE-ZONE TESTS













 §.§ Test of the adiabatic process


We would like to stress the importance of an accurate computation of the adiabatic process in hydrodynamical simulations. In compressive turbulent environments the gas cells frequently experience changes of short compression and expansion periods. Even without fully developed turbulence, pressure waves travelling through the simulation domain result in numerous small oscillations.

For this test we only consider adiabatic compression and expansion. Injection of CRs as well as other loss processes are switched off. We compare our new method with a standard finite volume method <cit.> for the piecewise constant representation of the particle distribution function.

As initial condition we use a simple analytic function that mimics a steady state spectrum, i.e. a flat slope at low momenta and a scaling of f∝ p^-4.5 at high momenta,

    f(t=0) = A_0 [(p/1GeV/c)^-a/c + (p/1GeV/c)^-b/c]^-c.

The parameter a is the approximate slope for low momenta (p→0), which we set to a=0. Parameter b is the high-momentum counterpart with a value of -4.5 and c≡2 determines the width of the transition region between the two powerlaw regimes. The overall amplitude of the spectrum, A_0, is set to unity.

Motivated by the strong dynamics in many astrophysical systems, we do not test the methods only for one compression or expansion step but rather for hundreds of periodic oscillations. The divergence of the velocity u is modelled as

    (u)(t) = A_percos(2π t/T).

We set A_per=7.35 and T=1.0. Following equation (<ref>) we find

    p/p = -A_per/3cos(2π t/T)  t,

which yields

    p(t) = p(t_0) exp[-A_perT/6π{sin(2π t/T)-sin(2π t_0/T)}],

where we set t_0=0. After the first quarter of the period, t=0.25, the momentum reduces to 68% of the initial value. After three quarters of the period the maximum momentum is reached with an increase to 148% with respect to the initial value. Over one period the ratio of maximum to minimum momentum p_max/p_min≈2.2, which corresponds to a density compression factor of (p_max/p_min)^3≈10. We note that after a full compression and expansion cycle the spectrum needs to take again the initial shape.

Figure <ref> shows the spectral evolution over 100 periodic cycles for the new method (left-hand panel) and a classical finite-volume method (right-hand panel) using a minmod limiter. The top panel shows the spectrum, the lower one indicates the slope. Colour-coded is the time in units of full cycles. For the classical method we plot every 10^th cycle. In the case of the new method the deviations are so small that we only plot the initial function with a thicker black line and the spectrum after 100 cycles. For the new method we use 12 spectral bins, for the conventional finite volume method we show an example with 500 bins. The flat part of the spectrum is accurately solved with the conventional finite volume method. However, at the transition to the high-momentum powerlaw numerical diffusion causes the spectrum to broaden towards higher momentum with an inaccurate representation of the slope.

In a more quantitative way, we investigate the relative error of the number and energy density as a function of periodic cycles, N_cyc, with respect to the initial value,

    δ_e = |e_CR(N_cyc)-e_CR(0)/e_CR(0)|,

which is plotted in Fig. <ref>. We vary the number of bins from 300 to 1000 for the piecewise constant representation of f. In order to reach a relative error of order unity the conventional method needs more than 1000 bins for a momentum ranging over 5 orders of magnitude. The relative error in energy for the new method is of order 10^-3 for only 12 bins, which demonstrates the superior performance of the new scheme. The error in number density is again several orders of magnitude lower than the energy error and is not shown.

Hadronic and Coulomb losses are best tested in a freely cooling test and a steady state solution. In the case of free cooling we start the simulation with a powerlaw f(p)=10^3(p/(m_pc))^-4.5 and let the spectrum cool including hadronic and Coulomb losses. We compute the spectrum using sub-cycling. Figure <ref> shows the freely cooling spectrum for different times using 12 spectral bins including buffer bins. Also shown are the approximate slopes at large cooling times indicating the accuracy of the method. As expected (Section <ref>) in the case of hadronic cooling the cooled spectrum scales as the original spectrum. The Coulomb losses approach a flat spectrum for f, i.e. a spectral scaling with p^1.9 using the scaling in the plot. There is no analytic solution for the amplitude of the spectrum, so we need to focus on the accuracy of the slopes.

For the steady state spectrum we use again the same initial spectrum as in the case for free cooling, f(p)=10^3(p/(m_pc))^-4.5. We then apply a constant injection j(p)=2.58/Myr(p/(m_pc))^-4.5, i.e. a fraction of 2.58×10^-3 of the energy of the initial spectrum. We use time steps of Δ t =2 Myr and in each time step we apply injection, Coulomb and hadronic losses using 12 bins including buffer bins. The spectrum is shown in the top panel of Fig. <ref> for different times. We overplot the numerical solution at 10 Gyr with the analytical steady state spectrum as in equation (<ref>). The converged spectrum agrees well with the analytic solution. More quantitatively, the bottom panel shows the error of the numerical solution in comparison to the analytic result using the L^1 norm. We adopt a linear measure,

    L^1_lin = N^-1∑_i |f_i,num(p_i)-f_i,ana(p_i)|

as well as its logarithmic counterpart,

    L^1_log = N^-1∑_i |log (f_i,num(p_i)/f_i,ana(p_i))|,

where in both cases we sample the spectrum with 250 data points per momentum decade, giving rise to N=624 sampling points. The error reduces with a scaling of approximately L^1∝ N_bin^-1.3.



§ SPATIAL DIFFUSION



So far we have discussed the spectral evolution without investigating the connection to the spatial time evolution. For advection, this is only a problem of the hydrodynamical code, which is beyond the scope of this paper. Spatial diffusion in contrast reveals an interaction of spatial and spectral changes. As spatial diffusion is energy conserving, the diffusion step itself does not transfer CRs in momentum space, i.e., ∂ f/∂ p=0. However, the amount of diffused number and energy density to neighbouring hydrodynamical cells will depend on the spatial derivatives of f and the possible different diffusion speeds for n and e will result in changes of f_i-1/2 and _i within one bin. This in turn will affect the other processes in the following time step. We have to look at diffusion separately for number density and energy density.



 §.§ Diffusion of number density


The spatial diffusion term of CRs in the momentum range [p_1, p_2] is

    ∂_t n_diff   = ∫_p_1^p_2 4π(D f)p^2  p
       =4π∫_p_1^p_2 p^2 D f  p,

where D is the spatial CR diffusion tensor,

    D   = [ D_11 D_12 D_13; D_21 D_22 D_23; D_31 D_32 D_33 ]

whose components reflect the orientation of the magnetic fields. Following <cit.>, we set the components

    D_ij = D_⊥δ_ij + (D_∥ - D_⊥)b_ib_j

with the normalised magnetic field components b_i=B_i/|B|. The diffusion parameters depend on momentum, so

    D_∥(p)    = D_∥,10p/p_10^α
    
      D_⊥(p)    = D_⊥,10p/p_10^α,

where p_10 = 10 GeV/c and D_∥,10 and D_⊥,10 are diffusion coefficients parallel and perpendicular to the magnetic field line for momentum p_10. We can solve equation (<ref>) directly by replacing the individual components, which results in 

    ∂_t n_diff =
       4π  ∫_p_1^p_2 p^2 [ D_11∂_xf + D_12∂_yf + D_13∂_zf; D_21∂_xf + D_22∂_yf + D_23∂_zf; D_31∂_xf + D_32∂_yf + D_33∂_zf ]  p,

where ∂_k=∂/∂ k is a shorthand notation for the partial derivative with k∈{x, y, z}. In order to solve this equation we need to compute the spatial derivatives of f. Alternatively, we can write the diffusion equation as

    ∂ n_diff/∂ t   =∫_p_1^p_2 4π(D_n f)p^2  p
    
         = D_n n_cr,

such that it formally takes the form of a simple diffusion equation with modified diffusion tensor D_n, see equations (<ref>).
We arrive at that equation in a simple intuitive way formally by multiplying with unity. The individual tensor components take the form

    D_n,ij = ∫_p_1^p_2 p^2 D_ij ∂_j f  p/∫_p_1^p_2 p^2 ∂_j f  p

Using D_ij=D'_ij(p/p_10)^α we get

    D_n,ij = D'_ij/p_10^α ∫_p_1^p_2 p^2+α∂_j f  p/∫_p_1^p_2 p^2 ∂_j f  p




 §.§ Diffusion of energy density


In a similar way the diffusion of CR energy density obeys

    ∂_t e_diff   = ∫_p_1^p_2 4π(D f)p^2T(p)  p
       =4π∫_p_1^p_2 p^2T(p) D f  p,

Using equation (<ref>) we obtain modified diffusion coefficients for the energy

    D_e,ij = ∫_p_1^p_2 p^2T(p) D_ij ∂_j f  p/∫_p_1^p_2 p^2T(p) ∂_j f  p,

which can be connected to the scaling of the diffusion coefficients with the momentum to yield

    D_e,ij = D'_ij/p_10^α ∫_p_1^p_2 p^2+αT(p)∂_j f  p/∫_p_1^p_2 p^2T(p) ∂_j f  p.



 §.§ Simplified bin-centred diffusion






In principle, the correct diffusion can be computed using the previously derived coefficients. However, there are two disadvantages of the computations. One is the computational cost of the spatial derivative because it includes both the spatial derivative of the amplitude f_i-1/2 as well as the slope _i. This expensive derivative has to be done for every cell and every momentum bin separately, which means that no part of the diffusion tensor can be reused for several spectral bins. A second complication is that the full spatial derivative of f has a large dynamic range, in particular near CR sources where the spectrum can differ by orders of magnitude. The strong anisotropy with approximately two orders of magnitude larger parallel than perpendicular diffusion further increases the numerical demands. The derivatives and a stable numerical diffusion are therefore difficult to control in this case. For conserved quantities, slope limiters could be used for the fluxes.

Instead of computing the spatial derivatives of f we can use a simplified approximation and consider the diffusion of number and energy density with the same momentum-dependent diffusion coefficients evaluated with diffusion coefficients for the bin centres of the momenta p_i. This equal treatment of n and e will result in a non-varying slope _i for the diffusion step. This is an approximation that is reasonably stable even in the case of strong shocks and strong CR injection combined with highly anisotropic diffusion.

A few conceptual properties of spatial momentum dependent diffusion can be illustrated in a one-dimensional toy model. We set up an initial Dirac δ distribution in CR number and energy density. We rescale the amplitude in every momentum bin to give a spectral power-law with index -4.5 corresponding to injection of shock-accelerated CRs with a shock of moderate strength. The CR are allowed to diffuse spatially as indicated in the left-hand panel of Fig. <ref> for the spectral bin at p=1 GeV/c, where we applied an arbitrary renormalisation for better readability and use a diffusion coefficient of D_∥,10 = 10^28 cm^2 s^-1 and a scaling of α=0.5, cf equation (<ref>). With this approximation all momentum bins follow the same spatial functional evolution over time, which is a Gaussian function. However, the different diffusion speeds result in faster diffusion for higher CR momenta. The initial power-law spectrum therefore changes its shape over time, which is shown in the right-hand panel of Fig. <ref> for the position at x=150 pc. The spectrum is scaled with p^4.5, i.e. the initial power-law would be a horizontal line. At early times high-energy CRs can reach the position at x=150 pc much faster, which results in larger amplitudes at high p. Over time, the amplitude of the low-energy CRs increases gradually, which reflects their longer diffusion times to the measurement position. At the same time the amplitude of the high-energy CRs starts to rapidly decrease because the high-energy component diffuses out of the domain and the total remaining high-energy component drops. Those two effects result in a temporal evolution of the spectrum that appears in a curved form with a maximum that shifts towards lower momenta. We note that locally in each bin the spectral slope remains the same as in the initial spectrum, =4.5, because the diffusion coefficient scales with the momentum of the bin centre and we use the same diffusion coefficient for the number and the energy density. We would like to emphasise that the changes in the amplitudes of individual bins are solely due to the spatial diffusion. There is no transfer of number and energy density in momentum space, i.e. the spatial integral over each momentum bin is conserved separately. This also means that the diffusion of each momentum bin can be computed independently.

We additionally illustrate this effect in Fig. <ref> for two different times in a three-dimensional representation. The left-hand spectrum is measured after 10 kyr, the right-hand counterpart after 100 kyr. The red lines indicate the spectrum at x=0, the blue lines at 150 pc. After 10 kyr the blue line clearly indicates that the spectrum is dominated by high-energy CRs – we note the large dynamic range in z. At the origin the spectrum still appears to be flat. The decrease of the amplitude at p≳1000 GeV/c is not visible (cf. black line in right-hand panel of Fig <ref>). After 10 kyr the significantly faster diffusion for high p is visible in the spectrum. By that time the diffusion of low-energy CRs results in a shallower distribution along x, which in turn yields a spectrum that develops a stronger curvature at higher momenta.



§ CONCLUSIONS



The spectral shape of CR protons is important to properly account for energy dependent losses, injection, spatial diffusion and CR ionisation and as a result an accurate modelling of the energetic impact of CR protons. Predictive comparisons to observations furthermore requires to follow the spectral shape of a CR population. As CR protons are dynamically relevant, we would like to model them with a full spectrum in every computational cell of a hydrodynamical simulation. This requires an efficient algorithm that adequately represents the CR proton spectrum spanning the relevant energy regime from low-energy CRs that cool via Coulomb losses up to high-energy CR that suffer hadronic losses with a minimum of spectral bins.

Here, we present a new method that solves the Fokker-Planck equation using a piecewise power-law representation of the particle distribution function. The implemented two-moment approach uses the number and energy density to compute the time evolution of the spectrum and does not rely on a continuous particle distribution function, which makes the code more versatile and stable. Because of the low number of spectral bins both the memory requirements and the computational cost are low in comparison to classical methods which require orders of magnitude more bins and adopt piecewise constant values. The method is therefore well suited to be coupled to hydrodynamics and solved together with the gas fluid dynamics in every cell of a three-dimensional hydrodynamical simulation.

For the adiabatic process the new scheme reveals orders of magnitude lower errors for the number and energy density and a very stable and accurate computation of the spectral slope. The combination of injection and cooling including Coulomb and hadronic losses shows very good agreement with theoretical steady state spectra. This method is also able to capture momentum-dependent CR diffusion, which causes a region outside the source to first see the highest-energy CRs before the low-energy CRs catch up.

Besides the dynamics a full spectral representation of CRs allows us to connect the CRs to observables like the chemical changes caused by low-energy CRs or the emission of γ-rays caused by hadronic interactions of high-energy CRs. In our follow-up paper, we explore the coupling of this new method to MHD to study the hydrodynamical impact of evolving the CR spectrum. 



§ ACKNOWLEDGEMENTS

The authors thank Andy Strong and Stefanie Walch for fruitful discussions. We also thank the anonymous referee for very constructive comments that helped to improve the manuscript.
PG and TN acknowledge support from the DFG Priority Program 1573 Physics of the Interstellar Medium.
PG and CP acknowledge funding from the European Research Council under ERC-CoG grant CRAGSMAN-646955.
MH acknowledges support of the (Polish) National Science Centre through the grant No. 2015/19/ST9/02959.
TN acknowledges support from the DFG cluster of excellence ORIGINS.






mnras
 







§ PROOF FOR THE ADIABATIC PROCESS





The adiabatic process is simply an advection in (logarithmic) momentum space plus a change in normalisation. The slope of a power-law distribution, , does not change under the impact of an adiabatic process. For our numerical scheme this means that the ratio e_i/n_i does not change during compression or expansion. We thus take f to be continuous with local slopes ≡_i-1=_i=_i+1 as illustrated in Fig. <ref>. We note that this simple mathematical proof does not work for discontinuous f. This implies that f_i+1/2 = f_i-1/2(p_i+1/2/p_i-1/2)^-. We make a few assumptions without loss of generality:


  * We use the relativistic regime, so T(p)=pc, but the results are similarly true for all other power laws.

  * _i≠4 and ≠3, so that we do not have to treat integrands scaling with 1/p separately.

  * We look at a compression (p_s-1/2 < p_i-1/2) with p_s-1/2 and p_s+1/2 for the left and right edge of bin i and note that the ratios of the shift momenta at both boundaries of the bin i are equal,
  
    p_s-1/2/p_i-1/2=p_s+1/2/p_i+1/2.


  * By construction of a piecewise power-law distribution function and a continuous function we can write

    f_i+1/2 = f_i-1/2(p_i+1/2/p_i-1/2)^-.


We then need to show that

    e_i(t+Δ t)/n_i(t+Δ t) = e_i(t)/n_i(t).

The above simplifications give for the total number density in bin i

    n_i = 4π f_i-1/2p_i-1/2^3/-+3p_i+1/2/p_i-1/2^-+3-1.

The integrated fluxes at p_i-1/2, Φ_n(p_i-1/2), and at p_i+1/2, Φ_n(p_i+1/2) lead to Δ n_i-1/2 and Δ n_i+1/2, which are given by

    Δ n_i-1/2   = 4π f_i-1/2p_i-1/2^3/-+31-p_s-1/2/p_i-1/2^-+3
    Δ n_i+1/2   = 4π f_i+1/2p_i+1/2^3/-+31-p_s+1/2/p_i+1/2^-+3
       = 4π f_i-1/2p_i+1/2/p_i-1/2^- p_i+1/2^3/-+3×
           1-p_s-1/2/p_i-1/2^-+3
       = 4π f_i-1/2p_i-1/2^3/-+3 p_i+1/2/p_i-1/2^-+3×
           1-p_s-1/2/p_i-1/2^-+3

and finally the new number density after compression can be reduced to 

    n_i(t+Δ t)    = n_i + Δ n_i-1/2 - Δ n_i+1/2
    
    
    
    
    
         = n_i p_s-1/2/p_i-1/2^-+3

For the energy density we have the total energy in bin i

    e_i = 4π c f_i-1/2p_i-1/2^4/-+4p_i+1/2/p_i-1/2^-+4-1,

The total changes in the energy can be computed in two steps. The first is to compute the integrals over the distribution function from the shift momenta to the cell boundaries in the same way as for the number density. This is the part of the energy that is transferred between the bins. In a second step we need to consider that the total energy is shifted along the momentum axis and each particle with an energy density e gains energy as

    T   → T p_i-1/2/p_s-1/2,

where we used the assumption that the particle has relativistic energies, i.e. T=pc. This applies to the entire energy bin. We therefore find the energy in bin i as

    e_i(t+Δ t)    = p_i-1/2/p_s-1/2e_i + Δ e_i-1/2 - Δ e_i+1/2

For the first step we compute the integrals as

    Δ e_i-1/2   = ∫_p_s-1/2^p_i-1/2 4π p^3 c f_i-1/2p/p_i-1/2^- p
       = 4π c f_i-1/2p_i-1/2^4/-+4p/p_i-1/2^-+4_p_s-1/2^p_i-1/2
       = 4π c f_i-1/2p_i-1/2^4/-+41-p_s-1/2/p_i-1/2^-+4

and for the integral at the upper bin boundary

    Δ e_i+1/2   = 4π c f_i+1/2p_i+1/2^4/-+41-p_s+1/2/p_i+1/2^-+4
       = 4π c f_i-1/2p_i+1/2/p_i-1/2^-p_i+1/2^4/-+4×
           1-p_s-1/2/p_i-1/2^-+4
       = 4π c f_i-1/2p_i+1/2/p_i-1/2^-+4p_i-1/2^4/-+4×
           1-p_s-1/2/p_i-1/2^-+4
       = 4π c f_i-1/2p_i-1/2^4/-+4p_i+1/2/p_i-1/2^-+4×
           1-p_s-1/2/p_i-1/2^-+4

where in the last equation we have used that f_2 = f_1(p_2/p_1)^- and p_s+1/2/p_i+1/2 = p_s-1/2/p_i-1/2. The difference of those two terms gives

    Δ e_i-1/2 - Δ e_i+1/2   = e_ip_s-1/2/p_i-1/2^-+4-1

Taken together, we find

    e_i + Δ e_i-1/2 - Δ e_i+1/2   =  e_i p_s-1/2/p_i-1/2^-+4.

and finally

    e_i(t+Δ t)    = e_i(t)  p_s-1/2/p_i-1/2^-+3

with the same scaling as the number density compared to initial value e(t). The ratio after an adiabatic compression thus remains the same as well as the slope, ,

    e_i(t+Δ t)/n_i(t+Δ t) = e_i(t)/n_i(t),

which is required for an adiabatic change.





