# Generalization of Auto-Regressive Hidden Markov Models to Non-Linear
Dynamics and Unit Quaternion Observation Space

Michele Ginesi1 and Paolo Fiorini2 *This project has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme, ARS (Autonomous Robotic Surgery) project,
grant agreement No. 742671. 1 Department of Computer Science, University of
Verona, Strada le Grazie 15, Verona, Italy<EMAIL_ADDRESS>Department
of Computer Science, University of Verona, Strada le Grazie 15, Verona, Italy
<EMAIL_ADDRESS>
###### Abstract

Latent variable models are widely used to perform unsupervised segmentation of
time series in different context such as robotics, speech recognition, and
economics. One of the most widely used latent variable model is the Auto-
Regressive Hidden Markov Model (ARHMM), which combines a latent mode governed
by a Markov chain dynamics with a linear Auto-Regressive dynamics of the
observed state.

In this work, we propose two generalizations of the ARHMM. First, we propose a
more general AR dynamics in Cartesian space, described as a linear combination
of non-linear basis functions. Second, we propose a linear dynamics in unit
quaternion space, in order to properly describe orientations. These extensions
allow to describe more complex dynamics of the observed state.

Although this extension is proposed for the ARHMM, it can be easily extended
to other latent variable models with AR dynamics in the observed space, such
as Auto-Regressive Hidden semi-Markov Models.

##  I INTRODUCTION

_Hidden Markov Models_ (HMMs) [1, 2, 3] are a type of graphical model widely
used in speech recognition [4, 5], hand-writing recognition [6], natural
language modeling [7], and to segment kinematics in the context of _Minimally
Invasive Surgery_ [8, 9, 10]. The model consists of a Markov chain governing
the evolution of an _hidden_ (or _latent_) _mode_. At each time $t$, the
hidden mode $z_{t}$ emits an observation $\mathbf{y}_{t}$ with a certain
probability distribution $p(\mathbf{y}_{t}|z_{t})$.

A well-known generalization of HMMs is the _Auto-Regressive_ HMM (ARHMM) [11,
4, 12, 13, 14, 15, 16], in which the observation $\mathbf{y}_{t}$ at time $t$
is given by the hidden mode $z_{t}$ and the previous observed state
$\mathbf{y}_{t-1}$ via a linear Auto-Regressive dynamic. This means that the
current latent mode $z_{t}$ does not emit an observation, but instead it
describes the (linear) vector field governing the evolution of the observed
state.

This type of model, and its generalizations, have been successfully used to
temporally segment robot kinematics. For instance, in [15] a _state-based
transition_ ARHMM (STARHMM) was developed. This model adds a conditional
dependence between the current observed state and the next hidden mode.
Similarly, in [17] this model was improved to learn skills by creating a
primitive library.

In the literature, most of the improvements to ARHMM have been implemented in
the definition of the latent or observation space (e.g. by adding hidden
variables [18, 17]), by adding conditional probabilities relation between
modes (e.g. by making the next hidden mode dependent on the current observed
state [15, 17]), or by increasing the number of previous states taken into
account in the definition of the AR dynamics. That is, most of the
generalizations of the ARHMM have focused on the _topology_ of the graphical
model.

On the other hand, at the best of the authors’ knowledge, no generalization of
the definition of dynamics of the observed state exists in the literature. For
this reason, in this work, we propose to generalize the ARHMM to allow Non-
Linear AR dynamics in Cartesian space and linear dynamics in Unit Quaternion
space. The main advantage of our proposed improvement lies in the fact that it
can be straightforwardly adapted to work in combination with the modification
of the topology of the latent variable model mentioned above.

The paper is structured as follows. In Section II we present the theory behind
ARHMM, recalling the _Expectation Maximization_ (EM) algorithm used to infer
the set of model parameters from data. In Section III we present our proposed
modification to the ARHMM model and the modification to the EM algorithm to
infer the model’s parameter. In Section IV we propose the experiments
validating our proposed approach. In Section V we present the conclusion of
the work and possible future extensions.

##  II AUTO REGRESSIVE HIDDEN MARKOV MODELS

In this section, we provide a formal definition of the Auto-Regressive Hidden
Markov Model.

In details, an ARHMM is a model defined as follows:

$\mathcal{H}\doteq\left\\{\mathcal{S},\mathcal{Y},\Theta=\left\\{\bm{\varpi},\mathbf{T},\left\\{\mathbf{A}_{s},\mathbf{b}_{s},\bm{\Sigma}_{s}\right\\}_{s\in\mathcal{S}}\right\\}\right\\}.$

$\mathcal{S}=\left\\{1,2,\ldots,S\right\\}$ is the set of hidden modes, and
the mode at time $t\in{1,2,\ldots,T}$ is denoted by $z_{t}$.
$\mathcal{Y}=\mathbb{R}^{d}$ is the observation space, and the observed state
at time $t$ is denoted by $\mathbf{y}_{t}\in\mathcal{Y}$. $\Theta$ is the set
of model parameters. Vector $\bm{\varpi}=[\varpi_{i}]_{i\in\mathcal{S}}$
defines the initial mode probability such that $\varpi_{i}=\Pr(z_{1}=i),$ for
$i\in\mathcal{S}$. Matrix
$\mathbf{T}=[\texttt{t}_{i\,j}]_{i\in\mathcal{S}}^{j\in\mathcal{S}}$ is the
transition probabilities matrix between hidden modes such that
$\texttt{t}_{i\,j}=\Pr(z_{t+1}=j|z_{t}=i).$ The emissions are modeled by the
following Auto-Regressive (AR) dynamics with Gaussian white noise:

$\mathbf{y}_{t}|z_{t}=s,\,\mathbf{y}_{t-1}\sim\mathcal{N}\left(\mathbf{y}_{t}|\mathbf{A}_{s}\mathbf{y}_{t-1}+\mathbf{b}_{s},\bm{\Sigma}_{s}\right).$
(1)

The graphical representation of the model is given in Figure 1. In this
representation, a grey dot represents an observed variable, while a white dot
denotes a hidden variable. An arrow between two nodes $a$ and $b$ means that
the probability of $b$ depends on the value of $a$. For a refresh on graphical
model representation, we suggest [19].

Two main algorithms characterize the ARHMM (and the latent-variable models in
general): the _Expectation Maximization_ (EM) algorithm and the _Viterbi_
algorithm. The former is used to infer the set of parameters $\Theta$ by
maximizing the likelihood of the observed data $p(\mathbf{Y}|\Theta).$ The
latter, instead, allows extracting the sequence of latent mode
$\hat{\mathbf{z}}$ that maximize the joint probability distribution
$p(\mathbf{Y},\mathbf{z}|\Theta)$ for a given observed sequence $\mathbf{Y}.$
In short, the EM algorithm allows to learn the model parameters from the data,
while the Viterbi algorithm allows to segment an observed sequence once the
model parameters are known.

$z_{1}$$z_{t-1}$$z_{t}$$z_{t+1}$$z_{T}$$\mathbf{y}_{0}$$\mathbf{y}_{1}$$\mathbf{y}_{t-1}$$\mathbf{y}_{t}$$\mathbf{y}_{t+1}$$\mathbf{y}_{T}$$\mathbf{y}_{0}$$\mathbf{y}_{1}$$\mathbf{y}_{t-1}$$\mathbf{y}_{t}$$\mathbf{y}_{t+1}$$\mathbf{y}_{T}$

Figure 1: Graphical representation of ARHMM.

###  II-A EM Algorithm in brief

The EM algorithm consists of an iterative procedure that, at each iteration,
maximizes, w.r.t. $\Theta$, the auxiliary function

$\medmath{\begin{multlined}Q(\Theta,\Theta^{\text{old}})\doteq
Q_{\text{init}}(\Theta,\Theta^{\text{old}})+Q_{\text{trans}}(\Theta,\Theta^{\text{old}})+Q_{\text{out}}(\Theta,\Theta^{\text{old}})\end{multlined}Q(\Theta,\Theta^{\text{old}})\doteq
Q_{\text{init}}(\Theta,\Theta^{\text{old}})+Q_{\text{trans}}(\Theta,\Theta^{\text{old}})+Q_{\text{out}}(\Theta,\Theta^{\text{old}})}$
(2)

where $\Theta^{\text{old}}$ is the current guess of the set of parameters, and
the three terms of the sum are defined as

$\displaystyle\medmath{\begin{multlined}Q_{\text{init}}(\Theta,\Theta^{\text{old}})\doteq\sum_{z_{1}\in\mathcal{S}}\log
p(z_{1}|\Theta)p(z_{1}|\mathbf{Y},\Theta^{\text{old}}),\end{multlined}Q_{\text{init}}(\Theta,\Theta^{\text{old}})\doteq\sum_{z_{1}\in\mathcal{S}}\log
p(z_{1}|\Theta)p(z_{1}|\mathbf{Y},\Theta^{\text{old}}),}$ (3b)
$\displaystyle\medmath{\begin{multlined}Q_{\text{trans}}(\Theta,\Theta^{\text{old}})\doteq\sum_{t=1}^{T-1}\sum_{z_{t},z_{t+1}\in\mathcal{S}}\log
p(z_{t+1}|z_{t},\Theta)p(z_{t},z_{t+1}|\mathbf{Y},\Theta^{\text{old}}),\end{multlined}Q_{\text{trans}}(\Theta,\Theta^{\text{old}})\doteq\sum_{t=1}^{T-1}\sum_{z_{t},z_{t+1}\in\mathcal{S}}\log
p(z_{t+1}|z_{t},\Theta)p(z_{t},z_{t+1}|\mathbf{Y},\Theta^{\text{old}}),}$ (3d)
$\displaystyle\medmath{\begin{multlined}Q_{\text{out}}(\Theta,\Theta^{\text{old}})\doteq\sum_{t=0}^{T-1}\sum_{z_{t+1}\in\mathcal{S}}\log
p(\mathbf{y}_{t+1}|z_{t+1},\mathbf{y}_{t},\Theta)p(z_{t+1}|\mathbf{Y},\Theta^{\text{old}}).\end{multlined}Q_{\text{out}}(\Theta,\Theta^{\text{old}})\doteq\sum_{t=0}^{T-1}\sum_{z_{t+1}\in\mathcal{S}}\log
p(\mathbf{y}_{t+1}|z_{t+1},\mathbf{y}_{t},\Theta)p(z_{t+1}|\mathbf{Y},\Theta^{\text{old}}).}$
(3f)

During the _Expectation Step_ , the quantities that depend on
$\theta^{\text{old}}$, that is
$p(z_{t},z_{t+1}|\mathbf{Y},\Theta^{\text{old}})$ and
$p(z_{t+1}|\mathbf{Y},\Theta^{\text{old}})$, are computed using the _forward-
backward algorithm_ [19]. During the _Maximization Step_ , the set of
parameters $\Theta$ that maximizes $Q_{\text{init}}$, $Q_{\text{trans}}$ and
$Q_{\text{out}}$ is computed. These two steps are repeated until convergence.
Convergence is reached when the changes in the log-likelihood $\log
p(\mathbf{Y},\mathbf{z})$ is below a given tolerance.

##  III GENERALIZATION OF ARHMMS TO DIFFERENT DYNAMICS

In this section, we present our proposed generalizations of the AR dynamic. In
particular, in Section III-A we discuss our proposed non-linear dynamics in
Cartesian space; while in Section III-B we present our proposed linear
dynamics in Unit Quaternion space.

###  III-A Cartesian Non-Linear ARHMM

When generalizing the AR dynamics in Cartesian space, we aim at modifying the
linear AR dynamics shown in (1) with a more general dynamic
$\mathbf{f}_{s}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ for each hidden mode
$s\in\mathcal{S}.$ To do so, we will use a _basis functions_ -based
formulation so that the non-linear function $\mathbf{f}_{s}$ is written as a
linear combination of non-linear basis functions
$\left\\{\varphi_{n}:\mathbb{R}^{d}\to\mathbb{R}\right\\}_{n=0,1,\ldots,N}:$

$\mathbf{f}_{s}(\mathbf{y}_{t})|_{i}=f_{i}^{(s)}(\mathbf{y}_{t})=\sum_{j=0}^{N}\omega_{i\,j}^{(s)}\varphi_{j}(\mathbf{y}_{t}),$
(4)

where $\omega_{i\,j}^{(s)}$ is the set of _weights_ that describes the vector
field. By defining the _weight matrix_
$\bm{\Omega}_{s}\doteq[\omega_{i\,j}^{(s)}]_{i=1,2,\ldots,d}^{j=0,1,\ldots,N}$
and the _non-linear map_
$\bm{\varphi}(\mathbf{y})\doteq[\varphi_{j}(\mathbf{y})]_{j=0,1,\ldots,N},$ we
have that formulation (4) can be written in matrix-vector notation as

$\mathbf{f}_{s}(\mathbf{y}_{t})=\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t}).$
(5)

Thus, (1) now reads

$\mathbf{y}_{t}|z_{t}=s,\,\mathbf{y}_{t-1}\sim\mathcal{N}\left(\mathbf{y}_{t}|\bm{\Omega}_{s}\bm{\varphi}\left(\mathbf{y}_{t-1}\right),\bm{\Sigma}_{s}\right).$
(6)

We remark that, since only the formulation of the AR dynamics changes from the
classical ARHMM to the non-linear case, the graphical model representation is
the same as in Figure 1.

During the _Expectation_ step, the quantities that depend on the current guess
of the set $\Theta^{\text{old}}$ of parameters have to be computed. Since the
conditional independence properties of the model depend exclusively on its
graphical representation [19], we have that the E-step for the _Non-Linear
ARHMM_ (NL-ARHMM) is identical to the linear case. The only difference is that
the emission probability $p(\mathbf{y}_{t}|z_{t},\mathbf{y}_{t-1})$ is given
by formula (6) instead of (1).

Similarly, in the _Maximization_ step, maximization of (3b) and (3d) is
achieved in the same way as the linear ARHMM case, since there is no
dependence on the emission probability. Thus, the only difference in the EM
algorithm between a linear ARHMM and the NL-ARHMM lies in the maximization of
quantity $Q_{\text{out}}.$

###### Proposition III.1

Denote by $\gamma(t)$ the quantity $p(z_{t}|\mathbf{Y},\Theta)$ so that
$\gamma_{s}(t)\doteq\Pr(z_{t}=s|\mathbf{Y},\Theta^{\text{old}}).$ Then,
maximization of $Q_{\text{out}}(\Theta,\Theta^{\text{old}})$ in (3f) is
achieved by setting

$\medmath{\bm{\Omega}_{s}=\\!\\!\begin{multlined}\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{y}_{t+1}\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\right)\\!\\!\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\bm{\varphi}(\mathbf{y}_{t})\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\right)^{-1}\end{multlined}\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{y}_{t+1}\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\right)\\!\\!\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\bm{\varphi}(\mathbf{y}_{t})\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\right)^{-1}}$
(7a) and
$\bm{\Sigma}_{s}=\frac{\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{e}^{(s)}_{t}\mathbf{e}^{(s)}_{t}{}^{\intercal}}{\sum_{t=0}^{T}\gamma_{s}(t+1)};$
(7b)

where the vector $\mathbf{e}_{t}^{(s)}\in\mathbb{R}^{d}$ is defined as

$\mathbf{e}_{t}^{(s)}=\mathbf{y}_{t+1}-\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t}).$

###### Proof:

We aim at maximizing (by dropping $\Theta$ for notation simplicity) the
quantity

$\widetilde{Q}(\Theta,\Theta^{\text{old}})=\sum_{t=0}^{T-1}\sum_{s=1}^{S}\gamma_{s}(t+1)\log
p(\mathbf{y}_{t+1}|z_{t+1}=s,\mathbf{y}_{t}).$

To maximize this quantity, we can maximize over each hidden mode
$s\in\mathcal{S}$ independently:

$\begin{multlined}{\bm{\Omega}_{s}^{\star},\bm{\Sigma}_{s}^{\star}}=\arg\max_{\bm{\Omega}_{s},\bm{\Sigma}_{s}}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log
p(\mathbf{y}_{t+1}|z_{t+1}=s,\mathbf{y}_{t}).\end{multlined}{\bm{\Omega}_{s}^{\star},\bm{\Sigma}_{s}^{\star}}=\arg\max_{\bm{\Omega}_{s},\bm{\Sigma}_{s}}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log
p(\mathbf{y}_{t+1}|z_{t+1}=s,\mathbf{y}_{t}).$

Since the emission probability is given by (6), the function to maximize reads

$\medmath{\begin{multlined}\widehat{Q}_{\text{out}}=-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log|\bm{\Sigma}_{s}|\\\
-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\varphi(\mathbf{y}_{t})\right)^{\intercal}\bm{\Sigma}_{s}^{-1}\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\varphi(\mathbf{y}_{t})\right).\end{multlined}\widehat{Q}_{\text{out}}=-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log|\bm{\Sigma}_{s}|\\\
-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\varphi(\mathbf{y}_{t})\right)^{\intercal}\bm{\Sigma}_{s}^{-1}\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\varphi(\mathbf{y}_{t})\right).}$
(8)

We start by discussing the maximization over the weight matrix
$\bm{\Omega}_{s}$. The gradient of $\widehat{Q}_{\text{out}}$ w.r.t.
$\bm{\Omega}_{s}$ is

$\begin{multlined}\nabla_{\bm{\Omega}_{s}}\widehat{Q}=\bm{\Sigma}_{s}^{-1}\Bigg{(}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{y}_{t+1}\bm{\varphi}(\mathbf{y})^{\intercal}\hskip
56.9055pt\\\
-\bm{\Omega}_{s}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\bm{\varphi}(\mathbf{y}_{t})\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\Bigg{)}.\end{multlined}\nabla_{\bm{\Omega}_{s}}\widehat{Q}=\bm{\Sigma}_{s}^{-1}\Bigg{(}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{y}_{t+1}\bm{\varphi}(\mathbf{y})^{\intercal}\hskip
56.9055pt\\\
-\bm{\Omega}_{s}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\bm{\varphi}(\mathbf{y}_{t})\bm{\varphi}(\mathbf{y}_{t})^{\intercal}\Bigg{)}.$

By setting it to zero, we get formula (7a).

We now discuss the maximization over the weight matrix $\bm{\Sigma}_{s}$. The
gradient of $\widehat{Q}_{\text{out}}$ w.r.t. $\bm{\Sigma}_{s}$ is

$\nabla_{\bm{\Sigma}_{s}^{-1}}\widehat{Q}=\begin{multlined}\frac{1}{2}\bm{\Sigma}_{s}\sum_{t=0}^{T-1}\gamma_{s}(t+1)-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\\\
\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t})\right)\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t})\right)^{\intercal},\end{multlined}\frac{1}{2}\bm{\Sigma}_{s}\sum_{t=0}^{T-1}\gamma_{s}(t+1)-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\\\
\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t})\right)\left(\mathbf{y}_{t+1}-\bm{\Omega}_{s}\bm{\varphi}(\mathbf{y}_{t})\right)^{\intercal},$

which, by setting it to zero, gives formula (7b).

∎

We remark that the matrix $\bm{\Omega}_{s}$ in (7b) is the “old” guess, and
not that given by the update formula (7a). Thus, when implementing the EM-
algorithm, the update of the covariance matrices $\bm{\Sigma}_{s}$ should be
performed before the update of the weight matrices $\bm{\Omega}_{s}$.

####  III-A1 Example of Basis Functions

In this section, we present some examples of basis functions.

The first set shows how to interpret the classical, linear ARHMM as a
particular case of NL-ARHMM. Indeed, by defining
$\bm{\varphi}^{\text{(lin)}}(\mathbf{y})\doteq\left[1,y_{1},y_{2},\ldots,y_{d}\right]^{\intercal},$
we have that the resulting weight matrix is a block matrix with the offset
$\mathbf{b}$ and the linear map $\mathbf{A}$ in (1) as
$\bm{\Omega}=\left[\mathbf{b},\mathbf{A}\right].$

The second well known set of basis functions is the family of _Gaussian Radial
Basis Functions_ (GRBFs). Given a set $\\{\bm{\mu}_{i}\\}_{i=1,2,\ldots,N}$ of
centers and a set $\\{\bm{\varSigma}_{i}\\}_{i=1,2,\ldots,N}$ of covariance
matrices, we define the basis functions as

$\displaystyle\varphi^{\text{(grbf)}}_{0}(\mathbf{y})$ $\displaystyle\doteq
1,$ $\displaystyle\varphi^{\text{(grbf)}}_{i}(\mathbf{y})$
$\displaystyle\doteq\begin{multlined}\exp\left(-(\mathbf{y}-\bm{\mu}_{i})\bm{\varSigma}_{i}^{-1}(\mathbf{y}-\bm{\mu}_{i})\right),i=1,\ldots,N.\end{multlined}\exp\left(-(\mathbf{y}-\bm{\mu}_{i})\bm{\varSigma}_{i}^{-1}(\mathbf{y}-\bm{\mu}_{i})\right),i=1,\ldots,N.$

Usually, the covariance matrices are set to be a multiple of the identity
matrix $\bm{\varSigma}_{i}=\varsigma_{i}\text{{{I}d}}_{d}.$ The main drawback
of this family of basis functions lies in the fact that, even for small values
of the dimension $d$ of the observed space, a high number of basis functions
is needed to ‘cover’ the space. Indeed, this family of basis functions is
usually adopted when learning functions in bounded domains (e.g. in _Dynamic
Movement Primitives_ , where the time domain is fixed [20]). For this reason,
we suggest to use this set of basis functions when dealing with bounded
observation spaces, such as $d-$dimensional cubes $[0,1]^{d}$.

Finally, a third family of basis functions is the set of polynomial functions
up to a degree $k$:

$\bm{\varphi}^{(\mathcal{P}_{k})}(\mathbf{y})\doteq\text{col}\left(\left\\{\prod_{i=1}^{d}y_{i}^{c_{i}}:c_{i}\in\mathbb{N},\,\sum_{i=1}^{d}c_{i}\leq
k\right\\}\right).$ (9)

For instance, assume that $d=2$ ($\mathbf{y}\in\mathcal{Y}=\mathbb{R}^{2}$)
and $k=3$. Then the basis functions are

$\bm{\varphi}^{(\mathcal{P}_{3})}(\mathbf{y})=\left[1,y_{1},y_{2},y_{1}^{2},y_{1}y_{2},y_{2}^{2},y_{1}^{3},y_{1}^{2}y_{2},y_{1}y_{2}^{2},y_{2}^{3}\right]^{\intercal}.$

As we will show in Section IV-B, this type of basis functions is particularly
useful in low-dimension spaces, where linear dynamics are not able to describe
complex evolutions of the observed state. On the other hand, they are less
useful in higher dimensional spaces for two main reasons: firstly, linear
dynamics are able to describe more complex evolutions; secondly, non-linearity
in higher dimensional spaces require a huge amount of parameters (i.e. the
number of columns of $\bm{\Omega}$ increases) risking over-fitting and
numerical inefficiency.

###  III-B Unit-Quaternion Linear ARHMM

Orientations can be modeled in different ways. There are two preferred spaces
to describe orientations: the space $SO(3)$ of _orthogonal $3\times 3$
matrices_, and the space $\mathbb{S}^{3}$ of unit quaternions. Both
representations are singularity-free. However, unit quaternions are preferred
since they require only four variable to be described, instead of the nine
parameters of $3\times 3$ matrices.

The main difficulty when dealing with dynamics in $\SS^{3}$ lies in the fact
that the resulting quaternion must still be of unitary norm. To solve this
problem, we take advantage of two properties of quaternions. The first is the
fact that the exponential of a quaternion with null real part always results
in a quaternion of unit norm:
$\left\|\text{Exp}(0+a\texttt{i}+b\texttt{j}+c\texttt{k})\right\|=1,\quad\forall
a,b,c,\in\mathbb{R}.$ The second is the fact that the norm of the product
between two quaternion $p$ and $q$ is the product of the norms:
$\left\|p*q\right\|=\left\|p\right\|\,\left\|q\right\|.$ Thus, we propose to
describe a dynamics in $\SS^{3}$ by 3 parameters $a,b,c\in\mathbb{R}$ such
that

$q_{t+1}=\text{Exp}(a\texttt{i}+b\texttt{j}+c\texttt{k})*q_{t}.$

From here, we will denote as ‘vec’ the function that maps a quaternion to its
vector formulation:

$\medmath{\text{vec}:\SS^{3}\to\mathbb{R}^{4},\qquad
q=q_{r}+q_{i}\texttt{i}+q_{j}\texttt{j}+q_{k}\texttt{k}\mapsto\mathbf{q}=[q_{r},q_{i},q_{j},q_{k}]^{\intercal}.}$

With this notation, we have that in this model the probability of the next
observation reads

$q_{t}|z_{t}=s,q_{t-1}\sim\\\
\mathcal{N}\big{(}\text{vec}(q_{t})|\text{vec}\left(\text{Exp}(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k})*q_{t-1}\right),\bm{\Sigma}_{s}\big{)}$
(10)

As in the case of Cartesian Non-Linear ARHMM, we need only to explain how to
maximize the _emission probability_. Therefore, in the _Maximization step_ we
have to maximize the quantity $Q_{\text{out}}$ in (3f) where the next-state
probability is given by (10). Differently from what we showed in Proposition
III.1, there is not a closed form solution to the maximization problem.
Indeed, by defining
$\gamma_{s}(t)\doteq\Pr(z_{t}=s|\mathbf{Y},\Theta^{\text{old}})$ we have that
the function to maximize reads, similarly to (8):

$\medmath{\begin{multlined}\widehat{Q}_{\text{out}}=-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log\left|\bm{\Sigma}_{s}\right|\\\
-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\text{vec}(q_{t+1}-\mu_{t+1}^{(s)})^{\intercal}\bm{\Sigma}^{-1}\text{vec}(q_{t+1}-\mu_{t+1}^{(s)})\end{multlined}\widehat{Q}_{\text{out}}=-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\log\left|\bm{\Sigma}_{s}\right|\\\
-\frac{1}{2}\sum_{t=0}^{T-1}\gamma_{s}(t+1)\text{vec}(q_{t+1}-\mu_{t+1}^{(s)})^{\intercal}\bm{\Sigma}^{-1}\text{vec}(q_{t+1}-\mu_{t+1}^{(s)})}$
(11)

where $\mu_{t+1}^{(s)}$ is the _expected orientation_ at time $t+1$ for hidden
mode $s$:

$\mu_{t+1}^{(s)}\doteq\text{Exp}(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k})*q_{t}.$

Similarly to Proposition III.1, maximization of (11) is achieved by setting

$\bm{\Sigma}_{s}=\frac{\sum_{t=0}^{T-1}\gamma_{s}(t+1)\mathbf{e}_{t+1}^{(s)}\cdot\mathbf{e}_{t+1}^{(s)}{}^{\intercal}}{\sum_{t=0}^{T-1}\gamma_{s}(t+1)}$

where $\mathbf{e}_{t+1}^{(s)}$ is the difference between the actual
orientation and the predicted one at time $t+1$ assuming that the hidden mode
$s$ is active:

$\mathbf{e}_{t+1}^{(s)}\doteq\text{vec}(q_{t+1}-\mu_{t+1}^{(s)}).$

On the other hand, maximization over the dynamics parameters can be formulated
as

$\medmath{\begin{multlined}(a_{s}^{\star},b_{s}^{\star},c_{s}^{\star})=\arg\max_{(a_{s},b_{s},c_{s})}\\\
\left(-\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left\|\text{vec}\left(q_{t+1}-\text{Exp}\left(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k}\right)*q_{t}\right)\right\|_{\bm{\Sigma}^{-1}}^{2}\right)\end{multlined}(a_{s}^{\star},b_{s}^{\star},c_{s}^{\star})=\arg\max_{(a_{s},b_{s},c_{s})}\\\
\left(-\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left\|\text{vec}\left(q_{t+1}-\text{Exp}\left(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k}\right)*q_{t}\right)\right\|_{\bm{\Sigma}^{-1}}^{2}\right)}$

where $\left\|\cdot\right\|_{\mathbf{A}}$ is the semi-norm induced by a
Symmetric Semi-Positive Definite (SPD) matrix $\mathbf{A}$:
$\left\|\mathbf{x}\right\|_{\mathbf{A}}\doteq\sqrt{\mathbf{x}^{\intercal}\mathbf{A}\mathbf{x}}.$
Unlikely the Cartesian formulation, this problem cannot be solved in closed
form. Therefore, to solve the Maximization step, we rewrite the maximization
problem as an equivalent minimization problem

$\medmath{\begin{multlined}(a_{s}^{\star},b_{s}^{\star},c_{s}^{\star})=\arg\min_{(a_{s},b_{s},c_{s})}\\\
\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left\|\text{vec}\left(q_{t+1}-\text{Exp}\left(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k}\right)*q_{t}\right)\right\|_{\bm{\Sigma}^{-1}}^{2}\right)\end{multlined}(a_{s}^{\star},b_{s}^{\star},c_{s}^{\star})=\arg\min_{(a_{s},b_{s},c_{s})}\\\
\left(\sum_{t=0}^{T-1}\gamma_{s}(t+1)\left\|\text{vec}\left(q_{t+1}-\text{Exp}\left(a_{s}\texttt{i}+b_{s}\texttt{j}+c_{s}\texttt{k}\right)*q_{t}\right)\right\|_{\bm{\Sigma}^{-1}}^{2}\right)}$

and solve it via any minimization algorithm (e.g. _gradient descent_).

The implementation (in Python 3.10) of our proposed generalizations is
available at https://github.com/mginesi/nl_arhmm.

###### Remark III.1

The use of our proposed generalizations can be combined with other extensions
to the ARHMM model. For instance, in the EM algorithm of _Auto-Regressive
Hidden semi Markov Models_ , the maximization formula for the emission
probabilities can be easily extended in a similar fashion to what we proposed
here.

###### Remark III.2

We do not discuss the generalization of the Viterbi algorithm [21] since its
most general formulation works with any latent-mode model with Markov dynamic,
independently of the dynamics of the observed state.

##  IV EXPERIMENTS

In this section, we will compare the results obtained with the NL-ARHMM
against the linear ARHMM, showing that our proposed approach results in higher
segmentation scores.

###  IV-A Validation Test

At first, we show that the EM algorithm for the NL-ARHMM allows learning the
parameters of the model. To do so, we define a NL-ARHMM with fixed parameters
and use it to generate data samples. In our test, we set both the number of
hidden modes and the dimension of the continuous state to 2: $S=2,d=2.$ The
initial mode distribution is set to $\bm{\varpi}=[0.5,0.5]^{\intercal}$ and
the transition probability matrix is

$\mathbf{T}=\begin{bmatrix}0.95&0.05\\\ 0.05&0.95\end{bmatrix}.$

The non-linear dynamics is written as follows. We start by defining the vector
field

$\mathbf{f}\left(\begin{bmatrix}y_{1}\\\
y_{2}\end{bmatrix}\right)=\begin{bmatrix}y_{1}{}^{3}+y_{2}{}^{2}y_{1}-y_{1}-y_{2}\\\
y_{2}{}^{3}+y_{1}{}^{2}y_{2}+y_{1}-y_{2}\end{bmatrix}.$

Then the dynamics of the two hidden modes are obtained using the Euler method
with time step $\delta_{t}=0.05$ for the dynamics $\mathbf{f}$ and its
opposite $-\mathbf{f}$, with Gaussian noise with standard deviation
$\varsigma=5\texttt{e}-03$:

$\displaystyle\mathbf{y}_{t+1}|z_{t+1}=1,\mathbf{y}_{t}\sim\mathcal{N}(\mathbf{y}_{t+1}|\mathbf{y}_{t}+\delta_{t}\mathbf{f}(\mathbf{y}_{t}),\varsigma\text{{{I}d}}_{2}),$
$\displaystyle\mathbf{y}_{t+1}|z_{t+1}=2,\mathbf{y}_{t}\sim\mathcal{N}(\mathbf{y}_{t+1}|\mathbf{y}_{t}-\delta_{t}\mathbf{f}(\mathbf{y}_{t}),\varsigma\text{{{I}d}}_{2}).$

Figure 2:  Results for the validation test presented in Section IV-A. The top
figure shows the components of the trajectory, the second plot shows the
“true” segmentation (the first mode is represented in purple, while the second
mode is represented in yellow), the third plot shows the segmentation obtained
with the learned ARHMM, and the fourth plot shows the segmentation obtained
with the learned NL-ARHMM.

(a) $d=1$,

(b) $d=2$

(c) $d=3$

Figure 3:  Comparison between linear and polynomial ARHMM. In all three
figures the first plot shows the components $y_{i}$ of the observed dynamics,
the second plot shows the true segmentation, the third, fourth, and fifth
plots shows the segmentations obtained by, respectively, the linear,
quadratic, and cubic formulations of the ARHMM (i.e. $k=1,2,3$ in (9)).

We use this NL-ARHMM to generate fifty samples with $T=100$ steps each. Next,
we initialize both a NL-ARHMM (with polynomial basis functions of degree
$k=2$) and a linear ARHMM, and apply the EM algorithm to learn the model
parameters. Finally, we use the Viterbi algorithm to segment a new trajectory
generated from the “true” NL-ARHMM. Before applying the EM algorithm we
standardize the dataset. That is, we translate it and multiply by a constant
to have null mean and unit variance.

Figure 2 shows one result for these tests. As can be observed, the NL-ARHMM is
able to properly segment the trajectory, being able to capture more complex
dynamics. On the other hand, the linear ARHMM results in a poorer
segmentation.

###  IV-B Comparison of Different Degrees Polynomial Functions

In this Section, we propose some tests to show that Non-Linear basis functions
are more effective when the observation space $\mathcal{Y}$ is low
dimensional. On the other hand, for high-dimensional observation space the
difference between the two is less noticeable.

To show this aspect we simulate three different ARHMM with non linear
dynamics. Each of them has two hidden modes but differs in the dimension $d$
of the observation space; $d=1,2,3$. To present a fair comparison, all the six
dynamics (two for each value of $d$) are non-polynomial: in this way none of
the model we will compare (linear, quadratic, and cubic) is able to perfectly
describe the dynamics.

Simularly to the previous test, for each value of $d$ we generate a total of
fifty trajectories of 100 steps. Figure 3 shows the results of these tests.

As it can be observed, with $d=1$ (Figure 3a) the linear model is not able to
properly describe the non-lienar dynamics, failing in the segmentation of the
obtained trajectory. On the other hand, for $d=2$ and $3$ (Figures 3b and 3c
respectively), the differences between the three dynamical model (linear,
quadratic, and cubic) significantly reduce.

These tests show that for lower-dimensional observation spaces, the
improvements in using the Non-Linear formulation for the ARHMM dynamics are
more significative. We argue that the motivation lies in the fact that for
high-dimensional data, linear vector fields can describe a wider range of
dynamics, since the number of parameters in the dynamics grows quadratically
with the space dimension (to be precise, the number of parameters for each
dynamics is $d^{2}+d$: $d^{2}$ elements for the matrix $\mathbf{A}_{s}$ and
$d$ for the vector $\mathbf{b}_{s}$ in (1)). Oh the other hand, for smaller
observation spaces, linear dynamics are too limited and the usage of
polynomial functions in the definition of the dynamics provides a broader set
of possible behaviors.

For this reason, in Section IV-C we will propose a model in which the
Cartesian position of the robot end-effector is modeled by a linear dynamics,
while the gripper angle (a 1-dimensional variable) will be modeled using a
polynomial dynamics.

###  IV-C Experiments on Real Setups

To prove the effectiveness of our generalizations, we propose a model
combining both Non-Linear Cartesian and Linear Unit-Quaternion dynamics. In
particular, we propose a pose+gripper model in which, for each hand of the
robot, we model the position with a linear Cartesian model, the end-effector
orientation with a linear Unit-Quaternion model, and the gripper angle with a
quadratic model. The decision to use a linear model for the position and a
quadratic model for the gripper angle are motivated by the tests presented in
Section IV-B: for 3-dimensional spaces, linear, quadratic, and cubic models
are almost indistinguishable (Figure 3c) and we thus choose the simplest one.
On the other hand, for 1-dimensional trajectories, linear models fail to
describe the complexity of the behavior, while a quadratic and a cubic one
give similar results (Figure 3c).

Additionally, we assume that position $\mathbf{x}^{h}_{t}$, orientation
$q^{h}_{t}$, and gripper angle $\vartheta^{h}_{t}$ of each arm $h$ are
independent of each others when the hidden mode is given:

$p\left(\begin{subarray}{c}\mathbf{x}^{1}_{t+1},q^{1}_{t+1},\vartheta^{1}_{t+1},\\\
\\\ \ldots,\\\ \\\
\mathbf{x}^{H}_{t+1},q^{H}_{t+1},\vartheta^{H}_{t+1}\end{subarray}\Bigg{|}\begin{subarray}{c}\mathbf{x}^{1}_{t},q^{1}_{t},\vartheta^{1}_{t},\\\
\\\ \ldots,\\\ \\\
\mathbf{x}^{H}_{t},q^{H}_{t},\vartheta^{H}_{t},z_{t+1}\end{subarray}\right)=\\\
\prod_{h=1}^{H}p(\mathbf{x}_{t+1}^{h}|z_{t+1},\mathbf{x}_{t}^{h})\,p(q_{t+1}^{h}|z_{t+1},q_{t}^{h})\,p(\vartheta_{t+1}^{h}|z_{t+1},\vartheta_{t}^{h}).$
(12)

A graphical model representation of the model is given in Figure 4.

$z_{t-1}$$z_{t}$$z_{t+1}$$\mathbf{x}_{t-1}^{1}$$\mathbf{x}_{t}^{1}$$\mathbf{x}_{t+1}^{1}$$q_{t-1}^{1}$$q_{t}^{1}$$q_{t+1}^{1}$$\vartheta_{t-1}^{1}$$\vartheta_{t}^{1}$$\vartheta_{t+1}^{1}$$\vdots$$\vdots$$\vdots$$\mathbf{x}_{t-1}^{H}$$\mathbf{x}_{t}^{H}$$\mathbf{x}_{t+1}^{H}$$q_{t-1}^{H}$$q_{t}^{H}$$q_{t+1}^{H}$$\vartheta_{t-1}^{H}$$\vartheta_{t}^{H}$$\vartheta_{t+1}^{H}$latentvariable1-st
e.e.$\Huge\vdots$$H$-th e.e.

Figure 4: Graphical representation of a Pose+Gripper ARHMM.

To demonstrate the improvement in the segmentation quality given by the
different dynamics, we compare our model against a Linear ARHMM with the same
topology, so that the only difference lies in the definition of the Auto
Regressive dynamics itself. This means that position, orientation, and gripper
angle for each end effector is independent of each other, and the conditional
probability of the next observed state follows the same law as in (12). The
only difference is that the Linear ARHMM uses linear dynamics for the
evolution of each observed state. Since the conditional independence (and,
thus, the topology) are identical, the graphical representation of the model
does not change.

To test our method, we use the _JIGSAW_ dataset [22]. It consists of three
different surgical tasks (‘SUTuring’, ‘Knot Tying’, and ‘Needle Passing’)
executed by eight surgeons of different skill levels. Data consists of
positions, velocities, orientations (in the form of rotation matrices),
angular velocities, and gripper angle. Each variable is given for both left
and right arms, and for both the patient-side and surgeon-side controllers.
Moreover, JIGSAW provides a segmentation of all the tasks in gestures. We use
only positions, orientations (converting from rotation matrices to unit
quaternions) and gripper angle for the patient-side end-effectors. Thus, in
our case, $H=2$.

To evaluate the quality of the algorithm, we proceed as follows. For a given
batch of data, we use the Expectation Maximization algorithm to infer the set
of parameters $\Theta$ for both models. Then, we apply the Viterbi algorithm
to segment one of the executions that were not in the training set, and we
compute different scores for segmentation. This test is repeated a given
number of times randomizing each time both the training and testing sets.

The scores we decide to adopt are the following: _Seg-score_ and _Silhouette
Index_ (SI). Seg-score (also known as _Jaccard index_) [23, 24] is a
_supervised score_ , that is, it compares the obtained segmentation to the
ground-truth; it is the size of similarity between segmentation result and
ground-truth. SI [25], instead, is an _unsupervised score_ and it evaluates
the goodness of the segmentation by treating it as a clustering problem; it
evaluates the clusters by comparing the average distance within a cluster with
the average distance to the points in the nearest clusters. We decided to use
these scores since they are already used for the evaluation of segmentation
algorithms in the context of robotic surgery [26, 27].

The main challenge in using this dataset is that it is not consistent: the
same task may contains different gestures in each trial (that is, a particular
gesture is present in one trial but not in the others). This makes the
segmentation a particularly hard problem: if the model has many hidden modes,
but some are rare, over-fitting and over-segmentation (i.e. a segment is
wrongly split in multiple smaller parts) will likely happen. On the other
hand, if we use fewer hidden modes, some different gestures will be described
by the same hidden mode, resulting in a lower Seg-score.

(a) Knot Tying - 6 modes, Seg-score. (b) Knot Tying - 6 modes, SI.

(c) Suturing - 5 modes, Seg-score. (d) Suturing - 5 modes, SI.

Figure 5:  Histograms showing the number of occurrences for different scores.
Color blue marks the ‘state of the art’ linear ARHMM, while orange marks our
proposed generalization.

In Figure 5 we show the results of this test. In particular, for each of 100
trials with a given task and a given number of hidden modes, we show the
histograms counting the number of occurrences of the scores in a given
interval. As it can be seen, the NL model is able to achieve, on average, a
higher score in both the supervised and unsupervised metric, and it is also
able to achieve higher scores (see, for instance, Figure 5d).

##  V CONCLUSION AND FUTURE WORKS

In this work, we proposed a generalization to the Auto-Regressive Hidden
Markov Model via modifications of the Auto-Regressive dynamics. In particular,
we proposed a Non-Linear dynamics in Cartesian space and a linear dynamics in
Unit-Quaternion space.

Experiments on real datasets show that adopting these new dynamics result in
an improvement in segmentation scores. This has been proved by comparing two
topologically identical models in which the only difference is the formulation
of the vector fields governing the evolution of the observed state.

As future work, we aim to further generalize the observed state dynamics to
dynamical systems used in trajectory learning for robotics such as _Dynamic
Movement Primitives_ [28, 29] and _Probabilistic Movement Primitives_ [30].
This would allow to simultaneously segment a robot trajectory and extract the
robot movements used to generate the trajectory.

Moreover, we aim at generalizing ARHMMs to deal with dynamics in generic
Riemannian Manifolds and non-Euclidean spaces, possibly extending our proposed
idea for Unit-Quaternions to the space $SO(3)$ of Rotation Matrices and the
space of Symmetric Positive Definite (SPD) matrices, both heavily used in
robotics.

## References

  * [1] S. Calinon, F. D’halluin, E. L. Sauser, D. G. Caldwell, and A. G. Billard, “Learning and reproduction of gestures by imitation,” _IEEE Robotics Automation Magazine_ , vol. 17, no. 2, pp. 44–54, June 2010.
  * [2] I. Visser, “Seven things to remember about hidden markov models: A tutorial on markovian models for time series,” _Journal of Mathematical Psychology_ , vol. 55, no. 6, pp. 403–415, 2011.
  * [3] A. Vakanski, I. Mantegh, A. Irish, and F. Janabi-Sharifi, “Trajectory learning for robot programming by demonstration using hidden markov model and dynamic time warping,” _IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)_ , vol. 42, no. 4, pp. 1039–1052, 2012.
  * [4] B.-H. Juang and L. Rabiner, “Mixture autoregressive hidden markov models for speech signals,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 33, no. 6, pp. 1404–1413, 1985.
  * [5] F. Jelinek, _Statistical methods for speech recognition_. MIT press, 1997.
  * [6] R. Nag, K. Wong, and F. Fallside, “Script recognition using hidden markov models,” in _ICASSP’86. IEEE International Conference on Acoustics, Speech, and Signal Processing_ , vol. 11. IEEE, 1986, pp. 2071–2074.
  * [7] C. Manning and H. Schutze, _Foundations of statistical natural language processing_. MIT press, 1999.
  * [8] C. E. Reiley, H. C. Lin, B. Varadarajan, B. Vagvolgyi, S. Khudanpur, D. D. Yuh, and G. D. Hager, “Automatic recognition of surgical motions using statistical modeling for capturing variability,” in _MMVR_ , vol. 132\. Citeseer, 2008, pp. 396–401.
  * [9] B. Varadarajan, C. Reiley, H. Lin, S. Khudanpur, and G. Hager, “Data-derived models for segmentation with application to surgical assessment and training,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2009, pp. 426–434.
  * [10] L. Tao, E. Elhamifar, S. Khudanpur, G. D. Hager, and R. Vidal, “Sparse hidden markov models for surgical gesture classification and skill evaluation,” in _International conference on information processing in computer-assisted interventions_. Springer, 2012, pp. 167–177.
  * [11] Y. Ephraim, D. Malah, and B.-H. Juang, “On the application of hidden markov models for enhancing noisy speech,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 37, no. 12, pp. 1846–1856, 1989.
  * [12] L. R. Rabiner, _A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition_. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1990.
  * [13] K. P. Murphy, “Switching kalman filters,” 1998.
  * [14] Y. Ephraim and W. J. Roberts, “Revisiting autoregressive hidden markov modeling of speech signals,” _IEEE Signal processing letters_ , vol. 12, no. 2, pp. 166–169, 2005.
  * [15] O. Kroemer, H. Van Hoof, G. Neumann, and J. Peters, “Learning to predict phases of manipulation tasks as hidden states,” in _Robotics and Automation (ICRA), 2014 IEEE International Conference on_. IEEE, 2014, pp. 4009–4014.
  * [16] X. Guan, R. Raich, and W.-K. Wong, “Efficient multi-instance learning for activity recognition from time series data using an auto-regressive hidden markov model,” in _International Conference on Machine Learning_ , 2016, pp. 2330–2339.
  * [17] O. Kroemer, C. Daniel, G. Neumann, H. Van Hoof, and J. Peters, “Towards learning hierarchical skills for multi-phase manipulation tasks,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 1503–1510.
  * [18] S. Chiappa and J. R. Peters, “Movement extraction by detecting dynamics switches and repetitions,” in _Advances in neural information processing systems_ , 2010, pp. 388–396.
  * [19] C. M. Bishop, _Pattern Recognition and Machine Learning (Information Science and Statistics)_. Secaucus, NJ, USA: Springer-Verlag New York, Inc., 2006.
  * [20] A. J. Ijspeert, J. Nakanishi, and S. Schaal, “Learning attractor landscapes for learning motor primitives,” in _Advances in neural information processing systems_ , 2003, pp. 1547–1554.
  * [21] G. D. Forney, “The viterbi algorithm,” _Proceedings of the IEEE_ , vol. 61, no. 3, pp. 268–278, 1973.
  * [22] Y. Gao, S. S. Vedula, C. E. Reiley, N. Ahmidi, B. Varadarajan, H. C. Lin, L. Tao, L. Zappella, B. Béjar, D. D. Yuh _et al._ , “Jhu-isi gesture and skill assessment working set (JIGSAWS): A surgical activity dataset for human motion modeling,” in _MICCAI Workshop: M2CAI_ , vol. 3, 2014, p. 3.
  * [23] G. Ivchenko and S. Honov, “On the Jaccard similarity test,” _Journal of Mathematical Sciences_ , vol. 88, pp. 789–794, 1998.
  * [24] M. Maška, V. Ulman, D. Svoboda, P. Matula, P. Matula, C. Ederra, A. Urbiola, T. España, S. Venkatesan, D. M. Balak _et al._ , “A benchmark for comparison of cell tracking algorithms,” _Bioinformatics_ , vol. 30, no. 11, pp. 1609–1617, 2014.
  * [25] P. J. Rousseeuw, “Silhouettes: a graphical aid to the interpretation and validation of cluster analysis,” _Journal of computational and applied mathematics_ , vol. 20, pp. 53–65, 1987.
  * [26] M. Fard, S. Ameri, R. Chinnam, and R. Ellis, “Soft boundary approach for unsupervised gesture segmentation in robotic-assisted surgery,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 1, pp. 171–178, 2016.
  * [27] A. Murali, A. Garg, S. Krishnan, F. Pokorny, P. Abbeel, T. Darrell, and K. Goldberg, “TSC-DL:Unsupervised trajectory segmentation of multi-modal surgical demonstrations with deep learning,” _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , 2016.
  * [28] A. J. Ijspeert, J. Nakanishi, and S. Schaal, “Movement imitation with nonlinear dynamical systems in humanoid robots,” in _Robotics and Automation, 2002. Proceedings. ICRA’02. IEEE International Conference on_ , vol. 2. IEEE, 2002, pp. 1398–1403.
  * [29] M. Ginesi, N. Sansonetto, and P. Fiorini, “Overcoming some drawbacks of dynamic movement primitives,” _Robotics and Autonomous Systems_ , vol. 144, p. 103844, 2021.
  * [30] A. Paraschos, C. Daniel, J. R. Peters, and G. Neumann, “Probabilistic movement primitives,” in _Advances in neural information processing systems_ , 2013, pp. 2616–2624.


# Large-scale anisotropies of extragalactic cosmic rays
below the ankle

S. Mollerach1, E. Roulet1 and O. Taborda2
1 Centro Atómico Bariloche, Comisión Nacional de Energía Atómica
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
Av. Bustillo 9500, R8402AGP, Bariloche, Argentina
2 Universidad Tecnológica de Pereira, Complejo Educativo La Julita
Pereira, Risaralda, Colombia

###### Abstract

We study the anisotropies on large angular scales which can be present in the
flux of cosmic rays reaching the Earth from a population of extragalactic
sources, focusing on the energy range between the second knee and the ankle.
In this energy range the particles are significantly affected by the Galactic
magnetic field, which then plays a relevant role in shaping the expected
anisotropies. The Galactic magnetic field deflects the cosmic-ray trajectories
and thus modifies the anisotropies present outside the halo of the Galaxy, in
particular the dipolar one associated with the translational motion of the
observer (Compton-Getting effect). Also, due to the Galactic rotation, in the
reference frame of an observer at Earth there is an electric component of the
Galactic field that produces a small change in the particles’ momentum. This
acceleration depends on the cosmic-ray arrival direction and it hence induces
anisotropies in the flux observed in a given energy range. We analyse the
expected amplitude and phase of the resulting dipolar component of the flux
and discuss the possibility to explain via these effects the change in the
phase of the right-ascension distribution which is observed at energies around
1 EeV.

##  1 Introduction

The study of the anisotropies in the arrival directions is expected to play a
fundamental role, together with the spectrum and composition measurements, for
understanding the origin and the propagation of high-energy cosmic rays (CRs).
For energies between $10^{10}$ eV and $10^{20}$ eV the overall CR spectrum has
a shape which is close to a power law, decreasing approximately as $E^{-3}$.
Its main features are the knee, which indicates the steepening at around 3 to
4 PeV, the second knee, corresponding to a further steepening at about 100 PeV
[1, 2], there is then a hardening at the ankle around 5 EeV, followed by a
strong falloff starting at around 50 EeV [3].

A widely accepted explanation for the knee features, also in agreement with
the composition trend from lighter to heavier nuclei which is measured in the
knee to second-knee energy range, is that the spectrum of Galactic protons
starts to fade at the knee while that of the Galactic irons starts to fade at
an energy 26 times higher, as a consequence of a rigidity dependent
phenomenon. Above the second knee and up to about 2 EeV the composition
becomes progressively lighter for increasing energies [4]. There are good
reasons to believe that this population of light cosmic rays, showing up
between the second knee and the ankle, has an extragalactic origin. On the one
hand, no Galactic source is expected to be able to accelerate light elements
up to those energies and, on the other hand, Galactic protons would give rise
to exceedingly large anisotropies in the distribution of arrival directions,
in disagreement with observations [5, 6]. Indeed, no significant anisotropies
have been measured between the second knee and the ankle, with the upper
bounds obtained on the equatorial component of the dipole amplitude being at
the level of 1 to 2% in this energy range. It is anyway worth to note that a
change in the right-ascension phase of the reconstructed dipole is observed,
going from directions close to the Galactic center one, $\alpha\simeq
270^{\circ}$, at energies below about 1 EeV [7, 8], to the nearly opposite
direction, $\alpha\simeq 90^{\circ}$, above 4 EeV [9, 6]. Regarding the
hardening observed at the ankle, it is probably associated with the emergence
of a new extragalactic component of powerful and closeby sources, while the
suppression at the highest energies is probably associated with a combination
of the effects of the maximum achievable acceleration at the sources and the
attenuation effects that the CR suffer during their propagation through
radiation fields.

Extragalactic cosmic rays in the energy range between the second knee and
about $2$ EeV can travel large distances without interacting significantly
with the extragalactic background radiation, with their energies being just
redshifted by the expansion of the Universe. Since they can reach the Earth
from cosmological distances, the distribution of their directions of arrival
to the Milky Way is expected to be close to isotropic, although some
anisotropies could arise from inhomogeneities in the distribution of their
sources, combined with the effects of the diffusion in the intergalactic
magnetic fields. Even in the case of a perfectly isotropic extragalactic
distribution, a dipolar anisotropy would be observed from Earth if there is a
relative motion between the observer and the reference frame in which the CR
distribution looks isotropic. Furthermore, after reaching the Galaxy, the
cosmic rays have to traverse the Galactic magnetic field before arriving to
Earth. This modifies the arrival direction distribution in two ways. On one
side, the trajectories are deflected by the magnetic field, and thus the
angular distribution of the flux reaching the Galaxy gets distorted upon
arrival to Earth. In particular, this will change the direction and amplitude
of a dipolar anisotropy present outside the Galaxy as seen from Earth, and
will also give rise to higher order multipoles. Note that this effect alone
cannot generate an anisotropy if the original flux is isotropic. On the other
side, in the reference frame of an observer at Earth there will be, as a
consequence of the Galactic rotation, an electric component of the Galactic
field, and hence the CRs arriving from extragalactic space will suffer a small
change in their momentum as they travel through the Galaxy [10]. This
acceleration is different for trajectories reaching the Earth from different
directions and it hence induces anisotropies in the flux observed in any given
energy range, even in the case in which the original flux outside the Galaxy
were isotropic. The detailed study of these different anisotropies will be the
subject of the present work, extending the study of ref. [10] to lower
energies, considering updated and more complete magnetic field models and
confronting with recent experimental results.

##  2 The Galactic magnetic field and its effects on the flux of extragalactic
cosmic rays

Let us consider the case in which there is a reference frame for which the
extragalactic CRs are isotropic outside the Galaxy, with a flux
$\Phi^{\prime}(p)$ depending only on the particle’s momentum $p$, but
independent from its direction $\hat{u}$ (where $\vec{p}=p\hat{u}$). The flux
of particles arriving to Earth with momentum $p_{0}$ from the direction
$\hat{u}_{0}$, $\Phi(p_{0},\hat{u}_{0})$, considering an energy spectrum
proportional to $E^{-\gamma}$, is related with that outside the Galaxy through
[10]

$\Phi(p_{0},\hat{u}_{0})\simeq\Phi^{\prime}(p_{0})\left[1+(\gamma+2)\frac{\hat{u}_{h}\cdot\vec{V}_{\rm
iso}}{c}-(\gamma+2)\frac{p_{h}-p_{0}}{p_{0}}\right],$ (1)

as a function of the associated directions ${\hat{u}_{h}}$ and the momenta
$p_{h}$ that the particles had before entering the Galactic magnetic field
region (i.e. just outside the Galactic halo), with $c$ being the speed of
light.

The second term in the parenthesis in Eq. (1) corresponds to the well-known
Compton-Getting effect [11], that leads to a dipolar anisotropy with an
amplitude proportional to the relative velocity $V_{\rm iso}$ of the observer
with respect to the frame in which CRs are isotropic. If this frame coincides
with the rest frame of the cosmic microwave background (CMB), the solar system
relative velocity will be $V_{\rm iso}\simeq V_{\rm CMB}\simeq 371$ km s-1, in
the direction with Galactic coordinates
$(l,b)\simeq(264.4^{\circ},48.4^{\circ})$ [12]. Ignoring the deflections
induced by the Galactic magnetic field (i.e. if $\hat{u}_{0}=\hat{u}_{h}$),
this motion would induce a dipole with amplitude $d\simeq 0.006$ for a
spectral index $\gamma\simeq 3$ [13]. However, the deflections of the
particles’ trajectories lead to a non-trivial mapping between the directions
of incidence at the halo ${\hat{u}_{h}}$ and the arrival directions at Earth
${\hat{u}_{0}}$, distorting the dipolar pattern.

The third term within the parentheses in Eq. (1) corresponds to the effect of
the change of momentum along the trajectory, which has to be evaluated through
the integral [10]

$\frac{p_{h}-p_{0}}{p_{0}}\simeq\frac{0.9\ {\rm EeV}}{E/Z}\int\frac{{\rm
d}\vec{l}}{\rm
kpc}\cdot\left[\frac{(\vec{V}-\vec{V}_{\odot})}{c}\times\frac{\vec{B}}{\mu{\rm
G}}\right].$ (2)

Note that due to the high conductivity of the Galactic medium, we can consider
that in the reference frame moving locally with the plasma there are only
magnetic fields present, with no electric fields. This ionised gas moves with
a velocity $\vec{V}$ which is approximately equal to the one of the rotation
curve of the Galaxy. To an observer moving with the solar system velocity
$\vec{V}_{\odot}$, CRs traveling through far away regions of the Galaxy will
then also experience an electric force coming from the relative motion of the
medium in which the field is purely magnetic, given by
$q\vec{E}=(q/c)\Gamma_{\Delta
V}(\vec{V}-\vec{V}_{\odot})\times\vec{B}^{\prime}$, with $\Gamma_{\Delta
V}\simeq 1$ being the Lorentz factor associated with the relative velocity
between the local system and the observer,
$\Delta\vec{V}\equiv\vec{V}-\vec{V}_{\odot}$, and $\vec{B}^{\prime}$ is the
magnetic field in the local frame (since $\Gamma_{\Delta V}\simeq 1$, we can
also consider that $\vec{B}\simeq\vec{B}^{\prime}$). To estimate
$|\vec{V}-\vec{V}_{\odot}|$ we will just consider the plasma to be moving
azimuthally with a flat rotation curve at 220 km s-1 in a clockwise direction
as seen from the North. Let us note that since the values of
$|\vec{V}-\vec{V}_{\odot}|$ are roughly of the same order of magnitude as
$V_{\rm CMB}$, for a magnetic field of strength $\sim\mu$G and a propagation
length of few kpc within the Galaxy, for values of $E/Z$ smaller than about 1
EeV the momentum change effect could give a contribution of similar or larger
size than the Compton-Getting one. Then, this effect is expected to play a
relevant role in shaping the anisotropies of extragalactic cosmic rays in the
energy range from the second knee to the ankle.

The resulting anisotropies in the arrival directions depend on the precise
model of the Galactic magnetic field, which is however not well known. For
definiteness, we will adopt as reference model the one by Jansson and Farrar
(JF12), which is constrained by Faraday rotation measures and the synchrotron
emission maps from the WMAP satellite [14, 15]. In order to exemplify the
dependence of the results on the magnetic field model adopted, we will also
consider the modifications to this model proposed in [16] (JF+Planck), which
provide a better fit to the synchrotron and dust emission maps determined with
the Planck satellite. In these models there is a coherent regular field, a
random striated (ordered) field and a random isotropic field. The coherent
field is described by a superposition of disk, toroidal halo and “X-field”
components [14]. The random striated field is aligned with the regular field,
has vanishing mean amplitude and root-mean-square (rms) strength $B_{\rm
str}^{2}=\beta|B_{\rm reg}|^{2}$. The random isotropic field consists of a
superposition of disk and halo components, with rms strength $B_{\rm
rand}^{2}=B_{\rm disk}^{2}+B_{\rm halo}^{2}$. For the coherence length of the
random fields we will consider a reference value of $l_{\rm c}=50$ pc [17,
18]. The effect of the striated field is taken into account by considering
$\vec{B}_{\rm str}=b\vec{B}_{\rm reg}$, with the parameter $b$ being changed
randomly after the particles travel a distance $l_{\rm c}$, sampling it from a
Gaussian distribution with variance $\beta$. Note that for the striated
component the JF12 model adopts $\beta=1.36$ while the JF+Planck model
considers $\beta=10$. The effect of the random isotropic field is taken into
account including in the Lorentz equation a stochastic deflection term, as
originally introduced in ref. [19]. This treatment is expected to give a good
estimate of the mean flux distribution over different realizations of the
random field, as well as of the resulting anisotropies on large angular
scales, while for the rigidities that are significantly affected by the random
field the small scale details depend on the actual realization of the random
field deflections that one may consider. Since the different components have
different associated uncertainties, we will separately display the effects
resulting from just the regular magnetic field component, the regular plus
striated components and the complete one including also the isotropic random
contribution, for both magnetic field models.

###  2.1 Anisotropies due to the acceleration by the electric force

Figure 1: Maps in Galactic coordinates of the anisotropies at Earth resulting
from the acceleration of CRs due to the electric force, considering only the
regular Galactic magnetic field and for different rigidities $E/Z$. The left
column corresponds to the magnetic field model JF12, while the right column
corresponds to JF+Planck. The direction of the dipolar component is indicated
by an asterisk and labeled as ‘dip’.

Let us first discuss the anisotropies resulting from the acceleration of the
CRs by the electric force described by the last term in Eq. (1), that can be
computed by integrating Eq. (2) along the different particle trajectories,
considering that the original extragalactic flux is isotropic. We show in Fig.
1 the anisotropy maps obtained by back tracking a regular grid of directions
at Earth up to the outskirts of the Galaxy (taken as a sphere with
galactocentric radius of 20 kpc) for different values of the CR rigidity at
arrival, $E/Z$, and considering a spectral index $\gamma=3.3$, as is measured
below the ankle [3]. In these plots only the regular component of the magnetic
field is included. The plots on the left correspond to the results obtained
using the JF12 magnetic field model, while those on the right correspond to
the results obtained with the modified field JF+Planck. Red regions in the
maps correspond to the arrival directions of trajectories along which cosmic
rays have been accelerated while traversing the Galactic magnetic field, and
thus a larger flux should reach the Earth with the quoted energy as a
consequence of the larger abundance of the particles which had lower energies
outside the Galaxy. Blue regions instead correspond to directions along which
particles have been decelerated, and thus a lower flux should reach the Earth
at that energy given the lower abundance of particles with higher energies
outside the Galaxy. We can see that the shape of the anisotropies depends on
the magnetic field model considered while the overall amplitudes of the
resulting anisotropies are comparable in both models.

At the highest value of $E/Z$ considered of $10^{19}$ eV, for which the
trajectories are less deflected by the magnetic field, we can identify more
easily the main magnetic field component that is responsible for the
anisotropy pattern observed. Since both the relative velocity
$\Delta\vec{V}\equiv\vec{V}-\vec{V}_{\odot}$ as well as the disk and toroidal
magnetic field components are parallel to the Galactic plane, these components
give rise to an electric field which is orthogonal to the Galactic plane.
Thus, particles reaching the Earth along directions close to the Galactic
plane are not significantly accelerated by these components, and the X-field
component is the one giving the dominant contribution, having the effect of
accelerating the particles coming through the inner Galactic directions, both
above and below the Galactic plane. At higher Galactic latitudes, both towards
the North and the South of the Galactic plane, it is the toroidal halo
component the one giving the largest contribution, involving a sign reversal
of the effect from positive to negative longitudes. On the other hand, in this
case there is no sign reversal from the North with respect to the South
because the toroidal field component has opposite directions in the two
hemispheres. For lower rigidities the CR trajectories experience larger
deflections, and hence the interpretation of these effects becomes less
direct.

Including the random striated and random isotropic magnetic field components
has in both cases the effect of smoothing the features that were present in
the maps including only the regular field. We show in Fig. 2 the effect of the
random field components for both the JF12 and JF+Planck, for a rigidity of
$10^{17.8}$ eV. As in Fig. 1, the left plots correspond to the JF12 model and
the right ones to the JF+Planck model, while the top row shows the results
including the regular and random striated components, and the bottom row the
complete field including also the isotropic random component. The details of
the maps actually depend on the particular realization of the random field
deflections that were considered, and the differences between different
realizations would become more significant when considering smaller
rigidities, for which the deflections are larger.

Figure 2: Maps in Galactic coordinates of the anisotropy at Earth from the
acceleration of CRs due to the electric force, considering the regular and
random striated magnetic field components (top row) and including also the
random isotropic component (bottom row), for a rigidity $E/Z=10^{17.8}$ eV.
The left column corresponds to the magnetic field model JF12 while the right
column corresponds to JF+Planck.

Of particular interest is the contribution due to the acceleration by the
electric force to the dipolar component of the flux distribution at Earth. The
direction of the associated dipole is shown with an asterisk in the examples
of Figs. 1 and 2 (and is labeled as ‘dip’). In Fig. 3 we show the amplitude of
the dipolar component as a function of the rigidity $E/Z$, obtained when
considering only the regular field, the regular and the random striated
components and the complete field including also the random isotropic
component. For the cases with random fields the colored bands correspond to
the dispersion around the mean values obtained from ten different simulations
for each energy. As already mentioned, for both models the dipolar amplitude
due to the CR acceleration tends to increase for decreasing rigidity, while on
the other hand the inclusion of the random field leads, due to the smearing
induced, to a reduction of the dipolar component. The amplitude resulting from
the regular magnetic field model is slightly smaller for the JF+Planck model,
mostly as a consequence of the reduced magnitude of the halo field. The effect
of the striated field on the dipolar component is more notorious in the
JF+Planck model due to the larger associated value of the factor $\beta$, and
it reduces the amplitude of the dipolar component at low energies.

Figure 4 shows with dots the mean dipole direction obtained in several
simulations with different realizations of the random deflections, and also
shows the circular sky region around it within which the dipole points in half
of the simulations, for different rigidities and considering the complete
(regular, striated and random isotropic) field. At the highest energies, for
which the dominant effect is that of the regular field, a small dispersion in
the directions (as well as in the amplitudes shown in Fig. 3) appears. The
trend is similar for the JF12 model, displayed in the left panel, and for the
JF+Planck one, displayed in the right panel. As the energy decreases and the
effect of the random magnetic field becomes relatively more important, the
dispersion in the dipolar directions increases, being this more notorious for
the JF+Planck model due to the larger amplitude of the halo random isotropic
component in this case. In particular, at the lowest energy displayed, the
random isotropic component is dominant for this model, resulting in an almost
random distribution of the dipolar directions in the different simulations.

Figure 3: Dipole amplitude as a function of the logarithm of the rigidity
$E/Z$. Left panel corresponds to the JF12 magnetic field model while the right
panel corresponds to JF+Planck.

Figure 4: Maps in Galactic coordinates showing the average dipole direction
for different values of log($E/Z/$eV) when the complete magnetic field
(regular + striated + random isotropic components) is considered. The circular
contours indicate the median of the distribution of the dipole directions.
Left panel corresponds to the JF12 magnetic field model, right panel
corresponds to JF+Planck.

From the plots in Figs. 1 and 2 it is apparent that the anisotropies have not
just a dipolar component, but that higher order multipole moments are also
significant. By expanding the flux reaching the Earth in terms of spherical
harmonics with coefficients $a_{\ell m}$, such that111We use the real
spherical harmonic basis, so that the $a_{\ell m}$ are real.

$\Phi(\hat{u})=\sum^{\infty}_{l=0}\ \sum^{\ell}_{m=-\ell}\ a_{\ell m}\ Y_{\ell
m}(\hat{u}),$ (3)

we can quantify the amplitude of the different multipoles through the angular
power spectrum $C_{\ell}$, which is obtained by averaging the $a_{\ell m}^{2}$
coefficients over the $2\ell+1$ possible values of $m$

$C_{\ell}=\frac{1}{2\ell+1}\ \sum_{m=-\ell}^{\ell}\ a_{\ell m}^{2}.$ (4)

Figure 5: Normalized multipole amplitude as a function of the rigidity for the
three lowest $\ell$ values of the anisotropies resulting from the acceleration
due to the electric force. The left panel corresponds to the regular component
of the JF12 magnetic field model, while the right panel corresponds to
JF+Planck.

We show in Fig. 5 the normalized contribution $\sqrt{\ell(\ell+1)C_{\ell}}$
for the three lowest order terms, which are the dipole, the quadrupole and the
octupole, as a function of the rigidity, when considering the regular magnetic
field component of the JF12 and JF+Planck models. The amplitude of all of them
are of similar size and they decrease with increasing energy, ranging from
${\rm few}\times 10^{-3}$ for rigidities close to $10^{17}$ eV down to $\sim
10^{-4}$ for rigidities above $10^{19}$ eV.

###  2.2 Modification of the anisotropies outside the Galaxy due to Galactic
magnetic field deflections

In addition to the anisotropies due to the acceleration by the electric force
just discussed, other sources of anisotropies of extragalactic CRs could be
the Compton-Getting effect resulting from the motion of the observer with
respect to the frame in which CRs are isotropic (second term in Eq. (1)) or
eventually some intrinsic anisotropies may already be present in the flux
reaching the halo of the Galaxy.

Setting aside possible intrinsic anisotropies and in the absence of the
Galactic magnetic field (so that $\hat{u}_{h}=\hat{u}_{0}$), the Compton-
Getting effect would give rise to a dipolar anisotropy pointing in the
direction of the velocity of the observer with respect to the isotropic frame.
For definiteness we will exemplify the distortion of this pattern resulting
from the propagation of CRs through the Galactic magnetic field for the case
in which the extragalactic CRs are isotropic in the rest frame of the CMB.
Figure 6 shows the resulting anisotropy maps at the Earth for different
rigidities, obtained considering the effects of the regular magnetic field
component of the JF12 model. At the highest rigidity displayed, the overall
distribution is close to dipolar, pointing in a direction close to the CMB
dipole one and with an amplitude $d\simeq 0.0064$. For decreasing energies the
amplitude of the dipolar component decreases to $d\simeq 0.0050$ for
$E/Z=10^{18.4}$ eV and $d\simeq 0.0020$ for $E/Z=10^{17.8}$ eV, while the
direction moves down towards the Galactic plane. At the smallest rigidity
displayed, $E/Z=10^{17.2}$ eV, the amplitude grows again to $d\simeq 0.0056$
and the direction flips to the almost opposite longitude.

Figure 6: Maps in Galactic coordinates of the anisotropies resulting from the
Compton-Getting effect, for different values of the rigidity, assuming that
the CRs are isotropic in the CMB rest frame and considering the regular
Galactic magnetic field model of JF12. The direction towards the CMB dipole is
indicated with a plus sign and that of the resulting CR dipole with a star.

The random striated and random isotropic components have the effect of
smoothing down the anisotropies and of reducing the dipolar amplitude. We show
in Fig. 7 the resulting Compton-Getting anisotropies when including the random
striated component (left panel) and for the complete field including also the
random isotropic component (right panel), for the JF12 model and for a
rigidity of $E/Z=10^{17.8}$ eV. Note that the amplitude of the dipolar
component drops to $d\simeq 0.0016$ when including the striated field and
$d\simeq 7\times 10^{-4}$ for the complete field case.222For simplicity we do
not display the maps for the JF+Planck model, but the features are
qualitatively similar to those shown for the JF12 model.

Figure 7: Same as Fig. 6 for a rigidity $E/Z=10^{17.8}$ eV, considering the
regular and striated magnetic field components (left panel) and the complete
field including also the isotropic random component (right panel) of the JF12
model.

Figure 8: Dipole amplitude resulting from the Compton-Getting effect after
traversing the Galactic magnetic field as a function of log($E/Z)$. The left
panel corresponds to the magnetic field JF12 while the right panel corresponds
to JF+Planck.

Figure 9: Maps in Galactic coordinates showing the average dipole direction
from the Compton Getting effect for different values of log($E/Z/{\rm eV})$
when the complete magnetic field is considered. The circular contours indicate
the median of the simulated dipole directions. Left panel corresponds to the
magnetic field JF12, right panel corresponds to JF+Planck.

In Fig. 8 we show the amplitude of the dipolar component as a function of the
rigidity $E/Z$, obtained when considering only the regular field, the regular
and random striated components and the complete field including also the
random isotropic component. For the cases including random fields, the colored
bands correspond to the dispersion around the mean values of ten simulations
performed for each energy. At the highest energies the amplitude is close to
that expected in the absence of magnetic field and for $\gamma=3.3$, which is
$d\simeq 0.0066$, and then decreases for decreasing energies. When only the
regular field or regular plus striated fields are included, the amplitude
rises at the lowest rigidities considered. The inclusion of the random
isotropic component leads to a smaller dipolar component, especially at the
lowest rigidities considered, for both the JF12 (left panel) and JF+Planck
(right panel) models.

Figure 9 shows the average dipole direction as well as the circular region
around it within which the dipolar component points in half of the
simulations, for different rigidities and for the case where the complete
field (regular, striated and random isotropic) is considered. At the highest
energy the dipole points close to the CMB dipole direction, while when the
energy decreases it moves down towards the Galactic plane, and then at the
lowest energies it moves up again, with an increasing dispersion being
obtained between the different realizations of the random field for decreasing
energies.

Figure 10: Normalized multipole amplitude as a function of the rigidity for
the three lowest $\ell$ values, assuming that the CRs are isotropic in the CMB
rest frame. Left panel is for the regular component of the JF12 model and
right panel for that of the JF+Planck model.

The deflections produced by the Galactic magnetic field not only change the
direction and magnitude of the dipole present outside the halo, but they also
produce significant distortions in the angular distribution pattern of the
flux, as is clearly seen in Figs. 6 and 7. In order to quantify this effect,
we show in Fig. 10 the normalized multipoles of order $\ell=1,\,2$ and 3 as a
function of the rigidity and considering just the regular component of the
JF12 or JF+Planck models. The amplitude of the dipolar component results
dominant with respect to the quadrupolar and octupolar ones both at the
highest and lowest rigidities considered, while for rigidities around
$E/Z=10^{17.5}$ eV the quadrupolar and octupolar amplitudes become comparable
to, or even larger than, the dipolar one. This rigidity range is where the
flip in the direction of the dipolar component occurs.

Figure 11: Direction in Galactic coordinates of the dipolar component observed
at the Earth for different dipolar directions outside the Galactic halo, after
propagation in the regular component of the JF12 (left panel) and the
JF+Planck (right panel) magnetic field models. Black circles indicate the
original dipole directions considered while the intermediate points along the
lines indicate, in the upper panels, the direction at Earth for rigidities
$E/Z=10^{19.3}$ eV, $10^{19}$ eV, $10^{18.7}$ eV and $10^{18.4}$ eV, this last
corresponding to the tip of the arrow. In the lower panels the intermediate
point corresponds to $E/Z=10^{17.2}$ eV, while the tip corresponds to
$E/Z=10^{16.9}$ eV. The color indicates the relative change in the amplitude
of the dipole at Earth with respect to that outside the halo.

If the reference frame in which the CRs are isotropic were different from the
CMB one, the solar system relative velocity $\vec{V}_{\rm iso}$ in Eq. (1)
would be different. The initial dipolar anisotropy, before being distorted by
the magnetic field deflections, would in this case thus point towards a
different direction and would have a different amplitude. The modification of
the direction and amplitude of the dipole by the effects of the Galactic
magnetic field depends on its initial direction. We summarize the effect
produced by the magnetic deflections, for different CR rigidities and initial
directions of the maximum of the flux outside the halo, in Fig. 11. The open
circles correspond to the set of the original dipole directions considered,
and the successive points along the arrows indicate the direction of the
corresponding dipolar component of the flux reaching the Earth for different
CR rigidities, as described in the caption. For rigidities above ${\rm
few}\times 10^{18}$ eV, the dipoles with originally positive Galactic
longitudes tend to be rotated towards directions closer to the inner spiral
arm, at $(\ell,b)\simeq(80^{\circ},0^{\circ}$), while those with original
negative longitudes tend to be rotated towards the outer spiral arm direction,
at $(\ell,b)\simeq(-100^{\circ},0^{\circ})$, as is shown in the upper panels
in Fig. 11. On the other hand, for rigidities below ${\rm few}\times 10^{17}$
eV, the dipoles with originally positive Galactic latitudes tend to be rotated
towards directions closer to the inner spiral arm, while those with originally
negative latitudes tend to be rotated towards the outer spiral arm direction,
as is shown in the lower panels of the figure. We then see that the change in
the dipolar phase from negative Galactic longitudes at high energies to
positive ones at low energies, obtained before for the case in which the
dipolar direction coincides with the CMB one, is actually expected to happen
more generically for original dipole directions roughly in the quadrant with
negative longitudes and positive latitudes. On the other hand, no drastic
change of phase is expected to appear for original dipole directions lying in
the quadrants with positive latitudes and longitudes or negative latitudes and
longitudes. Finally, if the original dipole direction were in the quadrant
with positive longitudes and negative latitudes, the dipole at Earth is
expected to display a change from directions closer to the inner spiral arm at
high energies to directions closer to the outer spiral arm at low energies.

If the CR flux reaching the Galactic halo were to have a significant intrinsic
dipolar anisotropy, as for example resulting from the diffusion in the
extragalactic magnetic field of CRs from some nearby powerful sources, this
anisotropy would be affected by the deflections in the Galactic magnetic field
in the same way as was just discussed for a Compton-Getting dipole. As a
result, the direction and the amplitude of the dipole observed at the Earth
would differ from the values outside the halo in a way similar to what is
shown in Fig. 11.

##  3 Combined results

In the previous sections we have discussed the effects of the Galactic
magnetic field upon the distribution of extragalactic CR arrival directions,
considering separately the acceleration due to the electric force arising from
the Galactic rotation and the modification of a Compton-Getting term due to
magnetic deflections. The former effect is expected to give the dominant
contribution at the lower rigidities. The latter one depends on the relative
velocity of the observer with respect to the reference frame in which CRs are
approximately isotropic, which is however unknown, but for definiteness we
will consider in this discussion the plausible case described in detail before
in which this frame coincides with the CMB one.

We show in Fig. 12 the total dipolar amplitude resulting from the combination
of both effects, as a function of $E/Z$, for the cases in which just the
regular magnetic field component is considered, the case in which the striated
field is also included and for the complete magnetic field with the random
isotropic component, both for the JF12 and JF+Planck models. At the highest
energies the amplitude tends to the value expected from the Compton-Getting
effect, $d\simeq 0.0066$, and it then decreases for decreasing energies as was
discussed in Section 2.2. However, when looking at the regular or regular plus
random striated cases, we see that for rigidities below $E/Z\sim 10^{18}$ eV
(for the JF12 model) or $E/Z\sim 10^{17.5}$ eV (for the JF+Planck model) the
amplitude rises again for decreasing energies. This rise in the amplitude
results from the acceleration due to the electric force. As was discussed in
Section 2.1, the acceleration effect is larger in the JF12 model than in the
JF+Planck and thus the crossover of the two effects occurs at higher
rigidities in the JF12 case. When the random isotropic magnetic field
component is also included, the amplitude for the lower rigidities gets
significantly suppressed.

Figure 13 shows the corresponding dipole directions in the sky for the case of
the complete field and for both magnetic field models. A significant change in
the dipole direction is apparent, specially for the JF12 model, from positive
Galactic longitudes for rigidities below $\sim 1$ EeV to negative Galactic
longitudes for larger rigidities, corresponding respectively to the
predominance of the acceleration or the Compton-Getting effects.

Figure 12: Total dipole amplitude resulting from both the acceleration due to
electric field associated with the Galactic rotation and the Compton-Getting
effect, after traversing the Galactic magnetic field and as a function of the
logarithm of the rigidity $E/Z$. The left panel corresponds to the magnetic
field model JF12, the right panel corresponds to JF+Planck model.

Figure 13: Maps in Galactic coordinates showing the average total dipole
direction for different values of log($E/Z$/eV) when the complete magnetic
field is considered. The circular contours indicate the median of the dipole’s
direction distribution. Left panel corresponds to the magnetic field JF12,
right panel corresponds to JF+Planck.

Since the anisotropies resulting from the two effects discussed in this paper
depend on the rigidity of the particles, $E/Z$, if the flux of particles
reaching the Earth has a mixed composition the anisotropy at a given energy
will result from the superposition of the anisotropies corresponding to the
different rigidities associated with each mass component. The Pierre Auger
Observatory measurements indicate that the composition for energies above
$10^{17.2}$ eV has indeed a mixed nature [20]. In particular, for energies
below the ankle a mixed composition of protons and intermediate mass nuclei is
inferred. We will hence consider, as a simplified approximation, that in this
energy range there is a constant fraction of 50% of H and 50% of N, both of
extragalactic origin, with a negligible contribution from the Galactic
component.

For a given mixture of elements in the flux of cosmic rays reaching the
Galaxy, we can compute the expected dipolar amplitude of the anisotropy at the
Earth for each energy $E$ by adding the dipole vectors for the corresponding
rigidities $E/Z$ weighted by the fraction of particles with each charge $Z$.
We show in Fig. 14 the total dipolar amplitude as a function of the energy in
the mixed composition model for the JF12 (left panel) and for the JF+Planck
(right panel) models. For high energies the amplitude tends to the Compton-
Getting value ($d\simeq 0.0066$) in both models, since CRs are not much
affected by the Galactic magnetic field. A more noticeable difference between
the two models appears at low energies, especially when looking at the regular
field component, with a larger amplitude being obtained for the JF12 model.
This is due to a partial cancellation at these energies between the
anisotropies of the H and N components in the JF+Planck model.

One should also note that at lower energies a transition to a regime dominated
by the Galactic CRs should start to take place, and hence approaching the
second-knee energy the anisotropies should progressively tend to those of the
Galactic component, not considered here. The fact that the effects of the
random magnetic fields become more notorious at low rigidities is also
indicating that a transition to a diffusive regime inside the Galaxy is taking
place for the incident extragalactic CRs.

Figure 14: Total dipole amplitude as a function of the energy in a mixed
composition model (50% of H and 50% of N) of extragalactic cosmic rays
assuming that they are isotropic in the CMB rest frame. Left panel corresponds
to the magnetic field JF12, right panel corresponds to JF+Planck.

##  4 Discussion and conclusions

The large angular scale modulation of the flux of cosmic rays can be better
measured as a function of right ascension (rather than declination), since
taking advantage of the uniform rotation of the Earth the experimental
spurious effects are better controlled. From the observed modulation of the
flux as a function of right ascension it is possible to estimate the
projection of the dipole in the plane orthogonal to the Earth rotation axis,
$d_{\perp}$, and the phase of the maximum of this modulation corresponds to
the right ascension of the dipole direction. We will then shift here to
consider equatorial coordinates so as to allow for a better comparison with
observations.

Figure 15:  Amplitude of the equatorial dipole (upper panels) and its right-
ascension phase (lower panels) in a mixed composition model (50$\%$ of H and
50$\%$ of N). Results are shown for the regular, regular + striated and
complete magnetic field scenarios. The points with error bars are the
measurements from the Pierre Auger Observatory [6]. For energies below 8 EeV,
where the measurements have chance probability larger than 1%, the downward
triangles in the $d_{\perp}$ plots indicate the 99% CL upper limits obtained.
Left panels correspond to the magnetic field JF12, right panels correspond to
JF+Planck.

We show in Fig. 15 the amplitude of $d_{\perp}$ and its right-ascension phase
as a function of the energy, as expected for the JF12 (left panel) and the
JF+Planck (right panel) models and considering only the regular component, the
regular plus the striated components or the complete field with also the
random isotropic component. We also included in the plot the results of the
measurements from the Pierre Auger Observatory [6].

The observed phase shows a transition from ${\rm RA}\simeq-90^{\circ}$ for
energies below 1 EeV to ${\rm RA}\simeq 90^{\circ}$ for energies above 4 EeV.
Note that a very similar shift in the phase as a function of energy is
obtained from the effects of the regular or the regular plus striated
components of the JF12 model. However, the addition of the isotropic random
component somewhat spoils this good agreement between model predictions and
observations. It should be mentioned however that several authors have pointed
out that the comparison of the predictions of the original JF12 model for the
synchrotron emission maps with the final WMAP data release, and also with data
from Planck, indicates that a significant reduction of the assumed amplitude
of the isotropic random component is required [18, 23]. A reduction in the
isotropic random field amplitude is on the other hand expected to lead to
predictions for the dipolar components in between those obtained for the
regular and for the complete field models that were shown in Fig. 15, and
hence should improve the agreement with the data. One can also see that the
transition in the phase occurs for lower energies in the JF+Planck model than
in the JF12 one.

In the energy bins below 8 EeV no significant measurement has been obtained
for the equatorial dipole amplitude $d_{\perp}$, and the 99% CL upper bounds
on the amplitude are also displayed in the plot for those bins. For both
magnetic field models considered the predicted amplitude is compatible with
the observations. For energies around 10 EeV a significant dipolar component
has been measured by the Pierre Auger Collaboration, with an amplitude of
about $d\simeq 7\%$ [9, 22], increasing approximately linearly with energy
above 4 EeV and up to the highest energy bin considered [21, 22]. This
anisotropy is expected to be associated with the sources of the highest energy
CRs that dominate the flux for energies above the ankle. Since at the highest
energies CRs are limited to arrive to the Earth from relatively closeby
sources, lying at distances smaller than about 100 Mpc, a plausible
explanation is that the observed anisotropy arises from the inhomogeneous
distribution of the CR sources in the local Universe. The distribution of
matter around us, as traced by galaxies in the 2MRS catalog, does in fact show
a significant dipolar component with a maximum not far from the CMB dipole
maximum [24]. This is in fact expected as it is the gravitational attraction
from the nearby matter what produces the Local Group peculiar velocity that is
ultimately responsible for the CMB dipole. The anisotropy of this high-energy
component is thus expected to have a direction not very different from the
Compton-Getting one but with a larger amplitude, making it the dominant
contribution above the ankle. This could in principle allow to also account
for the measurements at the two highest energy points shown in Fig. 15.

A significant uncertainty in this study is due to the poor knowledge of the
Galactic magnetic field. Different observations, like Faraday rotation
measures (RM) and the intensity of total and polarized diffuse synchrotron
emission, are used to constrain different magnetic field components.333As
described in [25], the regular field component contributes to both polarized
emission as well as to RM, while the isotropic random component contributes
only to the total intensity of the diffuse emission. The striated, or ‘random
ordered’, component contributes to the polarized emission but not to RM.
These observables depend on different projections of the magnetic field, along
or perpendicular to the line of sight, and other not well known quantities
like the electron density along the line of sight, what leads to different
uncertainties affecting the coherent, random striated and random isotropic
components (for a state of the art review of Galactic magnetic field modelling
see [26]). For this reason, we have presented separately the results due to
the regular component alone, for the regular and the striated components and
for the complete field including also the random isotropic component.

As a summary, we have explored in this paper the effects of the Galactic
magnetic field on the large angular scale distribution of the flux of
extragalactic cosmic rays reaching Earth, focusing on the energy range between
the second knee and the ankle. In this energy range anisotropies leading to
dipolar amplitudes at the few per mille level are predicted both from the
Compton-Getting effect and from the acceleration due to the electric force
induced by the rotation of the Galaxy. These amplitudes are consistent with
the existing bounds on the equatorial component of the dipole amplitude below
the ankle energy, and may lead to a change in the right-ascension phase in a
way compatible with the observations, depending on the details of the magnetic
field model. This could also imply that the transition in the phase at
energies of order 1 EeV suggested from observations may not be directly
related to the transition from a Galactic to an extragalactic predominance,
which may be taking place actually at lower energies, closer to the second
knee.

More precise measurements of the anisotropies in this energy range should give
precious input for a better understanding of the different contributions to
the CR fluxes and should also help to study the different processes
originating them.

## Acknowledgments

We are grateful to I. Allekotte for useful comments and to D. Harari for
discussions. This work was supported by CONICET (PIP 2015-0369 and PIP
2021-0565).

## References

  * [1] M.G. Aartsen et al. (IceCube Collaboration), Phys. Rev. D 100 (2019) 310082002
  * [2] D. Ivanov et al. (Telescope Array Collaboration), PoS ICRC2019 (2020) 298
  * [3] P. Abreu et al. (Pierre Auger Collaboration), Eur. Phys. J. C 81 (2021) 966; V. Novotný for the Pierre Auger Collaboration, PoS ICRC2021 (2021) 324
  * [4] A. Yushkov for the Pierre Auger Collaboration, PoS ICRC2019 (2020) 482
  * [5] P. Abreu et al. (Pierre Auger Collaboration), ApJL 762 (2013) L13
  * [6] A. Aab et al. (Pierre Auger Collaboration), ApJ 891 (2020) 142
  * [7] M.G. Aartsen et al. (IceCube Collaboration), ApJ 826 (2016) 220
  * [8] A. Chiavassa et al. (KASCADE-Grande Collaboration), PoS (ICRC2015) 281
  * [9] A. Aab et al. (Pierre Auger Collaboration), Science 357 (2017) 1266
  * [10] D. Harari, S. Mollerach and E. Roulet, JCAP 11 (2010) 033
  * [11] A.H. Compton and I.A. Getting, Phys. Rev. 47 (1935) 817
  * [12] D.J. Fixsen et al., ApJ 473 (1996) 576
  * [13] M. Kachelriess and P.D. Serpico, Phys. Lett. B 640 (2006) 225
  * [14] R. Jansson and G.R. Farrar, ApJ 757 (2012) 14
  * [15] R. Jansson and G.R. Farrar, ApJL 761 (2012) L11
  * [16] R. Adam et al. (Planck Collaboration), A$\&$A 596 (2016) A103
  * [17] M. Haverkorn, J.C. Brown, B.M. Gaensler and N.M. McClure-Griffiths, ApJ 680 (2008) 362
  * [18] M.C. Beck, A.M. Beck, R. Beck, K. Dolag, A.W. Strong and P. Nielaba, JCAP 05 (2016) 056
  * [19] A. Achterberg et al., astro-ph/9907060
  * [20] J. Bellido for the Pierre Auger Collaboration, PoS (ICRC2017) 506
  * [21] A. Aab et al. (Pierre Auger Collaboration), ApJ 868 (2018) 4
  * [22] R. Menezes for the Pierre Auger Collaboration, PoS (ICRC2021) 335
  * [23] M. Unger and G.R. Farrar, EPJ Web of Conferences 210 (2019) 04005, Proceeding of UHECR2018
  * [24] P. Erdogdu et al., MNRAS 368 (2006) 1515
  * [25] T.R. Jaffe, J.P. Leahy, A.J. Banday, S.M. Leach, S.R. Lowe and A. Wilkinson, MNRAS 401 (2010) 1013
  * [26] T.R. Jaffe, Galaxies 7 (2019) 52