Advection fields for the Monte Carlo Propagators

[ehlertdo]

Hey Background I am looking at scenarios where the cosmic rays propagate close to the source in environments where advective transport dominates over diffusion for $E\lesssim1~\mathrm{EeV}$. As far as I can tell, the advection fields currently implemented in CRPropa are only meant to be used with the DiffusionSDE solver. However, it appears to me that it's not possible to include interactions (photopion, photodisintegration, etc.) for that solver. Please correct me if I'm wrong on either of these two points. My current approach I have modified the PropagationBP code to accept and apply advection fields in a very approximate manner (simply move the particle by $t\mathrm{step} \cdot v\mathrm{adv}$ each step). See PropagationBP.txt ( I can also create a pull request if necessary). It's clear that this solution is not exact as the final propagation distance of particles (column 'D' in the output files) can become shorter than the system size. I assume that this occurs since the extra movement due to the advection is in addition to the regular step size. This also leads me to believe that the interactions are likely underestimated since the propagation distance assumed by CRPropa is smaller than the true path length. In addition, simply adding the normal/diffusive motion and advection does not respect $v_\mathrm{CR}\leq c$ by default. My question / request Is there a way to include the advection into either the Boris-Push or Cash-Carp solvers directly instead of adding it afterward as a separate step? (Alternatively my problem can also be solved if there is a way to include interactions with the DiffusionSDE module).

Best, Domenik

[ehlertdo]

Short update: I have now included the advection in a way that is more realistic by taking it into account in terms of the "effective" propagation direction of the cosmic ray instead of simply adding the advection on top of the usual propagation as an additional step. See also the attached slide for more explanation. If anyone has some comments on this that would be appreciated.

PropagationBP::Y PropagationBP::dY(Vector3d pos, Vector3d dir, double step,
            double z, double q, double m) const {
        // velocity of advection field [m/s]
        Vector3d vWind = getAdvFieldAtPosition(pos);
        // add advection vector [(vx vy vz) / c] to propagation direction
        Vector3d dir_tot = dir + vWind.getUnitVector() * (vWind.getR() / c_light);
        // renormalise to unit vector
        dir_tot = dir_tot.getUnitVector();

        // BORIS-PUSH ALGORITHM ========================
        // half leap frog step in the position
        pos += dir_tot * step / 2.;

        // get B field at particle position
        Vector3d B = getFieldAtPosition(pos, z);

        // Boris help vectors
        Vector3d t = B * q / 2 / m * step / c_light;
        Vector3d s = t * 2 / (1 + t.dot(t));
        Vector3d v_help;

        // Boris push
        v_help = dir + dir.cross(t);
        dir = dir + v_help.cross(s);

        // include advection for the second half step
        dir_tot = dir + vWind.getUnitVector() * (vWind.getR() / c_light);
        dir_tot = dir_tot.getUnitVector();

        // the other half leap frog step in the position
        pos += dir_tot * step / 2.;

        return Y(pos, dir);
    }```