.. highlight:: rest

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 <https://ui.adsabs.harvard.edu/abs/2022arXiv220713838K/abstract>`_ 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 <doxyhtml/classcriptic_1_1propagation_1_1_propagation.html>`_ class. This class returns a `PropagationData <doxyhtml/structcriptic_1_1propagation_1_1_propagation_data.html>`_ struct that contains all the propagation coefficients. This page first describes :ref:`ssec-propdata`, then covers how to implement a member of :ref:`ssec-propagation`, and describes some :ref:`ssec-prop-predef`. Finally, it discusses :ref:`ssec-field-quantities` and their significance.

.. _ssec-propdata:

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 <doxyhtml/structcriptic_1_1propagation_1_1_propagation_data.html>`_ struct. The struct includes the following fields:

* ``kPar``: the spatial diffusion coefficient parallel to the local magnetic field, in astrophysical units (i.e., cm\ :sup:`2` s\ :sup:`-1` if the code is in GCS mode, m\ :sup:`2` s\ :sup:`-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., :math:`(m_p c)^2` s\ :sup:`-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 :math:`(m_p c)` s\ :sup:`-1`)
* ``dvStr_dp``: derivative of ``vStr`` with with respect to momentum, where momentum is in internal code units (has units of cm s\ :sup:`-1` / :math:`(m_p c)` in CGS mode, m s\ :sup:`-1` / :math:`(m_p c)` in SI mode)
* ``kParGrad``: gradient of ``kPar`` with respect to position (units of cm s\ :sup:`-1` in CGS mode, m s\ :sup:`-1` in SI mode)
* ``kPerpGrad``: gradient of ``kPerp`` with respect to position (units of cm s\ :sup:`-1` in CGS mode, m s\ :sup:`-1` in SI mode)
* ``vStrGrad``: gradient of ``vStr`` with respect to position (units of s\ :sup:`-1`)

.. _ssec-propagation:

The ``Propagation`` class
-------------------------

The basic object used to describe a model for how cosmic rays propagate is the `Propagation <doxyhtml/classcriptic_1_1propagation_1_1_propagation.html>`_ 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 <doxyhtml/classcriptic_1_1propagation_1_1_propagation.html>`_ and return a pointer to it in the ``initProp`` routine called as part of problem setup -- see :ref:`ssec-init-propagation`.

The abstract base class `Propagation <doxyhtml/classcriptic_1_1propagation_1_1_propagation.html>`_ 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 <doxyhtml/classcriptic_1_1gas_1_1_gas_data.html>`_ at that position and time, ``packet`` is the `CRPacket <Doxygen/html/classcriptic_1_1_c_r_packet.html>`_ whose propagation coefficient is being computed, and ``qty`` and ``qtyGrad`` are objects of type `FieldQty <doxyhtml/classcriptic_1_1_field_qty.html>`_ and `FieldQtyGrad <doxyhtml/classcriptic_1_1_field_qty_grad.html>`_, respectively, that represent the :ref:`ssec-field-quantities` and their gradients. The implementation provided by a derived class must return a `PropagationData <doxyhtml/structcriptic_1_1propagation_1_1_propagation_data.html>`_ 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 :ref:`ssec-field-quantities` for more details.

.. _ssec-prop-predef:

Pre-defined propagation descriptors
-----------------------------------

Criptic provides two-predefined propagation models. The first is `criptic::propagation::NoPropagation <doxyhtml/classcriptic_1_1propagation_1_1_no_propagation.html>`_, which simply sets all the entries in the `PropagationData <doxyhtml/structcriptic_1_1propagation_1_1_propagation_data.html>`_ 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 <doxyhtml/classcriptic_1_1propagation_1_1_prop_powerlaw.html>`_. 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

.. math::

   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}}

The values of the various coefficients on the right hand sides of these expressions
are supplied in the parameter file -- see :doc:`parameters` -- 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 cm\ :sup:`2` s\ :sup:`-1` and the quantity ``cr.kPP0``
is in units of (GeV/:math:`c`)\ :sup:`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\ :sup:`-1`. Finally, the keyword ``cr.varStreamDir`` determines the
direction of streaming (i.e., the sign of :math:`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

.. math::

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

.. _ssec-field-quantities:

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 <https://ui.adsabs.harvard.edu/abs/2022arXiv220713838K/abstract>`_. For this reason, criptic allows implementors of cosmic ray propagation models to declare, via the ``fieldQtyNeed`` method (see :ref:`ssec-propagation`) of the `Propagation <doxyhtml/classcriptic_1_1propagation_1_1_propagation.html>`_ 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``.

.. _ssec-limits:

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 :math:`k_{\parallel}`, :math:`k_{\perp}`, and :math:`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 :doc:`parameters` 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)

.. math::

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

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

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

.. math::

   \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 :math:`k_{\parallel,\mathrm{Bohm}} > k_{\parallel,\mathrm{max}}`, then the
parallel diffusion coefficient will be set to :math:`k_{\parallel,\mathrm{max}}`,
not :math:`k_{\parallel,\mathrm{Bohm}}`.

The second physical limit is that Criptic can limit the drift speed to be
:math:`<c` -- or, if the user prefers, some multiple of :math:`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

.. math::

   \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 :math:`\hat{\mathbf{t}}`, :math:`\hat{\mathbf{n}}`, and 
:math:`\hat{\mathbf{b}}` are the TNB basis vectors for the local
magnetic field, :math:`v_\mathrm{str}` is the local streaming speed,
and :math:`p` is the CR momentum. The user can supply a numerical value
:math:`beta_\mathrm{max}` that limits the drift speed to
:math:`\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 :math:`|\nabla k_\parallel| > \beta_\mathrm{max} c`, multiply both
   :math:`\nabla k_\parallel` and :math:`\nabla k_\perp` by a factor
   :math:`\beta_\mathrm{max} c/ |\nabla k_\parallel|`.
2. If :math:`|\nabla k_\perp| > \beta_\mathrm{max} c`, multiply both
   :math:`\nabla k_\parallel` and :math:`\nabla k_\perp` by a factor
   :math:`\beta_\mathrm{max} c/|\nabla k_\perp|`.
3. If :math:`|v_\mathrm{str} + (p/3)(\partial v_\mathrm{str}/\partial p)| > \beta_\mathrm{max} c`,
   multiply :math:`v_\mathrm{str}` and :math:`\partial v_\mathrm{str}/\partial p`
   by a factor :math:`\beta_\mathrm{max} c/|v_\mathrm{str} + (p/3)(\partial v_\mathrm{str}/\partial p)|`

By default :math:`\beta_\mathrm{max}` is set to unity, but this can be overridden
by the user -- see :doc:`parameters`.