# Extending the Cooperative Dual-Task Space
in Conformal Geometric Algebra

Tobias Löw12 and Sylvain Calinon12  1Idiap Research Institute, Martigny,
Switzerland  2EPFL, Lausanne, Switzerland

###### Abstract

In this work, we are presenting an extension of the cooperative dual-task
space (CDTS) in conformal geometric algebra. The CDTS was first defined using
dual quaternion algebra and is a well established framework for the simplified
definition of tasks using two manipulators. By integrating conformal geometric
algebra, we aim to further enhance the geometric expressiveness and thus
simplify the modeling of various tasks. We show this formulation by first
presenting the CDTS and then its extension that is based around a cooperative
pointpair. This extension keeps all the benefits of the original formulation
that is based on dual quaternions, but adds more tools for geometric modeling
of the dual-arm tasks. We also present how this CGA-CDTS can be seamlessly
integrated with an optimal control framework in geometric algebra that was
derived in previous work. In the experiments, we demonstrate how to model
different objectives and constraints using the CGA-CDTS. Using a setup of two
Franka Emika robots we then show the effectiveness of our approach using model
predictive control in real world experiments.

###### Index Terms:

Geometric Algebra, Dual-Arm Manipulation, Optimal Control

##  I INTRODUCTION

This work was supported by the State Secretariat for Education, Research and
Innovation in Switzerland for participation in the European Commission’s
Horizon Europe Program through the INTELLIMAN project (https://intelliman-
project.eu/, HORIZON-CL4-Digital-Emerging Grant 101070136) and the SESTOSENSO
project (http://sestosenso.eu/, HORIZON-CL4-Digital-Emerging Grant 101070310).

With the increasing desire to deploy robots in human environemnts, the need
for robots to have human-like manipulation capabilities arises. One inherent
ability that humans have is to manipulate objects using both their hands and
arms, which is needed for example when objects that are either too large or
too heavy need to be manipulated. In order to match these capabilities and to
be able to mimic them, robotic systems also need to be able to cooperatively
control two arms in order to perform tasks in human environments.

Apart from dual-arm systems being more human-like in terms of form factor,
they also have some technical advantages. One can have the stiffness and
strength of parallel manipulators combined with the flexibility and dexterity
of serial manipulators [1]. Furthermore, since they increase the redundancy in
the task-space due to their high number of degrees of freedom, they are better
suited for intricate tasks that require a high manipulability such as screw
assembly [2] and dishwashing [3]. Other applications of bimanual systems are
manipulating articulated objects [4] or cables [5]. More advantages and
examples are listed in [6].

Since many problems in robotics boil down to optimization problems that can be
solved efficiently with various state-of-the-art solvers, it is of great
interest to facilitate the modeling and increase the expressiveness of the
formulations. Choosing the correct representation can make a huge difference
in terms of how much prior knowledge we can embed into the formulation of
those optimization problems. These are in robotics often very geometric, hence
it is very beneficial to choose representations that intuitively allow to
incoporate the geometry of the problem. For the case of dual-arm manipulation,
the cooperative dual-task space (CDTS) [7] was proposed. This approach uses
dual quaternion algebra (DQA), which not only unifies the treatment of
position and orientation, it also allows the representation of various
geometric primitives [8]. These primitives can then be used to simplify the
modeling of the tasks.

Figure 1: Cooperative dual-task space using conformal geometric algebra. The
figure shows the two manipulators as well as their individual, relative and
absolute motors. Additionally, it shows the cooperative pointpair.

Dual quaternions have a strong connection to geometric algebra, especially the
variants known as projective (PGA) and conformal (CGA) geometric algebra [9],
since dual quaternions are isomorphically embedded in their sub-algebras [10].
The geometric algebras, however, are richer algebras that offer more geometric
primitives and, more importantly, they offer the geometric construction of
primitives based on operations such as intersections [10]. This leads to new
possibilities when formulating objectives and constraints based on the
geometric primitives in the CGA-CDTS compared to the DQ-CDTS, which we will
show in the experiments.

In this article, we formulate the CDTS in conformal geometric algebra and show
how this formulation naturally extends the DQ-CDTS. The resulting CGA-CDTS
retains the same properties of a compact representation of two arm system as
the DQ-CDTS, while adding a useful geometric primitive, the cooperative
pointpair, that represents both end-effector positions simultaneously.
Furthermore, we demonstrate how the CGA-CDTS can be used in the optimal
control formulation with geometric primitives for manipulators that we
presented in [11]. Hence, this article aims to explain the basic mathematical
formulations of various optimization problems using the CGA-CDTS.

##  II BACKGROUND

###  II-A Cooperative Dual-Task Space

We will introduce the cooperative dual-task space here conceptually, and not
mathematically, since its original definition uses dual quaternion algebra. In
contrast to that, we will be defining and extending it using conformal
geometric algebra. So, for the sake of conciseness, we are keeping this
introduction on a high level.

The cooperative dual-task space was proposed as a compact and singularity-free
representation of a two-arm system [7]. It is defined using two poses that
depend on the end-effector poses of the two manipulators, one is the relative
and the other one is the absolute pose. All poses are represented using unit
dual quaternions, which have advantages compared to other representations,
such as coupled position and orientation, singularity-free representation and
efficient computation [12].

Based on the compact representation of the CDTS, various control strategies
have been proposed. In [13], a coupled task-space admittance controller was
presented, that allowed for a geometrically consistent stiffness term. A
reactive control strategy was developed in [14] that leveraged geometric
primitves for task relaxations and priorities. Task priorities for control in
the CDTS were also proposed in [15]. In order to exploit human demonstration
that allow the teaching of cooperative motions, motion primitives for bimanual
systems were presented based on the CDTS [16].

Note that the results of the research that is based on the DQ-CDTS can also be
used with the CGA-CDTS, albeit with mathematical changes due to the different
algebra. Furthermore, the mentioned advantages of DQA also apply to CGA, since
dual quaternions and the corresponding subalgebra in CGA, i.e. the motors, are
isomorphic [17].

###  II-B Geometric Algebra

Geometric algebra is a single algebra for geometric reasoning, alleviating the
need of utilizing multiple algebras to express geometric relations. In this
article, we are using the variant known as conformal geometric algebra (CGA).
We use the following notation: $x$ to denote scalars, $\bm{x}$ for vectors,
$\bm{X}$ for matrices, $X$ for multivectors and $\bm{\mathcal{X}}$ for
matrices of multivectors.

A general element for computation in geometric algebra is called a
multivector. There are three main products that can be used with multivectors:
the geometric product $XY$, the inner product $X\cdot Y$ and the outer product
$X\wedge Y$. The trivial vector case shows that the geometric product combines
the inner $\cdot$ and the outer $\wedge$ product

$\bm{a}\bm{b}=\bm{a}\cdot\bm{b}+\bm{a}\wedge\bm{b}.$ (1)

The outer product is a spanning operation that allows the creation of
geometric primitives from points $P_{i}$. For example, two points $P_{1}\wedge
P_{2}$ yield a point pair, three points $P_{1}\wedge P_{2}\wedge P_{3}$ a
circle and four points $P_{1}\wedge P_{2}\wedge P_{3}\wedge P_{4}$ a sphere.
There are more primitives, that we will introduce when we use them in the
experiments.

Rigid body transformations are represented by motors $M$. They can be applied
to multivectors by a sandwiching operation, similar to how quaternions rotate
vectors

$Y=MX\widetilde{M},$ (2)

where $\widetilde{M}$ stands for the reverse of a motor. Motors are
exponential mappings of so-called bivector (i.e. the subspaces spanned by the
outer product of two vectors), the inverse operation is the logarithmic map

$M=\exp(B)\hskip 8.53581pt\iff\hskip 8.53581ptB=\log(M).$ (3)

The motor $M(\bm{q})$ corresponding to the forward kinematics of a kinematic
chain can be computed as the product of the individual joint motors

$M(\bm{q})=\prod_{i=1}^{N}M_{i}(q_{i}).$ (4)

The analytic Jacobian can then be found as the derivative of the forward
kinematics motor defined in Equation (4), i.e.

$\bm{\mathcal{J}}^{A}(\bm{q})=\frac{\partial
M(\bm{q})}{\partial\bm{q}}=\begin{bmatrix}\frac{\partial M(\bm{q})}{\partial
q_{1}}\ldots\frac{\partial M(\bm{q})}{\partial q_{N}}\end{bmatrix}.$ (5)

At this point, we only introduced the most important concepts for
understanding the proposed method. For a complete introduction we refer
interested readers to [18] and [19].

##  III METHOD

The CDTS in conformal geometric algebra is an extension of the CDTS in dual
quaternion algebra. We first present the basic reformulation of the CDTS in
CGA and then its extension by using the additional geometric primitives.
Lastly, we present how the CGA-CDTS can be used within an optimal control
framework using geometric algebra that we previously proposed for manipulation
tasks.

###  III-A Conformal Geometric Algebra Cooperative Dual-Task Space

The geometric algebra equivalent of the relative and absolute dual quaternions
of the DQ-CDTS are defined using motors. Given the joint configurations of the
two manipulators, $\bm{q}_{1}$ and $\bm{q}_{2}$, respectively, we can easily
find their end-effector motors $M_{1}(\bm{q}_{1})$ and $M_{2}(\bm{q}_{2})$
using the forward kinematics in CGA. From this it is straightforward to
formulate the relative motor as

$M_{r}(\bm{q}_{1},\bm{q}_{2})=\widetilde{M}_{2}(\bm{q}_{2})M_{1}(\bm{q}_{1}),$
(6)

while its Jacobian, i.e. the relative analytic Jacobian, can be found as

$\bm{\mathcal{J}}^{A}_{r}(\bm{q}_{1},\bm{q}_{2})=\begin{bmatrix}\widetilde{M}_{2}(\bm{q}_{2})\bm{\mathcal{J}}^{A}_{1}(\bm{q}_{1})&\bm{\mathcal{\widetilde{J}}}^{A}_{2}(\bm{q}_{2})M_{1}(\bm{q}_{1})\end{bmatrix}.$
(7)

Similarly, the absolute motor can be found as

$\begin{split}M_{a}(\bm{q}_{1},\bm{q}_{2})&=M_{2}(\bm{q}_{2})M_{r/2}(\bm{q}_{1},\bm{q}_{2})\\\
&=M_{2}(\bm{q}_{2})\exp\left(\frac{1}{2}\log\Big{(}M_{r}(\bm{q}_{1},\bm{q}_{2})\Big{)}\right),\end{split}$
(8)

with its corresponding absolute analytic Jacobian

$\begin{split}\bm{\mathcal{J}}^{A}_{a}(\bm{q}_{1},\bm{q}_{2})=&M_{2}(\bm{q}_{2})\bm{\mathcal{J}}^{A}_{M_{r/2}}(\bm{q}_{1},\bm{q}_{2})\\\
&+\begin{bmatrix}\bm{0}&\bm{\mathcal{J}}^{A}_{2}(\bm{q}_{2})M_{r/2}(\bm{q}_{1},\bm{q}_{2})\end{bmatrix},\end{split}$
(9)

where $\bm{\mathcal{J}}^{A}_{M_{r/2}}(\bm{q}_{1},\bm{q}_{2})$ is the analytic
Jacobian of the motor $M_{r/2}(\bm{q}_{1},\bm{q}_{2})$. It can be found as

$\bm{\mathcal{J}}^{A}_{M_{r/2}}(\bm{q}_{1},\bm{q}_{2})=\bm{J}_{\mathbb{B}\to\mathcal{M}}(B_{r/2})\bm{J}_{\mathcal{M}\to\mathbb{B}}(M_{r})\bm{J}^{A}_{r}(\bm{q}_{1},\bm{q}_{2}).$
(10)

Here, $B_{r/2}$ is the logarithm of the motor $M_{r/2}$. The matrices
$\bm{J}_{\mathbb{B}\to\mathcal{M}}$ and $\bm{J}_{\mathcal{M}\to\mathbb{B}}$
are the Jacobians of the exponential and logarithmic mapping respectively. We
already showed the derivation of the Jacobian of the logarithmic mapping in
the appendix of [11]. The Jacobian of the exponential mapping can be found in
Appendix A.

Both the relative and the absolute analytic Jacobians are $1\times 2N$
multivector matrices that contain motors as their elements. Hence, when
expanding it to normal matrix algebra they become $8\times 2N$ matrices.

Since motors in CGA can be used to transform any geometric primitive that is
part of the algebra in a uniform way, it is easy to find cooperative geometric
primitives. Their definition can be trivially found using Equation (2), where
$M$ is either the relative $M_{r}(\bm{q}_{1},\bm{q}_{2})$ or absolute
$M_{a}(\bm{q}_{1},\bm{q}_{2})$ motor and $X$ can be any geometric primitive.
The corresponding Jacobians are then found using the respective Jacobians
$\bm{\mathcal{J}}^{A}_{r}(\bm{q}_{1},\bm{q}_{2})$ and
$\bm{\mathcal{J}}^{A}_{a}(\bm{q}_{1},\bm{q}_{2})$.

###  III-B Cooperative Pointpair

In extension to the CDTS that was defined using dual quaternion algebra, the
CDTS presented here using CGA also allows a geometric primitive that
corresponds to both end-effector positions simultaneously. This cooperative
pointpair is defined as the outer product of the two end-effector points, i.e.

$P_{cdts}=M_{1}(\bm{q}_{1})\bm{e}_{0}\widetilde{M}_{1}(\bm{q}_{1})\wedge
M_{2}(\bm{q}_{2})\bm{e}_{0}\widetilde{M}_{2}(\bm{q}_{2}).$ (11)

The Jacobian of the cooperative pointpair can be found as

$\bm{\mathcal{J}}_{P_{cdts}}=\begin{bmatrix}\bm{\mathcal{J}}_{P_{cdts},1}&\bm{\mathcal{J}}_{P_{cdts},2}\par\end{bmatrix},$
(12)

where

$\begin{split}\bm{\mathcal{J}}_{P_{cdts},1}=&\bm{\mathcal{J}}^{A}_{1}(\bm{q}_{1})\bm{e}_{0}\widetilde{M}_{1}(\bm{q}_{1})\wedge
M_{2}(\bm{q}_{2})\bm{e}_{0}\widetilde{M}_{2}(\bm{q}_{2})\\\
&+M_{1}(\bm{q}_{1})\bm{e}_{0}\bm{\mathcal{\widetilde{J}}}^{A}_{1}(\bm{q}_{1})\wedge
M_{2}(\bm{q}_{2})\bm{e}_{0}\widetilde{M}_{2}(\bm{q}_{2}),\end{split}$ (13)

and

$\begin{split}\bm{\mathcal{J}}_{P_{cdts},2}=&M_{1}(\bm{q}_{1})\bm{e}_{0}\widetilde{M}_{1}(\bm{q}_{1})\wedge\bm{\mathcal{J}}^{A}_{2}(\bm{q}_{2})\bm{e}_{0}\widetilde{M}_{2}(\bm{q}_{2})\\\
&+M_{1}(\bm{q}_{1})\bm{e}_{0}\widetilde{M}_{1}(\bm{q}_{1})\wedge
M_{2}(\bm{q}_{2})\bm{e}_{0}\bm{\mathcal{\widetilde{J}}}^{A}_{2}(\bm{q}_{2}).\end{split}$
(14)

Note that the cooperative pointpair is a direct representation of both points
and is not the same as stacking the two points. Therefore the Jacobian matrix
is also different, which will lead to different solutions of optimization
problems. An example of this is shown in Figure 2, where two Franka Emika
robots are tasked to reach a plane, once individually (i.e. by stacking their
end-effector points) and once cooperatively (i.e. by using the cooperative
pointpair that is presented here). It can be seen that the corresponding
solution configurations are not the same, which shows that the cooperative
pointpair representation lets the two robots influence each other.

Figure 2: Dual-Arm manipulator reaching a plane. The target plane is shown in
red. The white Franka Emika robots show the initial configurations, the green
ones are the result of individually reaching the plane, and the blue ones
cooperatively.

Evidently, this cooperative pointpair Jacobian has a singularity in the case
when both end-effector positions are equal. This will cause the outer product,
by definition, to be zero. In practice, this can be avoided by posing a
constraint on the distance between end-effectors, which usually is required
anyways.

###  III-C Using the CGA-CDTS for Optimal Control

Integrating the CGA-CDTS into the optimal control framework that we presented
in [11] can be achieved by replacing the single end-effector motor with the
relative and absolute motors, respectively. The objective of reaching a target
motor is then formulated as

$E_{M}(\bm{q}_{1},\bm{q}_{2})=\log\left(\widetilde{M}_{target}M(\bm{q}_{1},\bm{q}_{2})\right),$
(15)

where $M(\bm{q}_{1},\bm{q}_{2})$ can be either $M_{r}(\bm{q}_{1},\bm{q}_{2})$
or $M_{a}(\bm{q}_{1},\bm{q}_{2})$.

Alternatively, we can define a residual multivector for reaching a geometric
primitive $X_{d}$ as

$E_{X_{d}}(\bm{q}_{1},\bm{q}_{2})=X_{d}\wedge
M(\bm{q}_{1},\bm{q}_{2})X\widetilde{M}(\bm{q}_{1},\bm{q}_{2}),$ (16)

again using either the relative or absolute motor. Deriving the respective
Jacobians is straightforward using the definitions of the relative and
absolute analytic Jacobians in Equations (7) and (9). With this, these
residual multivectors of the CGA-CDTS can then directly be used to define
objectives or constraints for optimal control problems, which would mean for
example that a relative or absolute geometric primitive should be reached.

The cooperative pointpair can also be used in order to define tasks as
optimization problems. The first way to do so is using the outer product in a
similar way to the above relative and absolute residual multivectors, i.e.

$E_{cdts}(\bm{q}_{1},\bm{q}_{2})=P_{target}\wedge
P_{cdts}(\bm{q}_{1},\bm{q}_{2}),$ (17)

with the Jacobian

$\bm{\mathcal{J}}_{E_{cdts}}(\bm{q}_{1},\bm{q}_{2})=P_{target}\wedge\bm{\mathcal{J}}_{P_{cdts}}(\bm{q}_{1},\bm{q}_{2}),$
(18)

This objective means that the two manipulators should cooperatively reach a
single point. The implications and results of this objective are further
detailed in the experiment section.

Another common use-case are containment relationships for the cooperative
pointpair with respect to other geometric primitives. These can then be used
to define tasks where the dual-arm system should cooperatively reach a target.
The mathematical formulation is using a product called the commutator product
$\times$, i.e.

$\begin{split}E_{cdts}(\bm{q}_{1},\bm{q}_{2})&=X_{d}\times
P_{cdts}(\bm{q}_{1},\bm{q}_{2})\\\
&=\frac{1}{2}(X_{d}P_{cdts}(\bm{q}_{1},\bm{q}_{2})-P_{cdts}(\bm{q}_{1},\bm{q}_{2})X_{d}).\end{split}$
(19)

The Jacobian of this containment relationship can be found as

$\bm{\mathcal{J}}_{E_{cdts}}(\bm{q}_{1},\bm{q}_{2})=X_{d}\times\bm{\mathcal{J}}_{E_{cdts}}(\bm{q}_{1},\bm{q}_{2}).$
(20)

The interpretation of the residual multivector
$E_{cdts}(\bm{q}_{1},\bm{q}_{2})$ is that it should be reduced to zero if the
cooperative pointpair is contained within the desired geometric primitive
$X_{d}$, e.g. if both end-effector points lie on a circle.

The distance between the two end-effector can be constrained using the inner
product between the two end-effector points, i.e.

$E_{d}(\bm{q}_{1},\bm{q}_{2})=-2M_{1}(\bm{q}_{1})\bm{e}_{0}\widetilde{M}_{1}(\bm{q}_{1})\cdot
M_{2}(\bm{q}_{2})\bm{e}_{0}\widetilde{M}_{2}(\bm{q}_{2})-d^{2},$ (21)

where $d$ is the desired distance.

##  IV EXPERIMENTS

In this section, we are presenting various cooperative tasks that are defined
in the CGA-CDTS. For each task we are providing the mathematical definition of
the optimization problem. We are then solving those optimization problems
using standard solvers such as Gauss-Newton. For simplicity, the problems in
simulation are formulated essentially as inverse kinematics problems, the same
objectives and constraints can, however, be used in optimal control problems
that are for example then used for model predictive control. We are
demonstrating this in the real world experiments, where the problems are then
solved by using a variant of the iterative linear quadratic regulator (iLQR)
[20]. Both the simulation and the real-world experiments use the same setup of
two table-top mounted Franka Emika robots. Additional material for the
experiments as well as the videos of the real world experiments can be found
on the accompanying website 111https://geometric-algebra.tobiloew.ch/cdts/.
All geometric algebra computations are done using our open-source library
gafro222https://github.com/idiap/gafro.

###  IV-A Cooperatively Reaching a Point

As mentioned before, the cooperative pointpair models both end-effector
positions simultaneously. Hence, when formulating a simple optimization
problem such as reaching a single point, the system automatically chooses
which manipulator should perform the task. Since the information of both
manipulators is encoded in one geometric primitive, it alleviates the need to
construct the problem using conditional formulations. The problem can be
expressed compactly using the outer product, i.e.

$\begin{split}\bm{q}^{*}=\min_{\bm{q}}&\Big{\|}P_{target}\wedge
P_{cdts}(\bm{q}_{1},\bm{q}_{2})\Big{\|}_{2}^{2}.\end{split}$ (22)

An example of this task is shown in Figure 3. Notice how in both cases the
manipulator that is closer to the point performs the reaching task while the
other remains in its initial configuration.

(a) Left arm reaching for the point.

(b) Right arm reaching for the point.

Figure 3: Two Franka Emika robots reaching for a single point in the
cooperative dual-task space. The initial configurations are shown in white and
the final ones in gray. The target point is shown in red.

We also show this experiment in using the real-world setup and the
accompanying video can be found on our website.

###  IV-B Cooperatively Reaching a Circle

One of the geometric primitives that is available in CGA but not DQA is a
circle. Hence, we can use the CGA-CDTS to define a task where a dual arm
system should cooperatively reach a circle. An example application of reaching
a circle would be holding a filled bucket with two manipulators.
Mathematically, a circle is obtained by the outer product of three points.

The problem of cooperatively reaching a circle is formulated as a constrained
optimization problem using the cooperative pointpair

$\begin{split}\bm{q}^{*}=\min_{\bm{q}}&\Big{\|}\log\left(\widetilde{M}_{r}(\bm{q}_{0})M_{r}(\bm{q})\right)\Big{\|}_{2}^{2}\\\
&\text{s.t. }C\times P_{cdts}(\bm{q}_{1},\bm{q}_{2})=0.\end{split}$ (23)

This optimization problem formulates the task of reaching a circle, while
trying to maintain the initial relative motor. The result of this problem can
be seen in Figure 4.

Figure 4: Two Franka Emika robots reaching a circle in the cooperative dual-
task space while trying to maintain the same relative motor.

Formulating the task in this manner circumvents the mentioned singularity
issue in the Jacobian of the cooperative pointpair, since the system will
reach circle while keeping as close as possible to their initial relative
poses.

###  IV-C Aligning Orientation Axis

When two manipulators are cooperatively manipulating an object, often it is
necessary to partially constrain their relative orientation, which is
necessary in nearly all dual-arm grasping scenarios of rigid objects. One way
to do so is enforcing collinear lines in the desired direction at the end-
effector level. This is again achieved by using the commutator product, i.e.

$\begin{split}E_{L_{1}L_{2}}&(\bm{q}_{1},\bm{q}_{2})\\\
&=M_{1}(\bm{q}_{1})L_{1}\widetilde{M}_{1}(\bm{q}_{1})\times
M_{2}(\bm{q}_{2})L_{2}\widetilde{M}_{2}(\bm{q}_{2}).\end{split}$ (24)

In Figure 5, we show the results of minimizing
$E_{L_{1}L_{2}}(\bm{q}_{1},\bm{q}_{2})$, where
$L_{1}=L_{2}=\bm{e}_{0}\wedge(\bm{e}_{0}+\bm{e}_{1}+\frac{1}{2}\bm{e}_{\infty})\wedge\bm{e}_{\infty}$,
which corresponds to aligning two lines that pass through the x-axes of the
frames at the end-effectors of the two manipulators.

Figure 5: Aligning the x-Axis.

There are, of course, other ways to achieve this behavior, e.g. by
constraining the orientation part of the relative motor. We chose to
demonstrate this method of aligning orientation, however, because it provides
a lot of flexibility, since the two lines $L_{1}$ and $L_{2}$ can be chosen
arbitrarily.

A similar instance of this problem can be defined using the absolute motor of
the CGA-CDTS. By using the absolute translator $T_{a}(\bm{q}_{1},\bm{q}_{2})$
we can move a desired line to the absolute position and require it to be
identical to the line moved by the absolute motor. This defines a desired
orientation with respect to the arbitrary axis that is defined by the line.
The formulation hence is

$\begin{split}E_{T_{L}}&(\bm{q}_{1},\bm{q}_{2})\\\
&=T_{a}(\bm{q}_{1},\bm{q}_{2})L\widetilde{T}_{a}(\bm{q}_{1},\bm{q}_{2})\times
M_{a}(\bm{q}_{1},\bm{q}_{2})L\widetilde{M}_{a}(\bm{q}_{1},\bm{q}_{2}).\end{split}$
(25)

###  IV-D Balancing a Plate

In this real-world experiment we are combining several of the previous
constraints in order to implement the task of balancing a plate. First, the
robot lifts the plate to a height of 20cm above the table, then it will try to
keep it in that position. This is formulated as constrained optimization
problem, where the objective is to stay close to the initial configuration and
the constraint is to keep the $z$-axis of absolute motor perpendicular to the
$xy$-plane.

We formulate this task as an optimal control problem

$\begin{split}\bm{u}^{*}=&\min_{\bm{u}}\Big{\|}E_{N}(\bm{x}_{N})\Big{\|}_{2}^{2}+\sum_{k=0}^{N-1}\Big{\|}E_{k}(\bm{x}_{k})\Big{\|}_{2}^{2}+\|\bm{u}_{k}\|_{R}^{2}\\\
&\text{s.t. }\bm{x}_{k+1}=f(\bm{x}_{k},\bm{u}_{k})\end{split}$ (26)

where
$\bm{x}=\begin{bmatrix}\bm{q}^{{\scriptscriptstyle\top}}_{1},\bm{q}^{{\scriptscriptstyle\top}}_{2},\bm{\dot{q}}^{{\scriptscriptstyle\top}}_{1},\bm{\dot{q}}^{{\scriptscriptstyle\top}}_{2}\end{bmatrix}^{{\scriptscriptstyle\top}}$
and
$\bm{u}=\begin{bmatrix}\bm{\ddot{q}}^{{\scriptscriptstyle\top}}_{1},\bm{\ddot{q}}^{{\scriptscriptstyle\top}}_{2}\end{bmatrix}^{{\scriptscriptstyle\top}}$.
As the state dependent cost we choose the residual multivectors of Equations
(24) and (25), where the line $L$ is chosen to be the $z$-axis, i.e.
$L=\bm{e}_{0}\wedge(\bm{e}_{0}+\bm{e}_{3}+\frac{1}{2}\bm{e}_{\infty})\wedge\bm{e}_{\infty}$.
This means that the two manipulators should simultaneously keep their relative
grasping positions on the plate and to keep the plate horizontal.

This formulation is given to a model predictive controller that uses second
order system dynamics to compute desired accelerations. Using inverse dynamics
we then compute torque commands for the control of the two manipulators. We
chose a short horizon of 10 timesteps with $\Delta t=0.01$, in order to
achieve a very reactive behavior of the controller. The model predictive
controller is then run at 100Hz and the inverse dynamics controller at 1000Hz.

Figure 6: Results of balancing a plate in different configurations.

Since there are only relative constraints in this task, the robots have a
compliant behavior when one of them or the plate is moved by hand. The other
one then adapts its configuration accordingly. We show several different
configurations in Figure 6. If no external forces are applied, the robots stay
in their current configuration. This experiment can also be seen in the
accompanying video on our website.

##  V CONCLUSION

In this article we presented an extension of the cooperative dual-task space
in conformal geometric algebra, namely the CGA-CDTS. This extension keeps all
the benefits of the original formulation that is based on dual quaternions,
but adds more tools for geometric modeling of the dual-arm tasks. After
reformulating the CDTS, we showed how the cooperative pointpair can be used to
simultaneously represent both end-effector positions and how that can be
exploited for cooperative reaching tasks and how the additional geometric
primitives facilitate the modeling of dual-arm tasks. We then demonstrated the
integration of the CGA-CDTS into an existing framework for optimal control
using geometric algebra. In future work the ideas of the CGA-CDTS could be
used to facilitate the modeling of the task spaces of robotic hands. For a
3-finger hand, for example, the cooperative pointpair could become a
cooperative circle.

## References

  * [1] S. Lee and S. Kim  “A Self-Reconfigurable Dual-Arm System”  In _1991 IEEE International Conference on Robotics and Automation Proceedings_ , 1991, pp. 164–169 vol.1  DOI: 10.1109/ROBOT.1991.131573
  * [2] Ruhizan Liza Ahmad Shauri and Kenzo Nonami  “Assembly Manipulation of Small Objects by Dual-arm Manipulator”  In _Assembly Automation_ 31.3 Emerald Group Publishing Limited, 2011, pp. 263–274  DOI: 10.1108/01445151111150604
  * [3] Petter Ögren, Christian Smith, Yiannis Karayiannidis and Danica Kragic  “A Multi Objective Control Approach to Online Dual Arm Manipulation1”  In _IFAC Proceedings Volumes_ 45.22, 10th IFAC Symposium on Robot Control, 2012, pp. 747–752  DOI: 10.3182/20120905-3-HR-2030.00032
  * [4] Diogo Almeida and Yiannis Karayiannidis  “Cooperative Manipulation and Identification of a 2-DOF Articulated Object by a Dual-Arm Robot”  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 5445–5451  DOI: 10.1109/ICRA.2018.8460511
  * [5] Jihong Zhu et al.  “Dual-Arm Robotic Manipulation of Flexible Cables”  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ Madrid: IEEE, 2018, pp. 479–484  DOI: 10.1109/IROS.2018.8593780
  * [6] Christian Smith et al.  “Dual Arm Manipulation—A Survey”  In _Robotics and Autonomous Systems_ 60.10, 2012, pp. 1340–1353  DOI: 10.1016/j.robot.2012.07.005
  * [7] Bruno Vilhena Adorno, Philippe Fraisse and Sébastien Druon  “Dual Position Control Strategies Using the Cooperative Dual Task-Space Framework”  In _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2010, pp. 3955–3960  DOI: 10.1109/IROS.2010.5650218
  * [8] Bruno Vilhena Adorno  “Robot Kinematic Modeling and Control Based on Dual Quaternion Algebra — Part I: Fundamentals.”, pp. 48
  * [9] Charles Gunn  “Geometric Algebras for Euclidean Geometry”  In _Adv. Appl. Clifford Algebras_ 27.1, 2017, pp. 185–208  DOI: 10.1007/s00006-016-0647-0
  * [10] Charles G. Gunn  “Projective Geometric Algebra: A New Framework for Doing Euclidean Geometry”, 2020  arXiv: http://arxiv.org/abs/1901.05873
  * [11] Tobias Löw and Sylvain Calinon  “Geometric Algebra for Optimal Control With Applications in Manipulation Tasks”  In _IEEE Transactions on Robotics_ , 2023, pp. 1–15  DOI: 10.1109/TRO.2023.3277282
  * [12] Emre Sariyildiz and Hakan Temeltas  “A Comparison Study of Three Screw Theory Based Kinematic Solution Methods for the Industrial Robot Manipulators”  In _2011 IEEE International Conference on Mechatronics and Automation_ Beijing, China: IEEE, 2011, pp. 52–57  DOI: 10.1109/ICMA.2011.5985630
  * [13] Mariana de Paula Assis Fonseca, Bruno Vilhena Adorno and Philippe Fraisse  “Coupled Task-Space Admittance Controller Using Dual Quaternion Logarithmic Mapping”  In _IEEE Robotics and Automation Letters_ 5.4, 2020, pp. 6057–6064  DOI: 10.1109/LRA.2020.3010458
  * [14] Riddhiman Laha et al.  “Reactive Cooperative Manipulation Based on Set Primitives and Circular Fields”  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 6577–6584  DOI: 10.1109/ICRA48506.2021.9561985
  * [15] Cristiana Miranda Farias, Yuri Gonçalves Rocha, Luis Felipe Cruz Figueredo and Mariana Costa Bernardes  “Design of Singularity-Robust and Task-Priority Primitive Controllers for Cooperative Manipulation Using Dual Quaternion Representation”  In _2017 IEEE Conference on Control Technology and Applications (CCTA)_ , 2017, pp. 740–745  DOI: 10.1109/CCTA.2017.8062550
  * [16] Jonathan Vorndamme et al.  “Integrated Bi-Manual Motion Generation and Control Shaped for Probabilistic Movement Primitives”  In _2022 IEEE-RAS 21st International Conference on Humanoid Robots (Humanoids)_ , 2022, pp. 202–209  DOI: 10.1109/Humanoids53995.2022.10000149
  * [17] Eduardo Bayro-Corrochano and Detlef Kähler  “Motor Algebra Approach for Computing the Kinematics of Robot Manipulators”  In _Journal of Robotic Systems_ 17.9, 2000, pp. 495–516  DOI: 10.1002/1097-4563(200009)17:9<495::AID-ROB4>3.0.CO;2-S
  * [18] Christian Perwass  “Geometric Algebra with Applications in Engineering”, Geometry and Computing 4  Berlin: Springer, 2009
  * [19] Eduardo Bayro-Corrochano  “Geometric Algebra Applications Vol. II: Robot Modelling and Control”  Cham: Springer International Publishing, 2020  DOI: 10.1007/978-3-030-34978-3
  * [20] Yuval Tassa, Tom Erez and Emanuel Todorov  “Synthesis and Stabilization of Complex Behaviors through Online Trajectory Optimization”  In _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2012, pp. 4906–4913  DOI: 10.1109/IROS.2012.6386025


# The strong suppression of Galactic cosmic rays reaching AU Mic b, c and Prox
Cen b

A. L. Mesquita1, D. Rodgers-Lee2,3, A. A. Vidotto1 and R. D. Kavanagh1
1Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The
Netherlands
2School of Physics, Trinity College Dublin, The University of Dublin, D02
E8C0, Ireland
3School of Cosmic Physics, Dublin Institute for Advanced Studies, 31
Fitzwilliam Place, D02 XF86, Ireland<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The propagation of Galactic cosmic rays is well understood in the context of
the solar system but is poorly studied for M dwarf systems. Quantifying the
flux of cosmic rays reaching exoplanets is important since cosmic rays are
relevant in the context of life. Here, we calculate the Galactic cosmic ray
fluxes in AU Mic and Prox Cen planetary systems. We propagate the Galactic
cosmic rays using a 1D cosmic ray transport model. We find for Prox Cen b, AU
Mic b and AU Mic c that the Galactic cosmic ray fluxes are strongly suppressed
and are lower than the fluxes reaching Earth. We include in our models, for
the first time for a star other than the Sun, the effect of radial particle
drift due to gradients and curvatures in the stellar magnetic field. For Prox
Cen we find that the inclusion of particle drift leads to less suppression of
Galactic cosmic rays fluxes than when it is excluded from the model. In the
case of AU Mic we explore two different wind environments, with a low and high
stellar wind mass-loss rate. For AU Mic, the particle drift also leads to less
suppression of the Galactic cosmic ray fluxes but it is only significant for
the high mass-loss rate scenario. However, both wind scenarios for AU Mic
suppress the Galactic cosmic rays strongly. Overall, careful modelling of
stellar winds is needed to calculate the Galactic cosmic ray fluxes reaching
exoplanets. The results found here can be used to interpret future exoplanet
atmosphere observations and in atmospheric models.

###### keywords:

methods: numerical – stars: individual: AU Microscopii, Proxima Centauri –
stars: low-mass – planetary systems – cosmic rays.

††pubyear: 2022††pagerange: The strong suppression of Galactic cosmic rays
reaching AU Mic b, c and Prox Cen b–B

##  1 Introduction

In a series of works, we have investigated the propagation of Galactic cosmic
rays in M dwarfs (Mesquita et al., 2021, 2022), with the aim to derive the
flux of cosmic rays that reach the habitable zone of these stars. We also
calculated the radiation dose reaching an exoplanet surface (Mesquita et al.,
2022). Here, we continue this series of works, focusing on two nearby M dwarfs
whose orbiting planets are amenable to atmospheric characterisation: AU Mic
and Prox Cen. M dwarfs have habitable zones that are close-in due to their low
luminosity. Because of our current observational biases towards finding close-
in planets, this means that it is easiest to observe planets orbiting in the
habitable zone of M dwarfs. However, being in the habitable zone (i.e. where
surface liquid water can exist, Kasting et al., 1993; Selsis et al., 2007) is
not thought to be the only element necessary for a planet to sustain life.
Many factors can influence exoplanet habitability (see e.g. Meadows & Barnes,
2018). One important factor for habitability is cosmic rays.

Cosmic rays are a source of ionisation that can drive the production of
prebiotic molecules in Earth-like (Airapetian et al., 2016) and Jupiter-like
atmospheres (Barth et al., 2021). Hence, cosmic rays may have been important
for the beginning of life on Earth and other planets. Fingerprint ions, such
as H3O+, H${}_{3}^{+}$ and NH${}_{4}^{+}$ (Helling & Rimmer, 2019; Barth et
al., 2021), have been identified as a good indication of ionisation in hot gas
giant atmospheres due to both Galactic and stellar cosmic rays (also known as
stellar energetic particles). Spectroscopic observations with the James Webb
Space Telescope (JWST, Gardner et al., 2006) and the Atmospheric Remote-
sensing Infrared Exoplanet Large-survey (ARIEL, Tinetti et al., 2021) of
absorption features due to these fingerprint ions may be able to constrain the
cosmic ray fluxes in exoplanet atmospheres in the future. In order to
interpret future exoplanet atmosphere observations and also to provide input
for atmospheric chemistry models, it is necessary to know the flux of Galactic
cosmic rays reaching exoplanets.

As Galactic cosmic rays propagate through stellar systems they are suppressed,
in a energy-dependent way, by the existence of a magnetised stellar wind (e.g.
Grießmeier et al., 2005; Sadovski et al., 2018; Herbst et al., 2020; Mesquita
et al., 2021; Rodgers-Lee et al., 2021b; Mesquita et al., 2022). This is
analogous to what occurs in the solar system, known as the modulation of
Galactic cosmic rays (e.g. Potgieter, 2013). The Galactic cosmic ray fluxes
reaching an exoplanet are affected by diffusion, advection, adiabatic losses
and particle drift processes that depend on the stellar wind properties (the
magnetic field and velocity). They are also affected by the Galactic cosmic
ray spectrum in the interstellar medium (ISM) and the ISM properties (velocity
and density).

Particle drifts, caused by gradients and curvatures of the heliospheric
magnetic field, have been pointed out as an important ingredient in the
description of cosmic ray transport inside the heliosphere (Jokipii et al.,
1977). Many works have included the effect of particle drifts to describe the
transport of cosmic rays in the heliosphere (Strauss et al., 2012; Strauss &
Potgieter, 2014; Potgieter, 2013; Vos & Potgieter, 2015; Potgieter & Vos,
2017; Kopp et al., 2021), in good agreement with observations at Earth.
Particle drifts were also included to study the young Sun system (Cohen et
al., 2012). However, for other stars the effect of particle drift are usually
neglected (Grießmeier et al., 2005; Herbst et al., 2020; Mesquita et al.,
2021; Rodgers-Lee et al., 2021b; Mesquita et al., 2022). One reason for this
is that to implement particle drift it is necessary to know the 3D large-scale
stellar magnetic field geometry. At this point, the 3D large scale stellar
magnetic field geometry is known for more than a hundred stars (e.g. Vidotto
et al., 2014; See et al., 2015; Klein et al., 2021a, b), albeit the stellar
wind mass-loss rates have been determined for a smaller number of stars (e.g.
Wood, 2004; Wood et al., 2021). Since different stars have different winds due
to their different (and some times extreme) magnetic fields and high rotation
rates, 3D models of the stellar wind are an important tool to study the
propagation of cosmic rays and the effects of particle drift in stellar
systems.

In this work, we investigate the well-studied M dwarfs AU Mic and Prox Cen as
they both host detected exoplanets (Anglada-Escudé et al., 2016; Plavchan et
al., 2020; Martioli et al., 2021). Prox Cen b is an Earth-size planet orbiting
at 0.048 au (Anglada-Escudé et al., 2016). AU Mic b and AU Mic c are both
Neptune-size planets at 0.064 au and 0.11 au, respectively (Plavchan et al.,
2020; Martioli et al., 2021). In particular, Prox Cen b orbits in the
habitable zone and AU Mic is a JWST target which will be observed soon using
the Near InfraRed Camera (NIRCam) to search for undetected close-in planets.
In addition, the large scale magnetic field geometry has been reconstructed
with Zeeman-Doppler Imaging (ZDI) for both Prox Cen and AU Mic (Klein et al.,
2021a, b). In terms of their stellar wind properties, Prox Cen has a mass-loss
rate constrained by Lyman-$\alpha$ observations to be $<4\times
10^{-15}\,M_{\sun}$ yr-1 (Wood et al., 2001). While mass-loss rate constraints
for AU Mic are less stringent, models predict that it ranges from 10 to 1000
times the solar wind mass-loss rate (Plavchan et al., 2020; Chiang & Fung,
2017). Kavanagh et al. (2021) modelled the winds of the Prox Cen and AU Mic
planetary systems using these observational constraints for the stellar
magnetic fields and mass-loss rates.

In this paper, we use the 3D stellar wind model results from Kavanagh et al.
(2021) to compute the propagation of Galactic cosmic rays through the AU Mic
and Prox Cen astrospheres (analogous to the heliosphere). We use our 1D model
of cosmic ray transport (Rodgers-Lee et al., 2020; Mesquita et al., 2021;
Rodgers-Lee et al., 2021b; Mesquita et al., 2022) to calculate the spectrum of
Galactic cosmic rays within the astrospheres. Here we implement in this model
the effect of particle drift for the first time in an M dwarf system. This
paper is structured as follows: our modelling framework is explained in
Section 2. We then apply our model to Prox Cen (Section 3) and AU Mic (Section
4), where we compute the cosmic ray fluxes that reach the known planets of
these systems, as well as in their habitable zones. In the case of AU Mic, we
also explore different assumptions for the stellar wind properties, in order
to investigate how a poorly constrained mass-loss rate affects our results.
Section 5 presents a discussion of our results and conclusions.

##  2 Modelling Framework

Our modelling framework consists of the following ingredients. First, we adopt
a stellar wind model to derive the stellar wind properties (described in
Section 2.1). Second, using the stellar wind properties, we model the cosmic
ray transport as a diffusive-advective process as they propagate through the
stellar wind (see Section 2.3). The Galactic cosmic rays are injected at the
boundary between the stellar wind and the ISM. We determine this outer
boundary by calculating the size of the astrosphere (see Section 2.2).

###  2.1 Stellar wind properties

In order to study the propagation of Galactic cosmic rays through an
astrosphere it is essential to know the stellar wind velocity and magnetic
field strength. For AU Mic and Prox Cen, we use the 3D magnetohydrodynamic
(MHD) wind simulations from Kavanagh et al. (2021). These simulations are
conducted using the AWSoM model (van Holst et al., 2014; Jin et al., 2017),
which was originally derived for the study of the solar wind. Additionally, to
benchmark our cosmic ray transport code with observations of cosmic rays at
the Earth, we also use a model of the solar wind, which is shown in Appendix
A.

In their study, Kavanagh et al. (2021) used reconstructed surface magnetic
field maps for Prox Cen and AU Mic that were obtained using the ZDI method
(Klein et al., 2021a, b). In the AWSoM model, the stellar wind is heated and
dissipated by Alfvén waves injected at the base of the chromosphere. To
account for the range of mass-loss rates expected for AU Mic, Kavanagh et al.
(2021) computed two different stellar wind models: one with $\dot{M}=5.5\times
10^{-13}\,M_{\sun}$ yr-1 and the other with $\dot{M}=1.2\times
10^{-11}\,M_{\sun}$ yr-1. We refer to these models as ‘AU Mic low’ and ‘AU Mic
high’, respectively. Although these two stellar wind models have the same
stellar wind magnetic field map, they have different Alfvén wave fluxes,
giving rise to different stellar wind mass-loss rates and terminal velocities.

Table 1 summarises some of the properties for the stars studied here,
including their adopted stellar wind properties. Figure 1 shows the equatorial
plane cut of the stellar wind velocity (top row) and magnetic field intensity
(bottom row) for Prox Cen, ‘AU Mic low’ and ‘AU Mic high’. The values shown
for small orbital distances ($r<100\,R_{\star}$) are from the 3D stellar wind
model of Kavanagh et al. (2021). Beyond that, the values shown are
extrapolated values (see Appendix B for details). The stellar magnetic field
geometry shapes the stellar wind velocity distribution, i.e., a smoother
magnetic field geometry results in a smoother velocity. For instance, Prox Cen
shows a complex magnetic field geometry (more than simply a dipole component)
and streams with very different speeds. AU Mic instead shows a less complex
magnetic field geometry (relatively close to a dipole) and consequently less
streams with different speeds. Compared to our model for the Sun, with
$u_{\infty}=540$ km s-1, (see Appendix A), Prox Cen (810 km s-1) and ‘AU Mic
low’ (650 km s-1) have comparable terminal wind velocities while ‘AU Mic high’
has a much faster terminal velocity of 2460 km s-1. In relation to the stellar
surface magnetic field, both Prox Cen and AU Mic have a much higher magnetic
field intensity in comparison to our model for the Sun, $B\sim 2$ G (see
Appendix A).

Table 1: Properties of the Prox Cen and AU Mic stellar systems. The quantities
without reference are calculated in this work.

Physical parameter | Prox Cen | AU Mic | Unit
---|---|---|---
|  | ‘AU Mic low’ | ‘AU Mic high’ |
Stellar mass ($M_{\star}$) | 0.12a | 0.50b | $M_{\sun}$
Stellar radius ($R_{\star}$) | 0.14a | 0.75b | $R_{\sun}$
Distance ($d$) | 1.30a | 9.79b | pc
Habitable zone | 0.03 – 0.09c | 0.26 – 0.66c | au
Average stellar surface magnetic field ($B_{\star}$) | 200d | 475e | G
ISM/star relative velocity ($v_{\rm ISM}$) | 25f | 19 | km s-1
ISM ram pressure ($P_{\rm ISM}$) | $2.5\times 10^{-12}$ | $1.5\times 10^{-12}$ | dyn cm-2
Stellar wind mass-loss rateg ($\dot{M}$) | $5.0\times 10^{-15}$ | $5.5\times 10^{-13}$ | $1.2\times 10^{-11}$ |  $M_{\sun}$ yr-1
Stellar wind ram pressure ($P_{\rm ram}$) | $3.3\times 10^{-6}$ | $1.2\times 10^{-5}$ | $4.6\times 10^{-4}$ | dyn cm-2
Terminal velocityg ($u_{\infty}$) | 810 | 650 | 2460 | km s-1
Astrospheric size ($R_{\rm ast}$) | 75 | 980 | 6140 | au

aAnglada-Escudé et al. (2016); bPlavchan et al. (2020); ccalculated using the
prescription of Kopparapu et al. (2014); dKlein et al. (2021a); eKlein et al.
(2021b); fWood et al. (2001); gKavanagh et al. (2021).

Figure 1: Equatorial plane cut of the stellar wind properties, from left to
right, for Prox Cen, ‘AU Mic low’ and ‘AU Mic high’. The top row shows the
stellar wind velocity and the bottom row shows the magnitude of the total
magnetic field for the full radial extent of each star’s astrosphere. The
values shown for small orbital distances are from 3D stellar wind models
(Kavanagh et al., 2021) and the values shown at large distances are an
extrapolation (see Appendix B for details).

###  2.2 The astrosphere size

The astrosphere is the ‘bubble’ region around a star dominated by its stellar
wind. The astrosphere size can be very important for the propagation of
Galactic cosmic rays because outside this region cosmic rays are not
suppressed. The size of the astrosphere can be calculated as

$R_{\text{ast}}=\sqrt{\frac{P_{\text{ram}}(R)}{P_{\text{ISM}}}}R,$ (1)

where $P_{\text{ram}}$ is the ram pressure of the stellar wind expressed as

$P_{\text{ram}}=\rho u^{2}=\frac{\dot{M}u}{4\pi R^{2}},$ (2)

where $u$ is the wind velocity, $R$ is a given reference distance beyond which
the wind has reached its terminal velocity. The ISM ram pressure,
$P_{\text{ISM}}$, can be expressed as

$P_{\text{ISM}}=m_{p}n_{\rm ISM}\nu_{\rm ISM}^{2},$ (3)

where $m_{p}$ is the proton mass, $n_{\text{ISM}}$ is the total ISM number
density of hydrogen and $\nu_{\text{ISM}}$ is the ISM velocity as observed by
the star.

To derive the total ISM number density of hydrogen ($n_{\rm ISM}$), we use the
same ISM neutral hydrogen density ($n_{\rm n}=0.14\,\text{cm}^{-3}$) and
ionised hydrogen density ($n_{\rm i}=0.1\,\text{cm}^{-3}$) as given by Model
10 of Wood et al. (2000) which is a good fit for nearby systems in the local
ISM. Assuming that the ISM cloud along the line of sight towards AU Mic is
most consistent with the Mic cloud, we derive the heliocentric flow vector in
the direction of AU Mic and the radial velocity of the Mic cloud using the ISM
Kinematic Calculator (Redfield & Linsky, 2008). The heliocentric velocity of
AU Mic is calculated using its proper motion (Gaia Collaboration, 2020) and
its radial velocity (Fouqué et al., 2018). From these assumptions, we
calculate the ISM velocity and pressure as seen by AU Mic to be $v_{\rm
ISM}=19$ km s-1 and $P_{\rm ISM}=1.5\times 10^{-12}$dyn cm-2. If instead we
assume that the ISM cloud along the line of sight towards AU Mic is consistent
with the Local Interstellar Cloud (LIC), the ISM velocity does not vary
significantly ($v_{\rm ISM}=19.4$ km s-1 ). For Prox Cen, we use the ISM
velocity as given by Wood et al. (2001) and shown in Table 1. The astrospheric
size of Prox Cen and AU Mic (for the two stellar wind scenarios) are given in
Table 1. The astrospheric size is proportional to the stellar wind mass-loss
rate. Compared with the heliosphere, where $R_{\rm ast}=122$ au (as observed
by Voyager 1, Stone et al., 2013, 2019), Prox Cen’s astrosphere is 40% smaller
while $R_{\rm ast}$ for ‘AU Mic low’ is 8 times larger and for ‘AU Mic high’
it is 50 times larger.

###  2.3 Cosmic ray transport equation

When interacting with magnetised stellar winds, Galactic cosmic rays can
experience global changes in their intensity and energy, referred to as
modulation. The cosmic ray propagation through a magnetised stellar wind can
be described by the transport equation of Parker (1965). We use a 1D cosmic
ray transport model to calculate the cosmic rays fluxes within the M dwarf
systems. The model used here is based on the model described in Rodgers-Lee et
al. (2020) and was already used for M dwarfs (Mesquita et al., 2021, 2022).
Here, we add for the first time, for a star other than the Sun, the effect of
particle drift.

We numerically solve the time-dependent, spherically symmetric transport
equation for cosmic rays, as given by

$\frac{\partial f}{\partial t}=\nabla\cdot(\mathbf{\kappa}\nabla
f)-(\mathbf{u}+\langle\mathsf{\mathbf{v}}_{d}\rangle)\cdot\nabla
f+\frac{1}{3}(\nabla\cdot\mathbf{u})\frac{\partial f}{\partial\ln p},$ (4)

where $f$ and $p$ are the cosmic ray phase space density and momentum,
respectively. The diffusion of cosmic rays is represented by the first term on
the right-hand side of Equation 4 and depends on the diffusion coefficient,
$\kappa$. The stellar wind velocity, $\mathbf{u}$, and the particle drift
velocity, $\langle\mathsf{\mathbf{v}}_{d}\rangle$, appear in the second term
and represent advection of the cosmic rays. The last term represents momentum
advection, also known as adiabatic losses, which moves the cosmic rays to
lower energies. More details of the model can be found in Rodgers-Lee et al.
(2020). Our results in Sections 3 and 4 present the differential intensity of
cosmic rays, $j$, rather than the phase space density, $f$, where
$j(T)=p^{2}f(p)$. $T$ is the cosmic ray kinetic energy.

We use logarithmically spaced spatial and momentum grids. The inner spatial
boundary is set as 0.01 au and the outer boundary is set as the astrosphere
size, $R_{\rm ast}$, for each star (see Table 1). The spatial grid includes
the orbits of the known exoplanets and the habitable zones. The spatial grid
has 80 grid zones for Prox Cen, 100 grid zones for the Sun and ‘AU Mic low’
and 120 grid zones for ‘AU Mic high’. The spatial resolution varies for each
system because they have different astrosphere sizes. The inner spatial
boundary condition is reflective and the outer spatial boundary condition is
fixed as the Local Interstellar Spectrum (LIS) for Galactic cosmic rays. The
LIS that we adopt is given by the model fit to the Voyager 1 observations
using Eq. (1) of Vos & Potgieter (2015) 111As discussed in Mesquita et al.
(2021, 2022), this is a good assumption since according to $\gamma$-ray
observations the inferred cosmic ray spectrum within a region of 1 kpc in the
local Galaxy is consistent with Voyager measurements in the local ISM (Neronov
et al., 2017).. The LIS is considered to be constant as a function of time in
our simulations. The momentum grid ranges from $p_{\rm min}=0.15$ GeV$/c$ to
$p_{\rm max}=100$ GeV$/c$ with 60 momentum bins. The limits of the momentum
grid are chosen due to the fact that cosmic rays with higher energies than the
upper limit are quite uncommon and not important for planetary atmosphere
chemistry (Rimmer & Helling, 2013), while cosmic rays with lower energy than
the lower limit are not expected to reach the planetary surface (Atri, 2017).
Low-energy particles, however, can deposit energy in the planetary atmosphere
(Rodgers-Lee et al., 2020) and may be relevant for chemical models (e.g.
Grießmeier et al., 2015; Barth et al., 2021). Both, inner and outer, momentum
boundary conditions are outflow boundary conditions.

####  2.3.1 The diffusion coefficient

The diffusion coefficient, $\kappa$, in Equation 4, depends on the level of
turbulence in the stellar magnetic field and the cosmic ray momentum. From
quasi-linear theory (Jokipii, 1966; Schlickeiser, 1989), the diffusion
coefficient can be denoted as

$\frac{\kappa(r,p)}{\beta
c}=\eta_{0}\left(\frac{p}{p_{0}}\right)^{1-\gamma}r_{\text{L}},$ (5)

where $r_{\text{L}}$ is the cosmic ray Larmor radius, $p_{0}=3\,$GeV$/c$ and
$\beta=v/c$ is the ratio between the cosmic ray’s velocity and the speed of
light. $\eta_{0}$ represents the level of turbulence present in the system. In
our simulations we set $\eta_{0}=1$ to be consistent with other works
(Rodgers-Lee et al., 2020; Mesquita et al., 2021, 2022). $\gamma$ determines
the power-law dependency with momentum. The type of turbulence can be
described by different prescriptions, such as, Bohm-type ($\gamma=1$),
Kolmogorov-type ($\gamma=5/3$) and MHD-driven ($\gamma=3/2$) turbulence. We
note that the Galactic cosmic ray spectrum calculated at any given distance is
affected by the type of turbulence adopted. The type of turbulence present in
M dwarf systems is still unknown and here we simply adopt $\gamma=1$, i.e.
Bohm-type turbulence. This has been used frequently in many works (Svensmark,
2006; Cohen et al., 2012; Rodgers-Lee et al., 2020; Mesquita et al., 2021;
Rodgers-Lee et al., 2021b; Mesquita et al., 2022) and fits the present-day
observations at Earth quite well (Rodgers-Lee et al., 2020).

The magnetic field and velocity are important wind properties in the context
of Galactic cosmic ray propagation. For instance, the stronger the magnetic
field, the more the Galactic cosmic ray spectrum is suppressed (for fixed
values of $\gamma$ and $\eta_{0}$). This is because a strong magnetic field
implies smaller diffusion coefficients. This means that generally advective
processes become more important, resulting in more modulation of cosmic rays.
The same is valid for the wind velocity, however to a lesser extent (Mesquita
et al., 2022), i.e. a stronger stellar wind velocity results in more advection
which suppresses the cosmic ray fluxes. The stellar wind density does not
attenuate the cosmic rays, as the density is very low. However, it influences
the size of the astrosphere (see Equation 1). In cases where advection is the
dominant physical process, the size of the astrosphere can affect the Galactic
cosmic ray fluxes (Rodgers-Lee et al., 2021b).

####  2.3.2 The particle drift

The averaged particle drift velocity, $\langle\mathsf{\mathbf{v}}_{d}\rangle$,
is caused by gradients and curvatures in the stellar magnetic field. The drift
velocity can be expressed as a function of the magnetic field as (Jokipii et
al., 1993)

$\langle\mathsf{\mathbf{v}}_{d}\rangle=\frac{p\beta
c^{2}}{3e}\mathbf{\nabla}\times\frac{\mathbf{B}}{B^{2}},$ (6)

where $B$ is the stellar wind magnetic field. The drift velocity also depends
on the cosmic ray momentum $p$ and velocity $v$ (through $\beta$).

In our 1D propagation model, we only include the radial component of the drift
velocity, which can be expressed in spherical coordinates as

$\langle\mathsf{\mathbf{v}}_{d}\rangle=\frac{p\beta
c^{2}}{3e}\frac{1}{r\sin\theta}\left(\underbrace{\frac{\partial}{\partial\theta}\left(\sin\theta\frac{B_{\phi}}{B^{2}}\right)}_{\rm
i}-\underbrace{\frac{\partial}{\partial\phi}\left(\frac{B_{\theta}}{B^{2}}\right)}_{\rm
ii}\right)\mathbf{e}_{r},$ (7)

where $\theta$ is the polar angle ($\theta=90^{\circ}$ at the equator) and
$\phi$ is the azimuthal angle.

Figure 2 shows the absolute value of the drift velocity components that we
calculate for Prox Cen for cosmic rays with $p=\,$1 GeV/c. The top panel shows
term i of Equation 7 as a function of radial distance and polar angle,
$\theta$, for the value of $\phi$ where the stellar wind velocity have its
strongest intensity ($\phi=97^{\circ}$). The bottom panel shows term ii of
Equation 7 as a function of radial distance and the azimuthal angle, $\phi$,
for the equatorial plane ($\theta=90^{\circ}$). The dotted areas represent
regions with negative velocities. The momentum dependency of the drift
velocity is a multiplicative factor, i.e., the drift velocity is 100 times
higher for 100 GeV/c momentum cosmic rays than for 1 GeV/c cosmic rays. Both
terms of the drift are spatially non-linear and vary as a function of $r$,
$\theta$ and $\phi$ significantly, as shown in Figure 2. In the top panel of
Figure 2 we did not include $\theta=0^{\circ}$ because $1/\sin\theta$ in
Equation 7 implies that the drift value becomes infinite at this point. In our
simulations we only include the drift velocity as a function of $r$ and $p$
for a single value of $\theta$ and $\phi$ where i and ii are combined. We will
discuss the effect of the 2D/3D spatial variation of the drift velocity in
Section 5.

Figure 2: Components of the radial drift velocity for Prox Cen for cosmic rays
with 1 GeV/c momentum. Top panel: term i of Equation 7 as a function of $r$
and $\theta$ for $\phi$ direction where the stellar wind is fastest (see
Figure 1 top row). Bottom panel: term ii of Equation 7 as a function of $r$
and $\phi$ for the equatorial plane. In both plots the dotted regions
represent regions with negative velocities.

The effects of $\langle\mathsf{\mathbf{v}}_{d}\rangle$ have been neglected in
past studies of M dwarf systems. However, if
$\langle\mathsf{\mathbf{v}}_{d}\rangle$ is larger or comparable to the stellar
wind velocity, $u$, it can affect the modulation of cosmic rays. The particle
drift acts as an additional advective term. Hence, when the particle drift
velocity is positive, i.e. radially outward pointing, then the Galactic cosmic
ray fluxes are reduced. We will discuss this in more detail in Sections 3 and
4.

##  3 Prox Cen: the Galactic cosmic ray fluxes including particle drift
effects

Prox Cen is the closest star to the solar system and has a terrestrial planet
in its habitable zone. Here, we investigate the propagation of Galactic cosmic
rays within its astrosphere. Since we use a 1D cosmic ray transport model, we
select particular cuts of the 3D wind simulations to calculate the drift
velocity at any given orbital distance. We investigate the Galactic cosmic ray
fluxes at the equatorial plane ($\theta=90^{\circ}$) because this allows us to
investigate the effect of Galactic cosmic rays at the planet’s equatorial
orbit. Since the stellar wind velocity also varies with $\phi$, we further fix
the value of $\phi=97^{\circ}$ for Prox Cen. We note that at the chosen
$\theta$ and $\phi$ the stellar wind has already reached its terminal velocity
within the 3D stellar wind simulation domain.

Figure 3 shows the intensity of cosmic rays as a function of kinetic energy
for Prox Cen when particle drift is considered (dashed lines) and when it is
neglected (solid lines) at different distances. The shaded areas represents
the Galactic cosmic ray flux in the habitable zone when particle drift is
considered (dark salmon) and when it is neglected (light salmon). The black
line is the LIS. When particle drift is included, Prox Cen shows a
considerably higher flux of cosmic rays for particles with kinetic energy
below 1 GeV. At the orbit of Prox Cen b, the intensity of 0.1 GeV energy
cosmic rays is 7 times higher when drift is considered (dashed violet line) in
comparison with results when drift is neglected (solid violet line). A similar
trend is observed at 1 au (green lines) and in the habitable zone (shaded
areas).

Figure 3: Differential intensity of Galactic cosmic rays as a function of
cosmic ray kinetic energy for Prox Cen. The solid lines represents the
Galactic cosmic ray flux without considering the particle drift and the dashed
lines when the drift is considered. The black solid line is the LIS. The green
lines are the Galactic cosmic ray flux at 1 au and the violet lines are the
Galactic cosmic ray flux at Prox Cen b. The shaded areas are the Galactic
cosmic ray flux in the habitable zone, both accounting for (dark salmon) and
neglecting (light salmon) particle drift. In the case of Prox Cen, the
inclusion of particle drift leads to Galactic cosmic ray fluxes that are an
order of magnitude higher for cosmic ray energies $\lesssim 0.1$ GeV.

The effects we see with the inclusion of particle drifts can be understood by
comparing the drift velocity with the stellar wind velocity. Figure 4 shows
the stellar wind velocity and the drift velocity profiles (for
$\theta=90^{\circ}$, $\phi=97^{\circ}$) for Prox Cen. The green triangles are
$|\langle\mathsf{v}_{d}\rangle|$ for cosmic rays with $p=1$ GeV/c, while the
blue triangles are for $p=0.15$ GeV/c. The stellar wind velocity is lower than
$|\langle\mathsf{v}_{d}\rangle|$ for $r>0.3\,$au (green triangles) and, as a
consequence, the drift velocity affects the modulation of cosmic rays. The
Galactic cosmic ray flux within the Prox Cen system is higher when particle
drifts are included because the drift velocity is large and negative in value
(open triangles in Figure 4). If, on the other hand, the drift velocity were
positive, the cosmic ray flux would be reduced.

Figure 4: Stellar wind and radial drift velocity as a function of distance for
Prox Cen at $\theta=90^{\circ}$, $\phi=97^{\circ}$. The magenta line is the
stellar wind velocity from Kavanagh et al. (2021). The green and blue
triangles are the absolute value of the radial drift velocity for cosmic rays
with 1 GeV/c and 0.15 GeV/c momentum, respectively. Open triangles represent
negative values and filled triangles represent positive values. Cosmic rays
with higher momenta will have larger drift velocity values and vice-versa.

Our results likely depend on the direction we choose to study ($\theta$ and
$\phi$ values). In the case of Prox Cen for instance, it is possible that for
azimuthal angles where the wind velocity is lower the particle drift may play
an even more important role in modulating the cosmic rays. In Figure 2 we can
see that each term of the drift velocity varies significantly with $\theta$
and $\phi$. A 2D or 3D cosmic ray simulation would be necessary to investigate
the cosmic ray propagation in other directions.

##  4 AU Mic: the effect of an unconstrained mass-loss rate on Galactic cosmic
ray fluxes

AU Mic is a young M dwarf stellar system that hosts two Neptune-size
exoplanets. Both exoplanets have close-in orbits and are not within the
habitable zone. As discussed in Section 2.1, we investigate Galactic cosmic
ray propagation for the two stellar wind scenarios explored by Kavanagh et al.
(2021). Although their magnetic field profiles are similar, the mass-loss rate
and terminal velocity of ‘AU Mic high’ are, respectively, 21 and 4 times
larger than the values for ‘AU Mic low’. Because of their different wind
properties, the calculated sizes of the astrosphere differ for these two wind
models. Similar to Prox Cen, we investigate the flux of Galactic cosmic rays
for AU Mic in the equatorial plane and we set $\phi=133^{\circ}$ for ‘AU Mic
low’ and $\phi=126^{\circ}$ for ‘AU Mic high’.

Figure 5 shows the cosmic ray fluxes as a function of kinetic energy for the
two wind scenarios of AU Mic. Here, in this figure, we only show the results
of our calculations that consider particle drift. The orange shaded areas
represents the cosmic ray fluxes in the habitable zone, where the dark orange
is for ‘AU Mic low’ and the light orange for ‘AU Mic high’. Blue, purple and
salmon lines are the cosmic ray fluxes at 1 au, AU Mic b and AU Mic c. In
general, the Galactic cosmic ray fluxes for ‘AU Mic high’ are much reduced in
comparison to ‘AU Mic low’ (see for instance blue lines in Figure 5). The
difference is particularly sharp for cosmic rays with energy lower than 1 GeV.
Within each wind prescription, the intensity of cosmic rays at AU Mic b and AU
Mic c is similar since the planets have close orbits.

In the case of AU Mic, the unconstrained mass-loss rate of the stellar wind
and the different terminal velocities from each model strongly affect the
Galactic cosmic ray intensities within the system as shown in Figure 5. This
is opposite to what was derived for GJ 436 by Mesquita et al. (2021), in which
the wind properties did not strongly affect the cosmic ray fluxes. The reason
for this, is that AU Mic is dominated by advective processes (with or without
the inclusion of particle drift), while GJ 436 is dominated by diffusive
processes. When advection dominates it results in a stronger suppression of
cosmic rays, while when diffusion dominates it results in little (or no)
modulation. The diffusive/advective processes are correlated with the stellar
wind velocity and magnetic field (see Section 2.3.1 and the
diffusive/advective time-scales in Section 3.2 of Mesquita et al. 2021). In
addition, when advection dominates the size of the astrosphere can affect the
suppression of cosmic rays (Rodgers-Lee et al., 2021b).

Figure 5: Differential intensity of Galactic cosmic rays as a function of
cosmic ray kinetic energy for ‘AU Mic low’ (solid lines) and ‘AU Mic high’
(dash dotted lines). Blue, purple and salmon lines are the differential
intensity of cosmic rays at 1 au, AU Mic b and AU Mic c, respectively. The
orange shaded areas represent the fluxes in the habitable zone, where the dark
orange is for ‘AU Mic low’ and the light orange for ‘AU Mic high’. The stellar
wind properties for the ‘AU Mic high’ case are more efficient at suppressing
the Galactic cosmic ray fluxes.

In the case of ‘AU Mic low’, the drift velocity profile for cosmic rays with 1
GeV/c momentum is negative but almost negligible when compared with the
stellar wind velocity profile and, as a consequence, the cosmic ray fluxes are
not influenced by the inclusion of drift. For ‘AU Mic high’, the drift
velocity for cosmic rays with 1 GeV/c momentum is negligible ($\lesssim
30\,\rm km~{}s^{-1}$) for $r<1\,$au. For $1\,\rm au<r<10\,\rm au$, it is
positive and relatively slow, while for $r>10\,$au
$|\langle\mathsf{v}_{d}\rangle|$ is negative and fast (but slower than the
drift velocity intensities observed for Prox Cen). Thus, the Galactic cosmic
ray intensities are slightly affected by the inclusion of particle drift for
close-in distances ($r<1\,$au). For large distances ($r>1\,$au) neglecting the
drift can underestimate the cosmic ray fluxes. For instance, we find at 744 au
for 1 GeV energy cosmic rays, the intensity of Galactic cosmic rays when drift
is considered is two orders of magnitude higher when compared with no drift.

Overall, particle drift is an important ingredient in cosmic ray modulation,
particularly when $|\langle\mathsf{v}_{d}\rangle|>u$. Here, we fixed $\phi$
and $\theta$ values as previously discussed. However, we can expect that the
presence of particle drift will affect the cosmic ray modulation of Prox Cen
and ‘AU Mic high’ at any given location of the equatorial plane. This is due
to the fact that particle drift changes the modulation of cosmic ray even in
the direction where the wind have the strongest velocity.

##  5 Discussion & Conclusions

In this paper, we investigated the intensity of Galactic cosmic rays within
the astrospheres of two planetary systems: AU Mic and Prox Cen. AU Mic hosts
two known close-in Neptune-size exoplanets. Prox Cen is our closest star and
has a known Earth-size exoplanet orbiting within the habitable zone. We used a
1D cosmic ray transport model which includes the effect of particle drift to
calculate the cosmic ray fluxes reaching the exoplanets and the habitable
zones. The radial particle drift velocity is calculated using 3D MHD stellar
wind models for Prox Cen, AU Mic (with a ‘low’ and ‘high’ mass loss rate) and
the Sun for comparison (Kavanagh et al., 2021, see also Appendix A). This is
the first time that radial drift velocities have been quantified for M dwarf
systems.

Overall, the inclusion of particle drift velocities, due to gradients and
curvatures of the stellar magnetic field, has a significant effect on the
propagation of Galactic cosmic rays. In general, a negative drift velocity
results in a larger intensity of cosmic rays and a positive drift velocity
results in a reduction of cosmic rays. This is because the particle drift acts
simply as an extra advective process in the transport model. In the two
systems studied here, the drift velocity is negative at the cut where we
investigated the flux of Galactic cosmic rays. Prox Cen and ‘AU Mic high’ both
have a larger negative drift velocity while for ‘AU Mic low’ the drift
velocity is rather slow for cosmic rays with 1 GeV/c momentum. Consequently,
the modulation of cosmic rays is significantly affected by the inclusion of
drift for Prox Cen and ‘AU Mic high’ but only slightly affected in the case of
‘AU Mic low’.

For Prox Cen, at 1 au (green dashed line in Figure 3) the Galactic cosmic ray
fluxes are larger in comparison to Earth. This is because Prox Cen’s magnetic
field profile is weaker at larger distances (e.g., from 1 au to the
astrosphere edge) in comparison to the solar system and results in less
suppression of Galactic cosmic rays for the same distance. In contrast, both
AU Mic scenarios (low and high) have a much stronger stellar magnetic field
profile in comparison to the solar system, resulting in an intense suppression
of Galactic cosmic rays at the same distance.

Herbst et al. (2020) have also investigated the Galactic cosmic ray fluxes at
Prox Cen b, albeit without considering the effects of particle drift. They
found that cosmic rays were not effectively reduced in comparison to the LIS.
In contrast, our results show lower cosmic ray fluxes at Prox Cen b. For
instance, compared to Herbst et al. (2020), the cosmic ray fluxes at Prox Cen
b with energy 1 GeV is 7 times smaller in our simulations when particle drift
is not included and 5 times smaller when particle drift is included. The
differences in our results and Herbst et al. (2020) is likely to be related to
the different wind and ISM properties used and the cosmic ray transport
properties. They used a surface magnetic field strength of $\sim 600\,$G
measured from Zeeman broadening (Reiners & Basri, 2008) and a wind velocity of
$1500\,\rm km~{}s^{-1}$ based on MHD simulations (Garraffo et al., 2016). Both
quantities are higher than the values we used. We also calculate different
astrosphere sizes, 122 au (Herbst et al., 2020) versus 75 au in our work, as a
result of our different wind and ISM properties. If the drift velocities for
Prox Cen used in our work were included in Herbst et al. (2020) work, they
would probably find an even larger flux of Galactic cosmic rays than the
values they calculate for Prox Cen b.

In the case of AU Mic, we showed that our lack of a strong constraint for the
wind properties (mass-loss rate and wind velocity) strongly affects the
intensity of Galactic cosmic rays calculated for the system and the size of
the astrosphere (980–6140 au). This is in contrast with the results found by
Mesquita et al. (2021) for GJ 436, in which the lack of knowledge of $u$ and
$\dot{M}$ did not strongly affect the flux of Galactic cosmic rays in the
system. This is related to the physical process that is more dominant
(Rodgers-Lee et al., 2021b). While GJ 436 is dominated by diffusion which
leads to little (or no) modulation of cosmic rays, AU Mic is dominated by
advection processes which strongly suppress cosmic rays. Interestingly, AU Mic
is dominated by advective processes even without the inclusion of the particle
drift. Additionally, because ‘AU Mic high’ has a larger astrosphere size the
Galactic cosmic rays have further to travel through a region of the stellar
wind dominated by advection in comparison with ‘AU Mic low’.

In addition to drift effects, the stellar magnetic field and the wind velocity
are important ingredients to effectively modulate the Galactic cosmic ray
fluxes throughout the astrosphere. A strong stellar magnetic field, leading to
smaller diffusion coefficients, combined with a fast stellar wind velocity is
effective at suppressing the propagation of Galactic cosmic rays inside the
astrosphere. The magnetic field intensity can change if the star has an
activity cycle, similar to the solar cycle (see e.g. Hathaway, 2010). When the
Sun is at activity minimum the intensity of Galactic cosmic rays is higher
than at activity maximum (Potgieter, 2013), with more than an order of
magnitude difference in flux, depending on the cosmic ray energy (Vos &
Potgieter, 2015). Prox Cen, similarly to the Sun, has been suggested to have a
7 year activity cycle (e.g. Yadav et al., 2016). The magnetic map adopted in
our wind model was derived near its activity maximum (Klein et al., 2021a).
This could imply that during its activity minimum Prox Cen would likely
experience a higher flux of Galactic cosmic rays than the results presented
here. In addition, the particle drift will also vary with activity cycle since
it depends on the magnetic field geometry.

Another important aspect to consider is that stars themselves can accelerate
stellar cosmic rays during events such as flares and coronal mass ejections
(Rodgers-Lee et al., 2021a). Prox Cen and AU Mic are both active stars with
strong surface magnetic fields (Klein et al., 2021a, b) and frequent flaring
activity (Gilbert et al., 2021; Gilbert et al., 2022) and, as such, they
should be efficient at accelerating stellar cosmic rays. As a result, for
these stars, stellar cosmic ray fluxes will dominate over Galactic cosmic rays
up to a certain energy. Some works have explored the effect of such events for
the exoplanets in the AU Mic and Prox Cen systems. Scheucher et al. (2020)
showed that a strong stellar cosmic ray event would be able to heat the
otherwise cold planet Prox Cen b, playing an important role in the planet’s
habitability. Recently, Carolan et al. (2020) investigated the effects of
strong stellar winds in the atmospheric evaporation of AU Mic b. They
concluded that, even when atmospheric erosion by the stellar wind is not
significant, the geometry and ionisation of the escaping atmosphere can change
substantially, impacting the interpretation/prediction of spectroscopic
transit observations. In addition to strong winds, coronal mass ejections and
flares can also affect planetary atmospheres, by increasing atmospheric
erosion (Khodachenko et al., 2007; Hazra et al., 2020; Hazra et al., 2022).
Although these processes are episodic, they are expected to be more frequent
and more energetic in young and/or active stars (Aarnio et al., 2014), such as
AU Mic.

To better account for the full 3D nature of planetary systems, 2D/3D cosmic
ray transport models should be used, similar to models implemented for the
solar system. Although the radial drift velocity varies with the polar and
azimuthal angles, which we have neglected, it is possible with 2D/3D
simulations that diffusion could smooth these variations, resulting in a more
continuous distribution of cosmic rays. In this case, the differential
intensity might not be too different than the ones we calculated with our
model. The polar and azimuthal drift velocities may also act in a similar way.

Finally, the Galactic cosmic rays fluxes at the exoplanets’ orbits can be used
to further investigate the intensity of Galactic cosmic rays in planetary
atmospheres and in atmospheric chemistry models. Our models can also be used
to understand future observations of exoplanet atmospheres with instruments
such as JWST and ARIEL.

## Acknowledgements

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 817540, ASTROFLOW). DRL would like to acknowledge that
this publication has emanated from research conducted with the financial
support of Science Foundation Ireland under Grant number 21/PATH-S/9339. The
authors wish to acknowledge the SFI/HEA Irish Centre for High-End Computing
(ICHEC) for the provision of computational facilities and support.

## Data Availability

The data described in this article will be shared on reasonable request to the
corresponding author.

## References

  * Aarnio et al. (2014) Aarnio A. N., Stassun K. G., Matt S. P., 2014, Proceedings of the International Astronomical Union, 8, 318
  * Airapetian et al. (2016) Airapetian V. S., Glocer A., Gronoff G., Hébrard E., Danchi W., 2016, Nature Geoscience, 9, 452
  * Anglada-Escudé et al. (2016) Anglada-Escudé G., et al., 2016, Nature, 536, 437
  * Atri (2017) Atri D., 2017, MNRAS, 465, L34
  * Barth et al. (2021) Barth P., et al., 2021, MNRAS, 502, 6201
  * Boro Saikia et al. (2020) Boro Saikia S., Jin M., Johnstone C. P., Lüftinger T., Güdel M., Airapetian V. S., Kislyakova K. G., Folsom C. P., 2020, A&A, 635, A178
  * Carolan et al. (2020) Carolan S., Vidotto A. A., Plavchan P., Villarreal D’Angelo C., Hazra G., 2020, MNRAS, 498, L53
  * Chiang & Fung (2017) Chiang E., Fung J., 2017, ApJ, 848, 4
  * Cohen (2011) Cohen O., 2011, MNRAS, 417, 2592
  * Cohen et al. (2012) Cohen O., Drake J. J., Kóta J., 2012, ApJ, 760, 85
  * Fouqué et al. (2018) Fouqué P., et al., 2018, MNRAS, 475, 1960
  * Gaia Collaboration (2020) Gaia Collaboration 2020, VizieR Online Data Catalog, p. I/350
  * Gardner et al. (2006) Gardner J. P., et al., 2006, Space Sci. Rev., 123, 485
  * Garraffo et al. (2016) Garraffo C., Drake J. J., Cohen O., 2016, ApJ, 833, L4
  * Gilbert et al. (2021) Gilbert E. A., Barclay T., Kruse E., Quintana E. V., Walkowicz L. M., 2021, Frontiers in Astronomy and Space Sciences, 8, 190
  * Gilbert et al. (2022) Gilbert E. A., et al., 2022, AJ, 163, 147
  * Grießmeier et al. (2005) Grießmeier J. M., Stadelmann A., Motschmann U., Belisheva N. K., Lammer H., Biernat H. K., 2005, Astrobiology, 5, 587
  * Grießmeier et al. (2015) Grießmeier J. M., Tabataba-Vakili F., Stadelmann A., Grenfell J. L., Atri D., 2015, A&A, 581, A44
  * Hathaway (2010) Hathaway D. H., 2010, Living Reviews in Solar Physics, 7, 1
  * Hazra et al. (2020) Hazra G., Vidotto A. A., D’Angelo C. V., 2020, MNRAS, 496, 4017
  * Hazra et al. (2022) Hazra G., Vidotto A. A., Carolan S., Villarreal D’Angelo C., Manchester W., 2022, MNRAS, 509, 5858
  * Helling & Rimmer (2019) Helling C., Rimmer P. B., 2019, Philosophical Transactions of the Royal Society of London Series A, 377, 20180398
  * Herbst et al. (2020) Herbst K., et al., 2020, ApJ, 897, L27
  * Jin et al. (2017) Jin M., Manchester W. B., van der Holst B., Sokolov I., Tóth G., Vourlidas A., de Koning C. A., Gombosi T. I., 2017, ApJ, 834, 172
  * Jokipii (1966) Jokipii J. R., 1966, ApJ, 146, 480
  * Jokipii et al. (1977) Jokipii J. R., Levy E. H., Hubbard W. B., 1977, ApJ, 213, 861
  * Jokipii et al. (1993) Jokipii J. R., Kota J., Merenyi E., 1993, ApJ, 405, 782
  * Kasting et al. (1993) Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, Icarus, 101, 108
  * Kavanagh et al. (2021) Kavanagh R. D., Vidotto A. A., Klein B., Jardine M. M., Donati J.-F., Ó Fionnagáin D., 2021, MNRAS, 504, 1511
  * Khodachenko et al. (2007) Khodachenko M. L., et al., 2007, Astrobiology, 7, 167
  * Klein et al. (2021a) Klein B., Donati J.-F., Hébrard É. M., Zaire B., Folsom C. P., Morin J., Delfosse X., Bonfils X., 2021a, MNRAS, 500, 1844
  * Klein et al. (2021b) Klein B., et al., 2021b, MNRAS, 502, 188
  * Kopp et al. (2021) Kopp A., Raath J. L., Fichtner H., Potgieter M. S., Ferreira S. E. S., Heber B., 2021, ApJ, 922, 124
  * Kopparapu et al. (2014) Kopparapu R. K., Ramirez R. M., SchottelKotte J., Kasting J. F., Domagal-Goldman S., Eymet V., 2014, ApJ, 787, L29
  * Martioli et al. (2021) Martioli E., Hébrard G., Correia A. C. M., Laskar J., Lecavelier des Etangs A., 2021, A&A, 649, A177
  * McComas et al. (2000) McComas D. J., et al., 2000, J. Geophys. Res., 105, 10419
  * Meadows & Barnes (2018) Meadows V. S., Barnes R. K., 2018, in Deeg H. J., Belmonte J. A., eds, , Handbook of Exoplanets. Springer, Cham, Switzerland, p. 2771, doi:10.1007/978-3-319-55333-7_57
  * Mesquita et al. (2021) Mesquita A. L., Rodgers-Lee D., Vidotto A. A., 2021, MNRAS, 505, 1817–1826
  * Mesquita et al. (2022) Mesquita A. L., Rodgers-Lee D., Vidotto A. A., Atri D., Wood B. E., 2022, MNRAS, 509, 2091
  * Neronov et al. (2017) Neronov A., Malyshev D., Semikoz D. V., 2017, A&A, 606, A22
  * Parker (1958) Parker E. N., 1958, ApJ, 128, 664
  * Parker (1965) Parker E. N., 1965, Planet. Space Sci., 13, 9
  * Plavchan et al. (2020) Plavchan P., et al., 2020, Nature, 582, 497
  * Potgieter (2013) Potgieter M. S., 2013, Living Reviews in Solar Physics, 10, 3
  * Potgieter & Vos (2017) Potgieter M. S., Vos E. E., 2017, A&A, 601, A23
  * Potgieter et al. (2015) Potgieter M. S., Vos E. E., Munini R., Boezio M., Di Felice V., 2015, ApJ, 810, 141
  * Redfield & Linsky (2008) Redfield S., Linsky J. L., 2008, ApJ, 673, 283
  * Reiners & Basri (2008) Reiners A., Basri G., 2008, A&A, 489, L45
  * Rimmer & Helling (2013) Rimmer P. B., Helling C., 2013, ApJ, 774, 108
  * Rodgers-Lee et al. (2020) Rodgers-Lee D., Vidotto A. A., Taylor A. M., Rimmer P. B., Downes T. P., 2020, MNRAS, 499, 2124
  * Rodgers-Lee et al. (2021a) Rodgers-Lee D., Taylor A. M., Vidotto A. A., Downes T. P., 2021a, MNRAS, 504, 1519
  * Rodgers-Lee et al. (2021b) Rodgers-Lee D., Vidotto A. A., Mesquita A. L., 2021b, MNRAS, 508, 4696
  * Sadovski et al. (2018) Sadovski A. M., Struminsky A. B., Belov A., 2018, Astronomy Letters, 44, 324
  * Scheucher et al. (2020) Scheucher M., et al., 2020, ApJ, 893, 12
  * Schlickeiser (1989) Schlickeiser R., 1989, ApJ, 336, 243
  * See et al. (2015) See V., et al., 2015, MNRAS, 453, 4301
  * Selsis et al. (2007) Selsis F., Kasting J. F., Levrard B., Paillet J., Ribas I., Delfosse X., 2007, A&A, 476, 1373
  * Sokolov et al. (2013) Sokolov I. V., et al., 2013, ApJ, 764, 23
  * Stone et al. (2013) Stone E. C., Cummings A. C., McDonald F. B., Heikkila B. C., Lal N., Webber W. R., 2013, Science, 341, 150
  * Stone et al. (2019) Stone E. C., Cummings A. C., Heikkila B. C., Lal N., 2019, Nature Astronomy, 3, 1013
  * Strauss & Potgieter (2014) Strauss R. D., Potgieter M. S., 2014, Sol. Phys., 289, 3197
  * Strauss et al. (2012) Strauss R. D., Potgieter M. S., Büsching I., Kopp A., 2012, Ap&SS, 339, 223
  * Svensmark (2006) Svensmark H., 2006, Astronomische Nachrichten, 327, 871
  * Tinetti et al. (2021) Tinetti G., et al., 2021, arXiv e-prints, p. arXiv:2104.04824
  * van Holst et al. (2014) van Holst B., Sokolov I. V., Meng X., Jin M., Manchester W. B. I., Tóth G., Gombosi T. I., 2014, ApJ, 782, 81
  * Vidotto et al. (2014) Vidotto A. A., et al., 2014, MNRAS, 441, 2361
  * Vos & Potgieter (2015) Vos E. E., Potgieter M. S., 2015, ApJ, 815, 119
  * Wood (2004) Wood B. E., 2004, Living Reviews in Solar Physics, 1, 2
  * Wood et al. (2000) Wood B. E., Müller H.-R., Zank G. P., 2000, ApJ, 542, 493
  * Wood et al. (2001) Wood B. E., Linsky J. L., Müller H.-R., Zank G. P., 2001, ApJ, 547, L49
  * Wood et al. (2021) Wood B. E., et al., 2021, ApJ, 915, 37
  * Yadav et al. (2016) Yadav R. K., Christensen U. R., Wolk S. J., Poppenhaeger K., 2016, ApJ, 833, L28

##  Appendix A Galactic cosmic ray propagation in the solar system

In order to benchmark our code we calculate Galactic cosmic ray propagation
through the heliosphere. Here we explain the solar wind model and the results
we obtain at Earth.

To model the solar wind plasma, we follow the same methodologies as in
Kavanagh et al. (2021). Unlike AU Mic and Prox Cen, the mass-loss rate of the
solar wind has been measured in-situ to be $2\times 10^{-14}$ $M_{\sun}$ yr-1
(Cohen, 2011). The mass-loss rate of the wind is sensitive to the value used
for the Alfvén wave flux per unit magnetic field strength $S_{\textrm{A}}/B$
in the chromosphere (Boro Saikia et al., 2020; Kavanagh et al., 2021). For the
Sun, this value has been estimated directly from observations to be $1.1\times
10^{5}$ erg s-1 cm-2 G-1 (Sokolov et al., 2013).

Using this value, with the remaining Alfvén wave inputs being the same as
those in Kavanagh et al. (2021), we model the solar wind using a dipolar
magnetic field with a polar strength of 3 G, tilted 5$\degr$ with respect to
the rotation axis. Our 3D spherical grid extends out to 20 solar radii, and
contains around 4.5 million cells. The solar wind model used here is
representative of solar minimum with the heliospheric magnetic field pointing
outwards in the northern solar hemisphere (commonly referred to as A>0
cycles). The mass-loss rate we obtain from our model is $\sim 3\times
10^{-14}\,M_{\sun}~{}\textrm{yr}^{-1}$ and the solar wind terminal velocity is
$540\,$km s-1 (which is in the range of observed values, e.g. McComas et al.,
2000). By extrapolating the magnetic field strength of the 3D model from its
outer edge to 1 au, we found a magnetic field strength at Earth’s orbit of
$\sim 4\,$nT (which is in the range of observed values, e.g. Potgieter et al.,
2015). We set the size of the heliosphere to 122 au, as observed by Voyager 1
(Stone et al., 2013, 2019). Figure 6 shows the solar wind velocity (left) and
total magnetic field intensity (right) throughout the heliosphere. The total
magnetic field resembles an $\sim$ aligned dipole, with a bimodal wind
velocity distribution typical of solar wind at minimum (McComas et al., 2000).
Using this model, we calculate the propagation of Galactic cosmic rays as
discussed in Section 2.3.The flux of cosmic rays found at Earth is shown as
red dots in Figure 7 as a function of the kinetic energy.

Figure 6: Left: Equatorial plane cut of the solar wind velocity. Right:
Equatorial plane cut of the magnitude of the total magnetic field for the full
radial extent of the heliosphere. The values shown for small orbital distances
are from the 3D solar wind model and the values shown at large distances are
an extrapolation (see Appendix B for details). Figure 7: Differential
intensity of Galactic cosmic rays as a function of kinetic energy at Earth
representative of solar minimum conditions. The different linestyles represent
different models for comparison, namely: red dots for this study, green dashed
line for Strauss & Potgieter (2014), magenta dash-dotted line for Strauss &
Potgieter (2014) and blue dotted line for Potgieter & Vos (2017).

Figure 7 shows the differential intensity of Galactic cosmic rays at Earth as
a function of kinetic energy during the solar minimum activity and A>0 for
different works, namely: this study (red dots), Strauss & Potgieter (2014)
(green dashed line), Strauss & Potgieter (2014) (magenta dash-dotted line) and
Potgieter & Vos (2017) (blue dotted line). Overall, we get remarkably good
agreement between our model and the 3D models to which we compare our results.
The main discrepancies are most likely related to the different values of $B$,
$u_{\infty}$ and heliospheric sizes used in each model and the fact that we
use a 1D cosmic ray transport model while the other works use 3D models.

Another point to note is that Potgieter & Vos (2017) uses a ‘softening
parameter’ in the drift velocity equation so that the drift is reduced for
cosmic rays with momentum below 0.55 GeV/$c$. This reduction in the particle
drift is necessary to explain the observations at Earth. In our models, we did
not use a ‘softening parameter’ to reduce the particle drift velocities to
avoid including an extra free parameter in to our model.

##  Appendix B Extrapolation of stellar wind profiles to the edge of the
astrosphere

The outer boundary of the stellar wind models for AU Mic and Prox Cen from
Kavanagh et al. (2021) is at $100\,R_{\star}$, however the astrosphere extends
much further out. Hence, we extrapolate the quantities $u$, $B_{r}$,
$B_{\theta}$ and $B_{\phi}$ to account for the whole extent of the
astrosphere. Beyond $100\,R_{\star}$, the radial velocity is extrapolated as a
constant since the velocity has already reached its asymptotic value. To
generate a Parker spiral (Parker, 1958), $B_{r}$ continues to fall with
$r^{-2}$ and $B_{\phi}$ with $r^{-1}$. In the case of $B_{\theta}$, we use a
power law fit for each system. For Prox Cen, $B_{\theta}$ falls with
$r^{-1.7}$, for ‘AU Mic low’ with $r^{-2.2}$ and for ‘AU Mic high’ with
$r^{-4.7}$. For the Sun we extrapolate $B_{\theta}$ as $1/r^{4.6}$ beyond
$r>0.04\,$au. Figures 1 and 6 show the extrapolated values of the wind
velocity and total magnetic field for the stars studied here. Note, our
stellar wind model does not take into account the termination shock (region
where the stellar wind properties change due to the interactions with the
local interstellar medium).