Describing the cosmic ray propagation model

Criptic allows a very wide range of prescriptions for how cosmic rays propagate through gas; propagation is allowed to depend on arbitrary combinations of position, time, the state of the background gas, the properties of the cosmic ray particles (e.g., the particle energy), and integrals over the cosmic ray populatio (e.g., the cosmic ray pressure gradient) – see the criptic method paper for a full discussion.

In criptic, the way users provide a description for their preferred cosmic ray propagation model is by implementing a class derived from the Propagation class. This class returns a PropagationData struct that contains all the propagation coefficients. This page first describes The PropagationData struct, then covers how to implement a member of The Propagation class, and describes some Pre-defined propagation descriptors. Finally, it discusses Field quantities and their significance.

The PropagationData struct

Criptic describes cosmic ray propagation in terms of a spatial diffusion coefficient parallel to the local magnetic field, a spatial diffusion coefficient perpendicular to the local magnetic field, a momentum diffusion coefficient, a streaming speed, and derivatives of these quantities. The data structure it uses to hold this information is the criptic::propagation::PropagationData struct. The struct includes the following fields:

  • kPar: the spatial diffusion coefficient parallel to the local magnetic field, in astrophysical units (i.e., cm2 s-1 if the code is in GCS mode, m2 s-1 in SI mode)

  • kPerp: the spatial diffusion coefficient perpendicular to the local magnetic field, in astrophysical units

  • kPP: the momentum diffusion coefficient, in internal cosmic ray units (i.e., \((m_p c)^2\) s-1).

  • vStr: the streaming velocity in astrophysical units; the sign convention is that a positive sign corresponds to streaming in the direction parallel to the magnetic field, a negative sign to streaming in the direction anti-parallel

  • dkPP_dp: derivative of kPP with respect to cosmic ray momentum, where both kPP and the cosmic ray momentum are expressed in internal code units (so dKpp_dp has units of \((m_p c)\) s-1)

  • dvStr_dp: derivative of vStr with with respect to momentum, where momentum is in internal code units (has units of cm s-1 / \((m_p c)\) in CGS mode, m s-1 / \((m_p c)\) in SI mode)

  • kParGrad: gradient of kPar with respect to position (units of cm s-1 in CGS mode, m s-1 in SI mode)

  • kPerpGrad: gradient of kPerp with respect to position (units of cm s-1 in CGS mode, m s-1 in SI mode)

  • vStrGrad: gradient of vStr with respect to position (units of s-1)

The Propagation class

The basic object used to describe a model for how cosmic rays propagate is the Propagation class. This is an abstract base class, and user-defined (and pre-defined) models of cosmic ray propagation are implemented as child classes derived from it. The user must instantiate an object of type Propagation and return a pointer to it in the initProp routine called as part of problem setup – see Initializing cosmic ray propagation.

The abstract base class Propagation contains two methods that any derived class must implement. the first is the operator (), which has the call signature:

virtual PropagationData
operator()(const RealVec& x,
           const Real t,
           const gas::GasData& gd,
           const CRPacket& packet,
           const FieldQty& qty,
           const FieldQtyGrad& qtyGrad) const;

Here x is the position, t is the time, gd is the GasData at that position and time, packet is the CRPacket whose propagation coefficient is being computed, and qty and qtyGrad are objects of type FieldQty and FieldQtyGrad, respectively, that represent the Field quantities and their gradients. The implementation provided by a derived class must return a PropagationData struct with every field filled. This struct can depend in an arbitrary way on any of the arguments provided.

The second virtual method that a derived class must implement is:

virtual FieldQtyNeedType fieldQtyNeed() const;

This method returns an enum describing which field quantities the propagation model requires, which can take on the value noFieldQty, needFieldQty, or needFieldQtyGrad. If it returns noFieldQty, then the arguments qty and qtyGrad above will be filled with zeros, and if it returns needFieldQty the argument qtyGrad will be filled with zeros. Only if this method returns needFieldQtyGrad will both qty and qtyGrad be set; see Field quantities for more details.

Pre-defined propagation descriptors

Criptic provides two-predefined propagation models. The first is criptic::propagation::NoPropagation, which simply sets all the entries in the PropagationData struct to zero. This is used mainly for code tests, or for calculations that assume a thick target, so that cosmic ray spatial motion is irrelevant.

The second is criptic::propagation::PropPowerlaw. This class implements a propagation model in which the diffusion coefficients and streaming speed are all assumed to be powerlaw functions of the cosmic ray momentum. Specifically, this method implements a model whereby

\[\begin{split}k_{\parallel} &= k_{\parallel,0} \left(\frac{p}{m_p c}\right)^{k_{\mathrm{idx}, \parallel}} \\ k_{\perp} &= k_{\perp,0} \left(\frac{p}{m_p c}\right)^{k_{\mathrm{idx},\perp}} \\ k_{pp} &= k_{pp,0} \left(\frac{p}{m_p c}\right)^{k_{\mathrm{idx},PP}} \\ \left|v_\mathrm{str}\right| &= v_\mathrm{str,0} \left(\frac{p}{m_p c}\right)^{v_\mathrm{idx,str}}\end{split}\]

The values of the various coefficients on the right hand sides of these expressions are supplied in the parameter file – see Parameter files – via the keywords cr.kPar0, cr.kParIdx, cr.kPerp0, cr.kPerpIdx, cr.kPP0, cr.kPPIdx, cr.vStr0, and cr.vStrIdx. The quantities cr.kPar0 and cr.kPerp0 are in units of cm2 s-1 and the quantity cr.kPP0 is in units of (GeV/\(c\))2. The interpretation of the quantity cr.vStr0 depends on another keyword: cr.vAStream. If this keyword is set to an integer not equal to zero, then cr.vStr0 is dimensionless and gives the streaming velocity as a multiple of the local ion Alfvén speed; if cr.vAStream is zero, then cr.vStr0 is interpreted as an absolute streaming velocity, in units of cm s-1. Finally, the keyword cr.varStreamDir determines the direction of streaming (i.e., the sign of \(v_\mathrm{str}\)). If this is set to zero, streaming is always in the same direction relative to the magnetic field – along the field if cr.vStr0 is positive, anti-parallel to the field if it is negative. If cr.varStreamDir is non-zero, then cosmic rays stream in the direction opposite the cosmic ray pressure gradient. Mathematically, this is implemented as

\[\mathrm{sgn}\left(v_\mathrm{str}\right) = -\mathrm{sgn}\left(\mathbf{B} \cdot \nabla P_\mathrm{cr}\right)\]

Field quantities

As mentioned above, criptic allows cosmic ray propagation to depend on arbitrary combinations of position, time, gas properties, cosmic ray properties, and on three integrals over the cosmic ray field: the cosmic ray number density, the cosmic ray pressure, and the cosmic ray kinetic energy density. We refer to these last three quantities as field quantities, since they involve integrals over the cosmic ray field.

Problems in which propagation depends on field quantities are significantly more computationally expensive than those where it does not, because computing integrals over the cosmic ray field requires performing a kernel density estimation step, as described in the criptic method paper. For this reason, criptic allows implementors of cosmic ray propagation models to declare, via the fieldQtyNeed method (see The Propagation class) of the Propagation class, whether their particular propagation model needs to make use of field quantities. Propagation models can declare that they do not need field quantities at all (by having this method return the enum noFieldQty), that they need access to field quantities but not to their gradients (via returning needFieldQty), or that they need access to both field quantities and their gradients (needFieldQtyGrad). To save unnecessary computation time, criptic will only compute the field quantities if this function returns needFieldQty or needFieldQtyGrad, and will only compute the gradients of field quantities if it returns needFieldQtyGrad.

Limits on propagation parameters

Criptic allows the user to enforce both numerical and physical limits on the propagation parameters. The numerical limits are that the user can specify maximum values of \(k_{\parallel}\), \(k_{\perp}\), and \(k_{pp}\). Any value computed by the propagation model – either pre-supplied or user-specified – that is above the user-specified limits will be flattened to it, and the corresponding gradients and derivatives set to zero. See see Parameter files for details on how to set these limits.

The first physical limits is that Criptic can enforce the Bohm limit, which implies a minimum parallel diffusion coefficient (in CGS units)

\[k_{\parallel,\mathrm{Bohm}} = \frac{r_G v}{3} = \frac{v}{3} \left(\frac{p_\perp}{Z e B}\right),\]

where \(r_G\) is the particle gyro-radius, \(v\) is the particle velocity, \(Z\) is the particle charge in units of the elementary charge \(e\), and \(B\) is the magnetic field strength. The quantity \(p_\perp\) is the perpendicular momentum of the particle, which for simulations that track the pitch angle (see Simulating pitch angle evolution) is set to \(p_\perp = p \sqrt{1-\mu^2}\), where \(p\) is the scalar momentum and \(\mu\) is the pitch angle; for simulations that do not track the CR pitch angle, \(p_\perp = (\pi/4)p\), which is the pitch angle-averaged value.

If the value of \(k_\parallel\) returned by the CR propagation model is smaller than \(k_{\parallel\mathrm{Bohm}}\), the value of \(k_\parallel\) will be increased to \(k_{\parallel\mathrm{Bohm}}\), and the gradient of \(k_\parallel\) will likewise be set to

\[\nabla k_{\parallel,\mathrm{Bohm}} = -\frac{v}{3} \left(\frac{p_\perp}{Z e B^2}\right) \nabla B.\]

Note that the Bohm limit is applied before any numerical limits are enforced, so if \(k_{\parallel,\mathrm{Bohm}} > k_{\parallel,\mathrm{max}}\), then the parallel diffusion coefficient will be set to \(k_{\parallel,\mathrm{max}}\), not \(k_{\parallel,\mathrm{Bohm}}\).

The second physical limit is that Criptic can limit the drift speed to be \(<c\) – or, if the user prefers, some multiple of \(c\), usually chosen to be less than unity, in order to avoid the time step becoming too small. The drift speeds in the frame co-moving with the background gas are

\[\mathbf{v}_\mathrm{drift} = \left[\nabla k_\parallel \cdot \hat{\mathbf{t}} + v_\mathrm{str} + \frac{p}{3}\frac{\partial v_\mathrm{str}}{\partial p}\right] \hat{\mathbf{t}} + \left(\nabla k_\perp \cdot\hat{\mathbf{n}}\right) \hat{\mathbf{n}} + \left(\nabla k_\perp \cdot\hat{\mathbf{b}}\right) \hat{\mathbf{b}},\]

where \(\hat{\mathbf{t}}\), \(\hat{\mathbf{n}}\), and \(\hat{\mathbf{b}}\) are the TNB basis vectors for the local magnetic field, \(v_\mathrm{str}\) is the local streaming speed, and \(p\) is the CR momentum. The user can supply a numerical value \(beta_\mathrm{max}\) that limits the drift speed to \(\beta_\mathrm{max} c\). This limit is applied independently to each of the terms that contribute to the drift velocity and is applied so as to keep the ratio of parallel to perpendicular drift constant. Specifically, the limits applied are, in order:

  1. If \(|\nabla k_\parallel| > \beta_\mathrm{max} c\), multiply both \(\nabla k_\parallel\) and \(\nabla k_\perp\) by a factor \(\beta_\mathrm{max} c/ |\nabla k_\parallel|\).

  2. If \(|\nabla k_\perp| > \beta_\mathrm{max} c\), multiply both \(\nabla k_\parallel\) and \(\nabla k_\perp\) by a factor \(\beta_\mathrm{max} c/|\nabla k_\perp|\).

  3. If \(|v_\mathrm{str} + (p/3)(\partial v_\mathrm{str}/\partial p)| > \beta_\mathrm{max} c\), multiply \(v_\mathrm{str}\) and \(\partial v_\mathrm{str}/\partial p\) by a factor \(\beta_\mathrm{max} c/|v_\mathrm{str} + (p/3)(\partial v_\mathrm{str}/\partial p)|\)

By default \(\beta_\mathrm{max}\) is set to unity, but this can be overridden by the user – see Parameter files.