A Bag of Receptive Fields for Time Series Extrinsic Predictions

Francesco Spinnato Riccardo Guidotti Anna Monreale Mirco Nanni

Scuola Normale Superiore, Pisa, Italy University of Pisa, Pisa, Italy
University of Pisa, Pisa, Italy ISTI-CNR, Pisa, Italy

###### Abstract

High-dimensional time series data poses challenges due to its dynamic nature,
varying lengths, and presence of missing values. This kind of data requires
extensive preprocessing, limiting the applicability of existing Time Series
Classification and Time Series Extrinsic Regression techniques. For this
reason, we propose BORF, a Bag-Of-Receptive-Fields model, which incorporates
notions from time series convolution and 1D-SAX to handle univariate and
multivariate time series with varying lengths and missing values. We evaluate
BORF on Time Series Classification and Time Series Extrinsic Regression tasks
using the full UEA and UCR repositories, demonstrating its competitive
performance against state-of-the-art methods. Finally, we outline how this
representation can naturally provide saliency and feature-based explanations.

##  1 Introduction

High-dimensional time series data is becoming more widely available to
decision-makers, domain experts, and researchers. This kind of data, by
definition, captures dynamic changes over time and is used in many fields,
including finance, healthcare, and environmental science (Shumway et al.,,
2000). Sensors often record such data at different frequencies or for varying
lengths of time, resulting in signals with a different number of observations.
Sensors can also fail, leaving gaps and missing values in time series signals.
Because of these factors, real-world time series data is frequently “dirty”
and difficult to handle, necessitating extensive preprocessing (Kreindler and
Lumsden,, 2006). These complexities present significant challenges to most
existing time series analysis techniques, frequently developed and tested on
very clean datasets, limiting their applicability in real-world scenarios
(Bagnall et al.,, 2017; Ruiz et al.,, 2021).

Within time series analysis (TSA), Time Series Classification (TSC), and Time
Series Extrinsic Regression (TSER) are extremely relevant tasks. TSC is the
task of predicting a categorical value from a time series, for example,
distinguishing between a healthy and unhealthy patient given vital signals
(Bagnall et al.,, 2017). TSER is instead the task of predicting a scalar,
i.e., a continuous value (Tan et al.,, 2021). TSER differs from time series
forecasting (TSF), also called time series regression (TSR). The term
“extrinsic” emphasizes that, in TSER, contrary to TSF (or TSR), the predicted
value may not be the future of the analyzed series (Tan et al.,, 2021).
Research into TSER has received much less attention in the time series
research community, even if it models interesting problems (Tan et al.,,
2020), such as predicting energy consumption by monitoring rooms in houses
(Hebrail and Berard,, 2012). Further, TSER and TSC approaches should be as
interpretable as possible to ensure applicability in critical domains. In
contrast, state-of-the-art models for TSA are mainly black-box models
(Theissler et al.,, 2022). Additionally, when the number of observations in a
time series dataset is irregular or there are missing values, these models
require several preprocessing steps.

Given these challenges, we propose borf, a Bag-Of-Receptive-Fields. borf is a
deterministic time series transform that converts any time series dataset into
an easier-to-handle, understandable tabular representation. borf is very
flexible in the input type, allowing the processing of real-world multivariate
time series with irregular signals of different lengths and containing missing
values. This is achieved by processing signals independently and generalizing
pattern discretization through complete case analysis. Further, to ensure an
expressive representation, we introduce convolutional operators to extract
receptive fields from the time series. Unlike standard subsequences, receptive
fields can better describe the development of a time series at different
resolutions and highlight its local and global characteristics. We show that
the novel connection between convolutional operators and multi-resolution
symbolic patterns allows for great performance in TSC and TSER. We test borf
on the full UEA and UCR repositories, comprising $177$ datasets, against
several state-of-the-art approaches, $12$ for TSC and $6$ for TSER. We show
that borf challenges competitor models in TSC and outperforms them in TSER.
Finally, we propose a scalable heuristic in time and space complexity and
define a connection between feature and saliency-based explanations in a borf
representation, outlining its natural interpretability.

The rest of the paper is organized as follows. Sec. 2 discusses the related
work and main competitors, while Sec. 3 summarizes the background of our
proposal. borf is presented in Sec. 4, and tested in Sec. 5. Finally, we
present our conclusions in Sec. 6.

##  2 Related Work

The field of time series predictors encompasses a wide range of approaches
(Ruiz et al.,, 2021; Tan et al.,, 2021). Our proposal belongs to the
dictionary-based category, which extracts features from time series by
recording characteristics of discretized subpatterns (Bagnall et al.,, 2017).
These methods segment time series into subsequences, convert them into
symbolic words, and create histograms of feature counts (Baydogan et al.,,
2013). For this reason, they rely heavily on effective time series
approximation methods such as the Symbolic Aggregate approXimation (SAX) (Lin
et al.,, 2007) and Symbolic Fourier Approximation (SFA) (Schäfer and
Högqvist,, 2012). SAX uses Piecewise Aggregate Approximation (PAA) (Keogh et
al.,, 2001) for segmentation and binning, while SFA focuses on spectral
properties using Fourier coefficients (Schäfer and Högqvist,, 2012). Although
informative for certain applications, SFA loses temporal information and has
higher computational complexity than SAX (Schäfer and Högqvist,, 2012).
Notably, 1D-SAX was introduced to enhance SAX representation by using
Piecewise Linear Approximation (PLA) (Hung and Anh,, 2007) instead of PAA, to
infer segment slopes (Malinowski et al.,, 2013), but it has yet to be adopted
in dictionary-based approaches.

In general, even if subsequence extraction can be implemented in supervised
and unsupervised ways, it is mainly applied to TSC. The first dictionary-based
method is the traditional Bag-of-Patterns (bop) (Lin et al.,, 2012), which
uses SAX to generate discretized subsequences, counting their frequency. Other
methods rely on supervised techniques by introducing a tf-idf weighing of
features (sax-vsm) (Senin and Malinchik,, 2013), or by filtering subsequences
extracted with SAX and SFA through a Chi-squared bound (mr-sqm) (Nguyen and
Ifrim,, 2022). boss (Bag-of-SFA-Symbols) (Schäfer,, 2015) and its more recent
extension, weasel+muse (Schäfer and Leser,, 2017), focus on the sole usage of
SFA and achieve state-of-the-art performance in multivariate time series
classification. Being specifically developed for TSC, the main limitation of
these classifiers is that they cannot be easily extended to other tasks.
Specifically, to the best of our knowledge, none was ever tested on TSER.
Unlike most of the previously mentioned approaches, borf is unsupervised and,
to the best of our knowledge, is the first dictionary-based method capable of
effectively tackling both TSC and TSER tasks. We achieve superior performance
by proposing theoretical and practical extensions, introducing several
generalizations addressing the primary limitations of bop. Specifically, borf
allows for extracting multi-resolution 1D-SAX patterns from univariate and
multivariate time series data, accommodating signals of varying lengths and
those with missing values. Dictionary-based features, i.e., time series
patterns, seem particularly well-suited for explaining prediction outcomes via
feature importance or saliency maps. Yet, the interpretability of such
approaches is rarely evaluated (Theissler et al.,, 2022). For instance, in
(Spinnato et al.,, 2022), SAX words are employed to explain tree-based models,
but the explanation is limited to feature importance and lacks mapping to time
series as saliency maps. Le Nguyen et al., (2019) demonstrate how coefficients
assigned to SAX words in a linear model can be mapped back to the original
time series. However, this explanation is class-wise and insufficient for
interpreting the prediction outcomes of individual time series instances. We
emphasize how borf naturally enables local explanations based on saliency maps
and feature importance.

The most famous pattern-based TSC approaches are based on shapelets. These
methods focus on finding short patterns that define a class independently of
the position of the pattern (Ye and Keogh,, 2011). The discriminative feature
is usually a distance, which can be used to convert the time series dataset
into a tabular one through the so-called shapelet transform (Lines et al.,,
2012). Many approaches in this field use different ways of extracting
shapelets and building the subsequent classifier. For example, in the Random
Shapelet Forest (rsf) algorithm (Karlsson et al.,, 2016), the shapelets are
sampled randomly, building a forest of decision trees with these shapelets as
splitting points. According to (Bagnall et al.,, 2017), the Shapelet Transform
(st) (Hills et al.,, 2014) algorithm is still the most accurate shapelet
approach, with the drawback of the high computational complexity resulting in
rather long training and inference times. To the best of our knowledge, rsf is
the only shapelet-based approach that can be used for TSER.

Interval-based classifiers, such as Time Series Forest (tsf) (Deng et al.,,
2013) and bag-of-features framework (tsbf) (Baydogan et al.,, 2013), divide
time series into random intervals and extract summary statistics like mean,
standard deviation, and slope, training a random forest on the resulting
dataset. The state-of-the-art approach in this category is the Canonical
Interval Forest (cif) (Middlehurst et al.,, 2020), which combines tsf and
Random Interval Spectral Ensemble (rise) (Flynn et al.,, 2019) methods.
Notably, these approaches heavily rely on randomization, causing variations in
the areas of interest used by the classifier across runs. In contrast, borf is
entirely deterministic, generating repeatable representations.

Finally, rocket (Dempster et al.,, 2020) generates many random convolutional
kernels to transform time series into feature vectors, which can be used with
any linear classifier or regressor, often achieving high accuracy and
processing speed even with large and complex datasets. rocket is currently the
best-performing model in both TSC and TSER (Ruiz et al.,, 2021; Tan et al.,,
2021). Inspired by the great performance of such an approach, we also
introduce convolution operators inside borf. As a final note, deep learning
approaches often fall behind in this task (Ruiz et al.,, 2021), given that
they require ad-hoc fine-tuning, often dependent on the specific task and
dataset.

##  3 Background

This section provides all the necessary concepts to understand our proposal.
We begin by defining time series data.

######  Definition 1 (Time Series Data).

A time series dataset,
$\mathcal{{X}}=\\{{X}_{1},\dots,{X}_{n}\\}\in\dot{\mathbb{R}}^{n\times c\times
m}$, is a collection of $n$ time series. A time series, ${X}$, is a collection
of $c$ signals (or channels),
${X}=\\{\textbf{x}_{1},\dots,\textbf{x}_{c}\\}\in\dot{\mathbb{R}}^{c\times
m}$. A signal, x, is a sequence of $m$ real-valued observations sampled at
equal time intervals, $\textbf{x}=[x_{1},\dots,x_{m}]\in\dot{\mathbb{R}}^{m}$.
Associated with each signal, is a timestep vector
$\textbf{t}=[t_{1},\dots,t_{m}]=[1,\dots,m]\in\mathbb{N}^{m}$, which contains
the index of the recorded observation.

When $c=1$, the time series is univariate, else it is multivariate. For
simplicity of notation, we use a unique symbol $m$ to denote the lengths of
the time series. However, signals in each time series can generally have
different lengths. We use the overdot over a set, (e.g., $\dot{\mathbb{R}}$)
to indicate the extension to NaN values, representing a symbol for a missing
observation (e.g., $\mathbb{R}\cup\\{\textit{NaN}\\}$).

Time series datasets can be used in a variety of tasks. This work focuses on
supervised learning, particularly time series extrinsic predictions.

######  Definition 2 (Time Series Extrinsic Prediction).

Given a time series dataset, $\mathcal{{X}}$, Time Series Extrinsic Prediction
is the task of training a model, $f$, to predict an output, $y$, for each
input time series, ${X}$, i.e.,
$f(\mathcal{{X}})=[f({X}_{1}),\dots,f({X}_{n})]=\textbf{y}$.

If the values in y are categorical, i.e., $\textbf{y}\in\mathbb{N}^{n}$, the
task is called Time Series Classification (TSC). If the values are scalars,
i.e., $\textbf{y}\in\mathbb{R}^{n}$, the task is called Time Series Extrinsic
Regression (TSER).

Two of the most common XAI techniques for explaining extrinsic predictions are
feature importance and saliency maps (Bodria et al.,, 2023).

######  Definition 3 (Feature Importance).

Give an instance of a tabular dataset, $\textbf{p}\in\mathbb{R}^{h}$, a
feature importance vector,
$\boldsymbol{\phi}=[\phi_{1},\dots,\phi_{h}]\in\mathbb{R}^{h}$, contains a
score for each value in p.

######  Definition 4 (Saliency Map).

Given a time series, ${X}$, a saliency map,
$\Psi=[\psi_{1,1},\dots,\psi_{c,m}]\in\mathbb{R}^{c,m}$, contains a score,
$\psi_{k,j}$, for every observation, $x_{k,j}$, in ${X}$.

These are similar concepts that differ in the input type. Feature importance
highlights the relevance of features in a tabular dataset, while saliency maps
are adopted in computer vision to highlight the contribution of single
observations toward time series predictions. These explanations are strongly
connected in our proposal, as shown in Section 5.

##  4 Proposed Method

Input: $\mathcal{{X}}$ \- time series dataset, $w$ \- window size, $d$ \-
dilation, $s$ \- stride, $l$ \- word length, $a_{\mu}$ \- mean alphabet size,
$a_{\beta}$ \- slope alphabet size, $\theta$ \- standard deviation threshold

Output: ${P}$ \- Bag-of-Receptive-Fields

1 ${P}\leftarrow\emptyset$;

2for  _${X}_{i}\in\mathcal{{X}}$_ do

3 for  _$\textbf{x}\in{X}$_ do

4 ${R}\leftarrow\textit{{{windowize}}}(\textbf{x},w,d,s)$;

5 for  _$\textbf{r}\in{R}$_ do

6 $\textbf{r}\leftarrow\textit{{{normalize}}}(\textbf{r},\textbf{x},\theta)$;

7 $\tilde{\textbf{r}}\leftarrow\textit{{1d-sax-
na}}(\textbf{r},l,a_{\mu},a_{\beta})$;

8

9 if  _$[*,\tilde{\textbf{r}}]\in P$_ then

10

11 $P[i,\tilde{\textbf{r}}]\leftarrow P[i,\tilde{\textbf{r}}]+1$;

12 else

13 $P[i,\tilde{\textbf{r}}]\leftarrow 1$;

return _${P}$_

Algorithm 1 borf (training)

The Bag-Of-Receptive-Fields (borf) introduces a novel connection between
convolutional operators and pattern extraction by generalizing the concept of
subsequence to that of a receptive field.

######  Definition 5 (Receptive Field).

Given a signal x, a receptive field of length $w\geq 1$ and dilation $d\geq
1$, is an ordered sequence of values,
$\textbf{r}=[x_{j},x_{j+d},\dots,x_{j+d(w-1)}]$.

In convolutional neural networks, the receptive field refers to the spatial
extent of input data that influences the activation of a particular neuron in
the network. Applied to our setting, it can be viewed as a generalized time
series subsequence, i.e., a temporal pattern. borf is an unsupervised
transform that can convert any time series dataset into a tabular
representation in which entries represent the frequency of appearance of
temporal patterns in the time series.

######  Definition 6 (Bag-Of-Receptive-Fields).

Given a time series dataset $\mathcal{{X}}\in\dot{\mathbb{R}}^{n\times c\times
m}$ and a set of $h$ patterns of length $w$, ${R}\in\mathbb{R}^{h\times w}$, a
Bag-Of-Receptive-Fields is a dataset ${P}\in\mathbb{N}^{n\times h}$, where
$p_{i,j}$ is the number of appearances of temporal pattern $j$ in the time
series $i$.

borf can take any real-world univariate and multivariate time series as input,
allowing the presence of irregular signals and missing values. This
flexibility is given by the proposed extension of 1d-sax through complete case
analysis to be able to deal with missing observations, and the sparse borf
representation to manage the high amount of extracted patterns. The borf
representation can be used for TSA, having the advantage of more
comprehensible features, which naturally connect feature and saliency-based
explanations.

Figure 1: A simplified borf schema for a univariate time series ${X}_{1}$. A
receptive field with size $8$ ($w=8$), dilation $2$ ($d=2$), and stride $1$
($s=1$) is first extracted. Then, it is normalized and approximated using a
modified version of 1d-sax (1d-sax-na), to obtain a pattern with a word length
of 2 ($l=2)$, and using 3 symbols for the mean ($a_{\mu}=3$) and 2 for the
slope ($a_{\beta}=2$). Finally, the borf representation is updated by adding
$1$ to the number of times the pattern “$BB,AB$” was observed in ${X}_{1}$.

borf can be schematized into four steps: windowing, normalization,
approximation, and transform, depicted in Figure 1, and presented in Algorithm
1. borf takes as input a time series dataset, $\mathcal{{X}}$, and outputs a
Bag-of-Receptive-Fields ${P}$, using a few hyperparameters that define the
characteristics of the extracted patterns. The first step is to initialize the
empty Bag-Of-Receptive-Fields, ${P}$ (line 1). For building ${P}$, we use the
sparse coordinate format (COO). In this format, the non-zero elements of a
matrix are stored along with their corresponding row and column indices in a
three-column matrix. This format is ideal in our setting where ${P}$ needs to
be built incrementally. For ease of understanding, in Algorithm 1 and Figure
1, we represent this data structure as a single array, ${P}$. However, ${P}$
is never stored in its dense form to ensure space efficiency.

Once ${P}$ is initialized, borf iterates over each signal of each time series
(lines 2-3) to divide it into equally-sized receptive fields (line 4). This is
the windowing step, depicted in Figure 1 (top-left). Each receptive field is
then standardized in the normalization step, using a threshold to identify
constant segments (lines 5-6). The normalization ensures that the pattern’s
scale does not influence its final symbolic representation. In the
approximation step (line 7), the pattern is then discretized using a
generalized version of 1d-sax, (1d-sax-na), and thus converted into a symbolic
word (Figure 1, bottom). Finally, ${P}$ is updated in the transform step
(lines 8-11) by checking if the retrieved word is already in ${P}$ (line 8).
If that is the case, the Bag-of-Receptive-Fields entry is increased by $1$
(Figure 1, top-right, and line 9), else a new entry is created and set to $1$
(lines 10-11).

Algorithm 1 shows the training process for a single configuration of borf. We
iterate this procedure $g$ times, changing hyperparameters to obtain patterns
at multiple resolutions, i.e., receptive fields that capture different input
signal properties. If two runs with different configurations yield the same
pattern, they will be treated as two separate patterns with separate counts.
In simple words, patterns from each configuration are independent, and the
resulting borf representations are stacked horizontally:
${P}=\left[{P}_{1}\,|\,{P}_{2}\,|\,\ldots\,|\,{P}_{g}\right]$ . This allows
for running the algorithm with multiple configurations in a parallel fashion.
In the experimental section, we discuss parameters for the multiple
resolutions and the algorithm’s complexity with empirical evidence of its
scalability, showing its possible usage as an interpretable machine learning
approach. All the steps are thoroughly explained in the following sections.

### Windowing.

This step introduces a strided and dilated windowing approach to extract
receptive fields. Dilation plays an essential role in convolution, offering
control over receptive field size and determining the spacing between
convolutional filter elements, enabling the capture of long-range dependencies
and patterns. By introducing gaps between elements, dilation expands the
receptive field and facilitates the incorporation of distant information. Via
sliding window, we can extract a set of receptive fields, which capture
different aspects of a time series at various levels of resolution. Formally,
we define a windowize function that takes as input a signal x, and, by varying
$j$, extracts a set of receptive fields
${R}={\textit{windowize}}(\textbf{x},w,d,s)$, with $j=1+i\cdot s$ and $0\leq
i\leq\left\lfloor\nicefrac{{(m-w-(d-1)(w-1))}}{{s}}\right\rfloor$. The stride
$s$ determines the displacement of the filter as it moves along the input
data, allowing for downsampling and adjustment of output resolution. The
receptive field definition generalizes the concept of subsequence, i.e., a
subsequence obtained through a standard sliding window is a receptive field
having $d=1$ and $s=1$.

### Normalization.

For each extracted receptive field, $\textbf{r}\in{R}$, borf performs
normalization through a thresholded standardization in the normalize function.
Formally, given a single point $x$ of r:

${\textit{normalize}}(x)=\dfrac{x-\mu_{\textbf{r}}}{\sigma_{\textbf{r}}}\cdot\textbf{1}\left[\dfrac{\sigma_{\textbf{r}}}{\sigma_{\textbf{x}}}\geq\theta\right]$

where $\mu_{\textbf{r}}$ and $\sigma_{\textbf{r}}$ are the average and
standard deviation of r, $\sigma_{\textbf{x}}$ is the standard deviation of
the full signal, x, $0\leq\theta\leq 1$ is a threshold parameter, and 1 is the
indicator function. This threshold avoids blowing up the noise of (almost)
constant, receptive fields. If
${\sigma_{\textbf{r}}}/{\sigma_{\textbf{x}}}<\theta$, we assume the segment to
be constant, and we directly set it to $0$. At this step, we deal with missing
values through complete case analysis, ignoring them in the computation of the
mean and standard deviation. If the receptive field, or even the full signal,
lacks valid values, we leave the NaN vector as is.

### Approximation.

Once the receptive fields are extracted and normalized, they can be converted
to symbolic words through the 1d-sax-na function. We specifically chose 1d-sax
as it is a standard generalization of sax, and for its linear complexity and
easy interpretability. Given r, we first segment the receptive fields in $l$
segments containing $q=w/l$ points. Contrary to the original sax, in the case
where $w$ is not divisible by $l$, we allow segments to have a slightly
different number of observations. The edges indices are therefore computed as
follows:
$\textbf{e}=\left[1,\left\lfloor\nicefrac{{(w-1)}}{{{l}}}\right\rfloor+1,\left\lfloor\nicefrac{{(2w-1)}}{{l}}\right\rfloor+1,\ldots,\left\lfloor\nicefrac{{(lw-1)}}{{l}}\right\rfloor+1\right]\in\mathbb{N}^{l+1}$.
The receptive field is divided then into a set of segments,
$S=\\{\textbf{r}_{e_{1}:e_{2}},\dots,\textbf{r}_{e_{l}:e_{l+1}}\\}\in\dot{\mathbb{R}}^{l\times
q}$, with their respective original time index,
$T=\\{\textbf{t}_{e_{1}:e_{2}},\dots,\textbf{t}_{e_{l}:e_{l+1}}\\}\in{\mathbb{N}}^{l\times
q}$. Given that each segment can possibly contain missing values, we
generalize PLA through complete case analysis. For each of the (almost) equal-
length segments, $\textbf{r}_{e_{i}:e_{i+1}}$ and
$\textbf{t}_{e_{i}:e_{i+1}}$, (r and t for simplicity), we create an index
vector, containing only the indices of valid values, i.e., $\textbf{i}=[i\mid
r_{i}\neq\textit{NaN}]$. We can now perform a modified PLA via least squares.
We compute the sample mean, iterating over the index vector i:

$\mu=\begin{cases}|\textbf{i}|^{-1}\sum_{i\in\textbf{i}}x_{t_{i}}&\quad|\textbf{i}|\geq
1,\\\ \textit{NaN}&\quad\textit{otherwise}.\end{cases}$

Then, the slope can be computed as follows:

$\beta=\begin{cases}\dfrac{\sum_{i\in\textbf{i}}(t_{i}-|\textbf{i}|^{-1}\sum_{i\in\textbf{i}}t_{i})(x_{t_{i}}-\mu)}{\sum_{i\in\textbf{i}}(t_{i}-|\textbf{i}|^{-1}\sum_{i\in\textbf{i}}t_{i})^{2}}&\quad|\textbf{i}|\geq
2,\\\ 0&\quad|\textbf{i}|=1,\\\
\textit{NaN}&\quad\textit{otherwise}.\end{cases}$

In simple terms, the slope can be computed via PLA only if at least two valid
observations exist in the segment. If only one valid value exists, the slope
is set to zero. Finally, the slope and average are set to NaN if no valid
value exists. At this point, we have two vectors, one containing the averages
$\boldsymbol{\mu}=[{\mu}_{1},\dots,{\mu}_{l}]\in\dot{\mathbb{R}}^{l}$, and one
containing the slopes,
$\boldsymbol{\beta}=[\beta_{1},\dots,\beta_{l}]\in\dot{\mathbb{R}}^{l}$, of
each segment of the receptive field. Values in each vector are then quantized
using two separate alphabets, $\mathbb{{A}}_{\mu}$ and $\mathbb{{A}}_{\beta}$,
containing $a_{\mu}\geq 2$ and $a_{\beta}\geq 1$ symbols, respectively. In
practice, depending on their value, elements in the vectors are assigned a
symbol. Specifically, the average values are quantized using $a_{\mu}$
symbols, obtained through the quantiles of the Gaussian distribution
$\mathcal{N}(0,1)$, while the slope values are quantized on $a_{\beta}$ levels
according to the quantiles of the Gaussian distribution
$\mathcal{N}(0,\nicefrac{{0.03}}{{(w\cdot d)}})$, where $0.03$ is a scaling
factor proposed by Malinowski et al., (2013) in the original 1d-sax paper, and
$w\cdot d$ is the effective time-span adjusted by the dilation of the
receptive field. Once quantized, the symbols are concatenated to form the
final approximated receptive field $\textit{{1d-sax-
na}}(\textbf{r},l,a_{\mu},a_{\beta})=[\textit{prefix},\tilde{\mu}_{1}\tilde{\beta}_{1},\dots,\tilde{\mu}_{l}\tilde{\beta}_{l}]=\tilde{\textbf{r}}\in\mathbb{A}^{l+1}$
with $|\mathbb{A}|=a_{\mu}\cdot a_{\beta}+1$, where the $+1$ counts for the
symbol added for denoting $NaN$ average and slope values. The prefix contains
the configuration and signal id, i.e., words extracted from different
configurations and different time series signals are independent of each
other. If $a_{\beta}=1$, 1d-sax-na reduces to standard sax.

### Transform.

Given a newly created word, $\tilde{\textbf{r}}$, we need to store it. First,
we check if $\tilde{\textbf{r}}$ already exists in ${P}$. If it does, i.e.,
$[*,\tilde{\textbf{r}}]\in P$, we increase the count in the entry
${P}[i,\tilde{\textbf{r}}]$ by $1$, where $i$ is the index of the time series
we are processing. If, instead, the word was never observed, the update
process differs at training and inference time. At training time, we
initialize a new entry in ${P}[i,\tilde{\textbf{r}}]$ by setting it to $1$. At
inference time, we can not consider it, given that it was never observed
during training. In other words, at inference time, the algorithm only
performs updates of ${P}$ for words found at training time, i.e., no new word
is added in this phase. At the end of the algorithm, ${P}$ contains the
overall count of appearances of each word in each time series. Thus, any model
trained on this data will base its predictions on easily understandable
features.

##  5 Experiments

### Datasets.

We test borf on TSC and TSER on the full main repositories proposed in the
literature. For TSC, we use the 158 datasets from the UCR and UEA TSC
repositories (Bagnall et al.,, 2018; Dau et al.,, 2019), which contain 128
univariate and 30 multivariate datasets. For TSER, we use the UEA and UCR Time
Series Extrinsic Regression Repository (Tan et al.,, 2020), which contains 19
datasets, 4 univariate, and 15 multivariate. We adopt the default training and
test splits. Some datasets contain time series with a variable number of
observations, and some contain missing values. For methods that do not support
this kind of data, we apply missing value imputation by linear interpolation
and last value padding. Further, we concatenate all signals in a single axis
for algorithms that do not support multivariate time series.

### BORF Hyperparameters.

Following the current trend in TSC and TSER (Tan et al.,, 2021; Ruiz et al.,,
2021), instead of fine-tuning borf to maximize performance for specific
applications, we propose a model that achieves good results in a wide range of
problems. Since several hyperparameters exist, we provide a heuristic that
computes a vector of window sizes, dilations, and word lengths as follows:
$\textbf{w}=[2^{2},2^{3},\dots,2^{\lfloor\log_{2}m\rfloor}]$,
$\textbf{d}=[2^{0},2^{1},\dots,2^{\lfloor\log_{2}(\log_{2}m)\rfloor}]$,
$\textbf{l}=[2^{0},2^{1},2^{2},2^{3}]$. Once these vectors are computed, each
valid combination of elements is used as a configuration of borf, resulting in
a multi-resolution representation of the time series. The stride is always set
to $1$, $s=1$. We only tune the following $3$ parameters once for each
task111Code available at https://github.com/fspinna/borf. For TSC we set
$a_{\mu}=2$, $a_{\beta}=3$, and $\theta=0.15$. For TSER, $a_{\mu}=3$,
$a_{\beta}=1$ and $\theta=0.05$.

After the transformation, many models can be used to perform supervised tasks.
We focused on fast linear or tree-based approaches. From preliminary
benchmarks on a subset of datasets, we found that linear models tend to
perform best in TSC, while tree-based models are better suited for TSER.
Linear models usually benefit from more normally-shaped features. For this
reason, we apply an arcsinh function element-wise after the borf transform and
before a linear-svc classifier with default parameters. For TSER, we use a
LightGBM regressor with default parameters.

### Competitors.

We compare borf against several state-of-the-art approaches presented in
Section 2. For TSC, as a common baseline, we use KNN with DTW as distance and
a Sakoe-Chiba band of $0.1$ (knn-dtw). For competitor dictionary-based
approaches, we test bop, sax-vsm, mr-sqm, a boss ensemble, and muse. For
interval-base approaches, we use Time Series Bag of Features (tsbf), Time
Series Forest (tsf), and Canonical Interval Forest (cif). For shapelet-based,
we benchmark Random Shapelet Forest (rsf) and the standard Shapelet Transform
(st). Finally, we also compare against rocket. For TSER, there are fewer
choices, and we benchmark all of the previous that can be directly applied to
TSER, i.e., knn-dtw, bop, tsf, sf, and rocket. We also directly test a
LightGBM regressor (lgbm) on the raw time series. All models are trained with
the default library hyperparameters or values proposed in the respective
papers. Each model is allowed approximately one week ($10,000$ minutes) for
training and inference on each dataset, and is allocated $256GB$ of memory. A
run is considered failed if it exceeds the maximum time or crashes due to out-
of-memory errors. Each benchmark is performed $3$ times with different seeds
for models with a random component, and the average performance is taken.

Figure 2: CD diagram for TSC in terms of bACC (left), and for TSER in terms of
MAPE (right). Best models to the right.

### Results.

The comparison of the ranks of all methods against each other considering all
datasets is presented in Figure 2 with Critical Difference (CD) diagrams
(Demsar,, 2006). Two methods are tied if the null hypothesis that their
performance is the same cannot be rejected using the Nemenyi test at
$\alpha{=}0.05$. The performance metric used for TSC is Balanced Accuracy
(bACC) (Brodersen et al.,, 2010), which accounts for both false positives and
false negatives when dealing with imbalanced data. For TSER, we adopt Mean
Absolute Percentage Error (MAPE). In both plots, best values have lower ranks
and appear on the right-hand side. For TSC, borf scores second place and is
statistically tied to rocket. Notably, borf outperforms significantly SAX-
based methods such as mr-sqm and sax-vsm, and is the best-performing
dictionary-based approach. For TSER, borf is the overall best approach,
outperforming even state-of-the-art black-box competitors. It is interesting
to see that lgbm, which does not consider the sequentiality of data,
outperforms every other method, except borf and tsf. A similar observation was
made in (Tan et al.,, 2021), highlighting the need for more research for this
task. Another issue is the lack of big repositories for TSER, which produces a
CD plot that cannot definitively separate the top-performing approaches.

Figure 3: Mean performance and training time between borf and competitors for
TSC (top) and TSER (bottom). Best models to the right.

For this reason, we also present a performance versus time comparison in
Figure 3 (less competitive approaches are shaded). The x-axis shows the mean
performance of the models, in terms of bACC for TSC, and R2 for TSER.
Therefore, the best-performing models are on the right. On the y-axis is the
training time in minutes. Thus, faster approaches are at the bottom. Regarding
TSC, the overall best approach is rocket, while cif and st are close to borf
but are also slower. We highlight that while rocket performs better in TSC, it
is also a black-box, completely opaque from a human standpoint. In TSER, borf
is the best-performing method. lgbm, tsf, and rocket are the only comparable
ones, with a noticeable performance gap. In that group lgbm performs
comparably to tsf and rocket, while being notably faster. In general, while
borf exhibits average training runtimes, it is also the only approach that
offers competitive performance while being deterministic, which is a crucial
feature for XAI purposes. Additionally, it is worth noting that because not
all approaches support parallel computation, runtimes are recorded using a
single thread to ensure better comparability. These results represent worst-
case scenarios, and many of the benchmarked approaches, including borf, can be
significantly accelerated through multiprocessing.

### Complexity and Scalability.

The complexity of borf is dominated by that of 1d-sax-na, which for a sequence
of $w$ points is linear, i.e., $O(w)$ (Malinowski et al.,, 2013)222Notice that
lines 8–11 of Algorithm 1, potentially expensive if implemented with naive COO
matrix operations, are performed in constant time through hashmaps.. For a
single signal, 1d-sax-na has to be performed one time for each receptive
field, i.e., $\left\lfloor\nicefrac{{(m-w-(d-1)(w-1))}}{{s}}\right\rfloor+1$
times. The worst case scenario is when $d=1$, $s=1$, and
$w=\nicefrac{{(m+1)}}{{2}}$. Applied to the full dataset and repeated for $g$
configurations, the time complexity of borf is $O(g\cdot n\cdot c\cdot
m^{2})$.

Regarding space complexity, there are two aspects to consider. First, the
number of possible extracted words for a given configuration depends on the
word length, $l$, and alphabet size, $a$, and is exponential, $O(a^{l})$.
However, there cannot be more words than receptive fields. Thus, in the worst
case scenario, when there are no duplicate words, the complexity is $O(g\cdot
c\cdot\min(a^{l},n\cdot m^{2}))$. Given the proposed heuristic, the effective
number of receptive fields is at least one order of magnitude smaller than
this bound, as shown in Figure 4 (left). For each dataset, Figure 4 (left)
shows the number of time series points (x-axis) versus the number of extracted
patterns (y-axis, blue) and the theoretical maximum of extractable patterns
(y-axis, orange). In particular, the empirical space complexity is linear
w.r.t. the number of points and very close to $n\cdot c\cdot m$. Further, as
shown in Figure 4 (right), the density (x-axis) of the borf representation
naturally decreases as the number of features increases (y-axis), with the
density, $\rho$, being defined as the ratio between the number of non-zero
entries w.r.t. the number of entries in ${P}$, $\rho=\textit{nnz}({P})/|{P}|$.
Models developed for this kind of representation tend to be more efficient as
the number of features increases, exploiting the increasingly lower proportion
of non-zero elements in the input matrix (Pedregosa et al.,, 2011).

Figure 4: Number of points against extracted patterns (left). Density against
number of patterns (right).

### Interpretability.

Here we show how a feature-based explanation can be partially mapped to a
saliency map in an borf representation. Given a time series, ${X}$, its borf
conversion, p, and a feature importance vector, $\boldsymbol{\phi}$, obtained
through some local explainer, such as SHAP (Lundberg and Lee,, 2017), we
propose a way to build a saliency map, $\Psi$. A single feature importance,
$\phi_{i}$, for a single symbolic pattern, $p_{i}$, can be mapped back to the
original timesteps only if $p_{i}$ is contained in the time series, i.e.,
$p_{i}>0$. If that is the case, we can retrieve the multiset
$T_{i}=\\{(k_{1},j_{1}),(k_{2},j_{2}),\dots\\}$, containing all the alignment
indices, possibly duplicate, for $p_{i}$. The saliency matrix is then defined
as $\Psi^{*}=[\psi_{1,1},\dots,\psi_{c,m}]\in\mathbb{R}^{c\times m}$, where
$\psi_{k,j}=\sum_{i=1}^{h}\phi_{i}\cdot|\\{(k,j)\in T_{i}\\}|.$ To obtain the
final saliency map, $\Psi^{*}$ is multiplied by a scaling factor:

$\Psi=\Psi^{*}\cdot\frac{\sum_{i=1}^{h}{\phi_{i}\cdot\textbf{1}[p_{i}>0]}}{\sum_{k=1}^{c}\sum_{j=1}^{m}{\psi_{k,j}}}$

In practice, the sum of the saliency values is scaled to equal the sum of the
feature importances of contained patterns. Hence, the final explanation
comprises a saliency map explaining the contained patterns, and a feature
importance vector, explaining the contribution of not-contained patterns.

In Figure 5, we show an example of such an explanation on one of the most
challenging and novel tasks from an interpretability standpoint, i.e.,
multivariate TSER. The time series is from the HouseholdPowerConsumption1
dataset, containing $5$ signals of the minutely voltage, current and sub-meter
energy of a house near Paris, France (Hebrail and Berard,, 2012). The task is
to predict the global active power in kiloWatt (kW), and the prediction of
borf ($1460.4kW$) is extremely close to the real one ($1460.9kW$). The feature
importance vector of p is obtained through SHAP. In Figure 5 (left), the
multivariate time series is shown, and each point is colored based on its
importance. Positive SHAP values (red) indicate that the observation “pushes”
the prediction toward higher active power, while negative values indicate
lower power consumption. Values close to $0$ are not relevant to the
prediction. On the right, we show each signal’s most important not-contained
patterns, i.e., $p_{i}=0$, colored based on their contribution toward the
prediction. These shapes are obtained by averaging all the receptive fields
corresponding to a specific SAX word, for each instance of the training set.
At first glance, it is easy to see that the most important signals for the
regression are sub-meter energy 2 and 3 ($SubMt_{2}$ and $SubMt_{3}$),
followed by the current ($Amp$). Voltage ($Volt$) and $SubMt_{1}$ are almost
all gray, so they do not contribute significantly toward the prediction. From
a layman’s point of view, it seems that rapid increases in sub-meter energy
($SubMt_{2}$, $SubMt_{3}$) are indicators of an increase in consumption.
Further, the most important pattern for $SubMt_{2}$ is “U-shaped”, and its
absence strongly indicates increased active power. $SubMt_{3}$ mainly pushes
the prediction toward lower power consumption, and the not-contained patterns
are harder to interpret, having a higher frequency. In this case, the
explanation is more technical and requires a domain expert.

Figure 5: A saliency and feature-based explanation on a multivariate time
series of the HouseholdPowerConsumption1 dataset. (left) original time series;
(right) most important not-contained patterns. The more intense the color, the
more important the feature.

##  6 Conclusion

We have presented borf, a Bag-Of-Receptive-Fields deterministic transform
that, by incorporating time series convolution operators and 1d-sax, addresses
the limitations of existing TSA techniques, offering improved performance,
interpretability, and applicability in real-world scenarios. We evaluated borf
on TSC and TSER tasks, showcasing its superior performance in TSER and
competitive results in TSC. Additionally, we outlined how borf can be used to
achieve saliency and feature-based explanations. A limitation of borf for very
big datasets is the quadratic time complexity w.r.t. the number of
observations. For future research directions, we plan on exploring the stride
parameter to achieve lower running times without sacrificing performance.
Further, we plan on extending univariate receptive fields to multivariate ones
to generate symbolic words that can span more than one signal. Regarding data,
the next natural obstacle is being able to deal with unequally sampled time
series in which the time dimension is not constant. Finally, we plan on
extending borf to achieve interpretable time series forecasting and test it on
unsupervised tasks like clustering and anomaly detection.

## References

  * Bagnall et al., (2018) Bagnall, A., Dau, H. A., Lines, J., Flynn, M., Large, J., Bostrom, A., Southam, P., and Keogh, E. (2018).  The uea multivariate time series classification archive, 2018.  arXiv preprint arXiv:1811.00075.
  * Bagnall et al., (2017) Bagnall, A., Lines, J., Bostrom, A., Large, J., and Keogh, E. (2017).  The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances.  Data mining and knowledge discovery, 31:606–660.
  * Baydogan et al., (2013) Baydogan, M. G., Runger, G., and Tuv, E. (2013).  A bag-of-features framework to classify time series.  IEEE transactions on pattern analysis and machine intelligence, 35(11):2796–2802.
  * Bodria et al., (2023) Bodria, F., Giannotti, F., Guidotti, R., Naretto, F., Pedreschi, D., and Rinzivillo, S. (2023).  Benchmarking and survey of explanation methods for black box models.  Data Mining and Knowledge Discovery, pages 1–60.
  * Brodersen et al., (2010) Brodersen, K. H., Ong, C. S., Stephan, K. E., and Buhmann, J. M. (2010).  The balanced accuracy and its posterior distribution.  In 2010 20th international conference on pattern recognition, pages 3121–3124. IEEE.
  * Dau et al., (2019) Dau, H. A., Bagnall, A., Kamgar, K., Yeh, C.-C. M., Zhu, Y., Gharghabi, S., Ratanamahatana, C. A., and Keogh, E. (2019).  The ucr time series archive.  IEEE/CAA Journal of Automatica Sinica, 6(6):1293–1305.
  * Dempster et al., (2020) Dempster, A., Petitjean, F., and Webb, G. I. (2020).  Rocket: exceptionally fast and accurate time series classification using random convolutional kernels.  Data Mining and Knowledge Discovery, 34(5):1454–1495.
  * Demsar, (2006) Demsar, J. (2006).  Statistical comparisons of classifiers over multiple data sets.  J. Mach. Learn. Res., 7:1–30.
  * Deng et al., (2013) Deng, H., Runger, G., Tuv, E., and Vladimir, M. (2013).  A time series forest for classification and feature extraction.  Information Sciences, 239:142–153.
  * Flynn et al., (2019) Flynn, M., Large, J., and Bagnall, T. (2019).  The contract random interval spectral ensemble (c-rise): the effect of contracting a classifier on accuracy.  In Hybrid Artificial Intelligent Systems: 14th International Conference, HAIS 2019, León, Spain, September 4–6, 2019, Proceedings 14, pages 381–392. Springer.
  * Hebrail and Berard, (2012) Hebrail, G. and Berard, A. (2012).  Individual household electric power consumption.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C58K54.
  * Hills et al., (2014) Hills, J., Lines, J., Baranauskas, E., Mapp, J., and Bagnall, A. (2014).  Classification of time series by shapelet transformation.  Data mining and knowledge discovery, 28:851–881.
  * Hung and Anh, (2007) Hung, N. Q. V. and Anh, D. T. (2007).  Combining sax and piecewise linear approximation to improve similarity search on financial time series.  In 2007 International Symposium on Information Technology Convergence (ISITC 2007), pages 58–62. IEEE.
  * Karlsson et al., (2016) Karlsson, I., Papapetrou, P., and Boström, H. (2016).  Generalized random shapelet forests.  Data mining and knowledge discovery, 30:1053–1085.
  * Keogh et al., (2001) Keogh, E., Chakrabarti, K., Pazzani, M., and Mehrotra, S. (2001).  Dimensionality reduction for fast similarity search in large time series databases.  Knowledge and information Systems, 3:263–286.
  * Kreindler and Lumsden, (2006) Kreindler, D. M. and Lumsden, C. J. (2006).  The effects of the irregular sample and missing data in time series analysis.  Nonlinear dynamics, psychology, and life sciences, 10(2):187–214.
  * Le Nguyen et al., (2019) Le Nguyen, T., Gsponer, S., Ilie, I., O’reilly, M., and Ifrim, G. (2019).  Interpretable time series classification using linear models and multi-resolution multi-domain symbolic representations.  Data mining and knowledge discovery, 33:1183–1222.
  * Lin et al., (2007) Lin, J., Keogh, E., Wei, L., and Lonardi, S. (2007).  Experiencing sax: a novel symbolic representation of time series.  Data Mining and knowledge discovery, 15:107–144.
  * Lin et al., (2012) Lin, J., Khade, R., and Li, Y. (2012).  Rotation-invariant similarity in time series using bag-of-patterns representation.  Journal of Intelligent Information Systems, 39:287–315.
  * Lines et al., (2012) Lines, J., Davis, L. M., Hills, J., and Bagnall, A. (2012).  A shapelet transform for time series classification.  In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 289–297.
  * Lundberg and Lee, (2017) Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  Advances in neural information processing systems, 30.
  * Malinowski et al., (2013) Malinowski, S., Guyet, T., Quiniou, R., and Tavenard, R. (2013).  1d-sax: A novel symbolic representation for time series.  In Advances in Intelligent Data Analysis XII: 12th International Symposium, IDA 2013, London, UK, October 17-19, 2013. Proceedings 12, pages 273–284. Springer.
  * Middlehurst et al., (2020) Middlehurst, M., Large, J., and Bagnall, A. (2020).  The canonical interval forest (cif) classifier for time series classification.  In 2020 IEEE international conference on big data (big data), pages 188–195. IEEE.
  * Nguyen and Ifrim, (2022) Nguyen, T. L. and Ifrim, G. (2022).  Fast time series classification with random symbolic subsequences.  In International Workshop on Advanced Analytics and Learning on Temporal Data, pages 50–65. Springer.
  * Pedregosa et al., (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., and Duchesnay, E. (2011).  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830.
  * Ruiz et al., (2021) Ruiz, A. P., Flynn, M., Large, J., Middlehurst, M., and Bagnall, A. (2021).  The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances.  Data Mining and Knowledge Discovery, 35(2):401–449.
  * Schäfer, (2015) Schäfer, P. (2015).  The boss is concerned with time series classification in the presence of noise.  Data Mining and Knowledge Discovery, 29:1505–1530.
  * Schäfer and Högqvist, (2012) Schäfer, P. and Högqvist, M. (2012).  Sfa: a symbolic fourier approximation and index for similarity search in high dimensional datasets.  In Proceedings of the 15th international conference on extending database technology, pages 516–527.
  * Schäfer and Leser, (2017) Schäfer, P. and Leser, U. (2017).  Multivariate time series classification with weasel+ muse.  arXiv preprint arXiv:1711.11343.
  * Senin and Malinchik, (2013) Senin, P. and Malinchik, S. (2013).  Sax-vsm: Interpretable time series classification using sax and vector space model.  In 2013 IEEE 13th international conference on data mining, pages 1175–1180. IEEE.
  * Shumway et al., (2000) Shumway, R. H., Stoffer, D. S., and Stoffer, D. S. (2000).  Time series analysis and its applications, volume 3.  Springer.
  * Spinnato et al., (2022) Spinnato, F., Guidotti, R., Nanni, M., Maccagnola, D., Paciello, G., and Farina, A. B. (2022).  Explaining crash predictions on multivariate time series data.  In International Conference on Discovery Science, pages 556–566. Springer.
  * Tan et al., (2020) Tan, C. W., Bergmeir, C., Petitjean, F., and Webb, G. I. (2020).  Monash university, uea, ucr time series extrinsic regression archive.  arXiv preprint arXiv:2006.10996.
  * Tan et al., (2021) Tan, C. W., Bergmeir, C., Petitjean, F., and Webb, G. I. (2021).  Time series extrinsic regression: Predicting numeric values from time series data.  Data Mining and Knowledge Discovery, 35:1032–1060.
  * Theissler et al., (2022) Theissler, A., Spinnato, F., Schlegel, U., and Guidotti, R. (2022).  Explainable ai for time series classification: A review, taxonomy and research directions.  IEEE Access.
  * Ye and Keogh, (2011) Ye, L. and Keogh, E. (2011).  Time series shapelets: a novel technique that allows accurate, interpretable and fast classification.  Data mining and knowledge discovery, 22:149–182.


pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# Scale-invariant dynamics in the Solar System

Indranil Banik1 and Pavel Kroupa1,2
1Helmholtz-Institut für Strahlen und Kernphysik (HISKP), University of Bonn,
Nussallee 14$-$16, D-53115 Bonn, Germany
2Charles University, Astronomical Institute, Faculty of Mathematics and
Physics, V Holešovičkách 2, CZ-180 00 Praha 8, Czech Republic Email:
<EMAIL_ADDRESS>(Indranil Banik)
<EMAIL_ADDRESS>(Pavel Kroupa)

###### Abstract

The covariant scale-invariant dynamics (SID) theory has recently been proposed
as a possible explanation for the observed dynamical discrepancies in galaxies
(Maeder & Gueorguiev 2020). SID implies that these discrepancies $-$ commonly
attributed to dark matter $-$ arise instead from a non-standard velocity-
dependent force that causes two-body near-Keplerian orbits to expand. We show
that the predicted expansion of the Earth-Moon orbit is incompatible with
lunar laser ranging data at $>200\sigma$. Moreover, SID predicts that the
gravitating mass of any object was much smaller in the past. If true, a low-
mass red giant star must be significantly older than in standard theory. This
would make it much older than the conventional age of the Universe, which
however is expected to be similarly old in SID. Moreover, it is not completely
clear whether SID truly contains new physics beyond General Relativity, with
several previous works arguing that the extra degree of freedom is purely
mathematical. We conclude that the SID model is falsified at high significance
by observations across a range of scales, even if it is theoretically well
formulated.

###### keywords:

gravitation – dark matter – ephemerides – Astrometry and celestial mechanics –
Moon – space vehicles

††pubyear: 2020††pagerange: Scale-invariant dynamics in the Solar
System–References

##  1 Introduction

The Large Synoptic Survey Telescope (Ivezić et al., 2019) has recently been
renamed the Vera C. Rubin Observatory in honour of her observational work
establishing what is still one of the great mysteries of astronomy (Rubin &
Ford, 1970). Their work revealed very large dynamical discrepancies between
the actual rotation curves of galaxies and the predictions of Newtonian
gravity applied to their luminous matter distributions, as also shown by
several other authors (e.g. Babcock, 1939; Rogstad & Shostak, 1972). Such
acceleration discrepancies are usually attributed to dark matter haloes
surrounding each galaxy (Ostriker & Peebles, 1973). The dark matter would have
to consist of hypothetical particles not in the well-tested standard model of
particle physics. This interpretation is challenged by continued null
detection of any dark matter particles in sensitive searches, for instance in
11 years of Fermi data on dwarf spheroidal satellites of the Milky Way (Hoof
et al., 2020) and in sensitive terrestrial experiments (Liu et al., 2017).

Regardless of the hypothetical particle physics that might give rise to dark
matter, assuming it holds galaxies together leads to inconsistencies on a
variety of scales ranging from hundreds to millions of parsecs (Kroupa, 2012,
2015). The longest studied problem is that the galactic acceleration
discrepancies follow some remarkable regularities (Famaey & McGaugh, 2012)
that can be summarized as a unique relation between the acceleration $g$
inferred from the rotation curve and the Newtonian gravitational field
$g_{{}_{N}}$ generated by the baryonic distribution (McGaugh et al., 2016).
Such a radial acceleration relation (RAR) was predicted several decades
earlier using Modified Newtonian Dynamics (MOND, Milgrom, 1983). In this
model, the dynamical effects usually attributed to dark matter are instead
provided by an acceleration dependence of the gravity law $-$ the
gravitational field strength $g$ at distance $r$ from an isolated point mass
$M$ transitions from the Newtonian $GM/r^{2}$ law to

$\displaystyle g~{}=~{}\frac{\sqrt{GMa_{{}_{0}}}}{r}~{}~{}~{}\text{for
}~{}g,g_{{}_{N}}\ll a_{{}_{0}}\,.$ (1)

MOND (or Milgromian dynamics) introduces $a_{{}_{0}}$ as a fundamental
acceleration scale of nature below which the deviation from Newtonian dynamics
becomes significant. Empirically, $a_{{}_{0}}=1.2\times{10}^{-10}$ m/s2 to
match galaxy rotation curves (Begeman et al., 1991; McGaugh, 2011).

Remarkably, this is approximately where the classical energy density of a
gravitational field (equation 9 of Peters, 1981) becomes comparable to the
dark energy density $u_{{}_{\Lambda}}\equiv\rho_{{}_{\Lambda}}c^{2}$ that
conventionally explains the accelerating expansion of the Universe (Efstathiou
et al., 1990; Ostriker & Steinhardt, 1995; Riess et al., 1998; Perlmutter et
al., 1999).

$\displaystyle\frac{g^{2}}{8\mathrm{\pi}G}~{}<~{}u_{{}_{\Lambda}}~{}~{}\Leftrightarrow~{}~{}g~{}\la~{}2\mathrm{\pi}a_{{}_{0}}\,.$
(2)

In regions of space with such a weak gravitational field, the dominant
contribution to the energy density is $u_{{}_{\Lambda}}$, which is often
identified with the quantum-mechanical zero point energy density of the
vacuum. If this is correct, poorly understood quantum gravity effects could
well be rather important to such regions $-$ but these are totally neglected
in General Relativity. MOND could be an empirical way to include such effects,
as suggested by the coincidence of scales in Equation 2 (e.g. Milgrom, 1999;
Pazy, 2013; Verlinde, 2016; Smolin, 2017).

Regardless of its underlying microphysical explanation, MOND can accurately
match the rotation curves of a wide variety of spiral galaxies across a vast
range in mass, surface brightness, and gas fraction using only the
distribution of luminous matter (figure 5 of Lelli et al., 2017). Fits to
individual rotation curves show that intrinsic scatter about its predictions
must be ${<13\%}$ and is consistent with 0 (Li et al., 2018). A few discrepant
galaxies were claimed by Rodrigues et al. (2018), but it was later shown that
the actual discrepancies are either very mild or arise when the distance is
particularly uncertain and could plausibly be outside the range they allow
(Kroupa et al., 2018a). MOND can also explain the rotation curves of
elliptical galaxies, at least when these can be measured in a sub-dominant
rotating disc of neutral hydrogen (figure 8 of Lelli et al., 2017).

The successes of MOND extend beyond near-circular orbital motion in the non-
relativistic regime. A relativistic version of MOND has recently been
developed in which gravitational waves propagate at the speed of light
(Skordis & Złośnik, 2019), consistent with the near-simultaneous detection of
the gravitational wave event GW170817 and its electromagnetic counterpart
(Virgo & LIGO Collaborations, 2017). Galactic globular clusters can be used to
test MOND in dispersion-dominated systems (Baumgardt et al., 2005; Haghi et
al., 2011). In this regard, NGC 2419 was claimed to be problematic for MOND
(Ibata et al., 2011a, b). It was later shown that the disagreement could be
caused by plausible systematic uncertainties like a mild departure from
spherical symmetry, rotation within the sky plane, or a radially-dependent
polytropic index (Sanders, 2012a, b). In MOND, we also expect NGC 2419 to be
affected by the ‘external field effect’ (EFE) whereby the internal gravity
binding a system is weakened when it is embedded in an external gravitational
field, a consequence of the non-linear MOND equations (Milgrom, 1986).
However, including the Galactic external field did not much improve the
agreement with observations (Derakhshani & Haghi, 2014). MOND can also explain
the velocity dispersion of the ultra-diffuse galaxy Dragonfly 2 (DF2, Kroupa
et al., 2018b), DF4 (Haghi et al., 2019a), and DF44 (Bílek et al., 2019), with
the latter galaxy having a measured dispersion profile that is consistent with
MOND at ${2.40\sigma}$ (Haghi et al., 2019b). The calculations for DF2 and DF4
are significantly affected by the EFE. DF44 is more isolated than the almost
gas-free DF2 (Chowdhury, 2019; Sardone et al., 2019) and indeed has a higher
internal velocity dispersion despite a similar baryonic distribution. In a
MOND context, the EFE was recently confirmed at high significance based on the
relative velocities of wide binary stars in the Solar neighbourhood $-$ these
are inconsistent with MOND without the Galactic EFE (Pittordis & Sutherland,
2019). The more relevant case with the EFE is quite similar to the Newtonian
case, though a more careful analysis could distinguish them in the near future
(Banik & Zhao, 2018; Banik, 2019). Recent progress with hydrodynamic MOND
simulations indicates that it can naturally form exponential disc galaxies out
of a collapsing gas cloud (Wittenburg et al., 2020) and produce realistic
morphologies for the interacting Antennae galaxies (Renaud et al., 2016).

Most of these observations only became apparent long after the MOND field
equation was first published (Bekenstein & Milgrom, 1984), making these
achievements successful a priori predictions. It is difficult to explain the
success of these predictions in a conventional gravity context, even with the
observational facts in hand (Desmond, 2017a, b; Ghari et al., 2019). In
particular, section 4.2 of the latter work showed that it is still rather
difficult to get a tight RAR despite a diversity of rotation curve shapes at
fixed peak velocity (Oman et al., 2015).

Despite its successes, the theoretical underpinnings of MOND remain unclear.
An important clue may be that in the deep-MOND limit, the dynamics of
gravitational systems become scale invariant $-$ scaling the spacetime co-
ordinates by some factor $\lambda$ yields another valid solution to the
equations of motion (Milgrom, 2009). In principle, $\lambda$ could be time-
dependent. This is at the basis of the fully relativistic scale invariant
dynamics theory (SID, Maeder, 1978; Maeder & Bouvier, 1979). SID applies the
Weyl Integrable Geometry (Canuto et al., 1977) and reproduces some of the MOND
phenomenology, in particular a tight RAR (Maeder & Gueorguiev, 2020). There
are some differences in its predicted form when $g_{{}_{N}}\la
0.01a_{{}_{0}}-$ while MOND predicts that $g\to\sqrt{g_{{}_{N}}a_{{}_{0}}}$,
SID predicts that $g$ flattens out at $\approx 0.1a_{{}_{0}}$ (see their
figure 1). This leads to potentially testable consequences in low surface
brightness dwarf galaxies.

In this contribution, we avoid discussing possible theoretical issues with
SID, which may in fact be identical to classical General Relativity as the
extra degree of freedom may be a mathematical artefact (Tsamis & Woodard,
1986; Jackiw & Pi, 2015). We assume that there is genuinely new physics in SID
and focus on its high-acceleration limit, where SID has some unusual
consequences.

After introducing the SID model and its context (Section 1), we discuss how it
can be constrained by Solar System ephemerides (Section 2). We also consider
what SID implies for stellar evolution (Section 3). Our conclusions are given
in Section 4.

##  2 Solar System ephemerides

SID has some success reproducing the RAR (section 5 of Maeder & Gueorguiev,
2020). This comes about because of a non-standard velocity-dependent force
which causes a slow outwards expansion of a two-body near-Keplerian orbit (see
their equation 26). This has always been an important part of the SID theory
(Maeder, 1978). As shown in equation 68 of Maeder & Bouvier (1979), the extra
term causes a near-circular orbit of radius $r$ to expand as

$\displaystyle\dot{r}_{{}_{SID}}~{}=~{}\frac{r}{t}\,,$ (3)

where $t$ is the time since the Big Bang and $\dot{q}\equiv\frac{dq}{dt}$ for
any quantity $q$. The SID subscript indicates a theoretical expectation. In
this contribution, it will be important to measure time from the Big Bang
since the dynamical equations of SID are not invariant under a translation of
the time variable.

Since the expansion rate of the SID universe is similar to the Lambda-Cold
Dark Matter ($\Lambda$CDM) standard cosmological paradigm (section 4.1 of
Maeder & Gueorguiev, 2020), we assume that currently ${t=13.8}$ Gyr (Planck
Collaboration XIII, 2016). Thus, the 9.58 AU orbit of Saturn should expand at
$\dot{r}_{{}_{SID}}=104$ m/yr. Over a decade, the expansion should easily be
detectable in Cassini radio tracking data given that its accuracy is $\approx
32$ m (Viswanathan et al., 2017). However, they found no significant deviation
of Saturn from its conventionally calculated trajectory. It is difficult to
see how a constant outwards expansion could be masked by e.g. changing the
masses of other planets.

###  2.1 Lunar laser ranging

Applying Equation 3 to the Earth-Moon distance of ${r=3.84\times 10^{8}}$ m,
we get that $\dot{r}_{{}_{SID}}=28$ mm/yr. Maeder & Bouvier (1979) predicted
that this effect should be observable in “sufficiently accurate observations
with laser reflectors” (see their section 7). Nowadays, the Earth-Moon
distance is constrained to mm accuracy by lunar laser ranging (LLR, Adelberger
et al., 2017). The Moon is indeed receding from the Earth, but at a rate of
$\dot{r}_{obs}={38.05\pm 0.04}$ mm/yr (section 3c of Folkner et al., 2014).
This is caused by the tides it raises in Earth’s oceans. These tidal bulges
are carried ahead of the sub-lunar point because Earth rotates much faster
than the Moon orbits it. Folkner et al. (2014) estimated that the effect of
these tides is known to an accuracy of $\approx 0.5\%$ or 0.2 mm/yr, making
this the limit to any unconventional contributions to $\dot{r}$.

Nonetheless, we consider the possibility that this calculation is seriously in
error such that the tidal contribution to $\dot{r}_{obs}$ is only
$\dot{r}_{tide}=10$ mm/yr, with SID contributing the remaining 28 mm/yr in
line with Equation 3. The total $\dot{r}$ would then be consistent with the
$\dot{r}_{obs}$ measured by LLR.

An important aspect of the SID contribution is that orbital velocities $v$ do
not change even though the radius $r\propto t$ (Equation 3). Since
$v^{2}=GM/r$ for a circular orbit around any object with mass $M$, its
gravitational parameter $GM\propto t$ in this model.111This partially cancels
the orbit expansion due to the non-standard velocity-dependent force, which
would otherwise be twice as much. Thus, in the absence of tides, SID predicts
that the Moon’s orbital angular frequency $\Omega\propto t^{-1}$ (equation 71
of Maeder & Bouvier, 1979). However, if the lunar orbit expands at the same
rate due to tides in a conventional context, Kepler’s Third Law implies that
$\Omega\propto t^{-3/2}$. Combining the tidal and SID contributions, we get
that

$\displaystyle\frac{\dot{\Omega}}{\Omega}~{}=~{}-\frac{\dot{r}_{{}_{SID}}+\frac{3}{2}\dot{r}_{tide}}{r}\,.$
(4)

Model | $\dot{r}_{tide}$, mm/yr | $\dot{r}_{{}_{SID}}$, mm/yr | $\dot{\Omega}$, $\arcsec$/century2
---|---|---|---
Standard | 38 | 0 | $-25.74$
SID | 28 | 10 | $-19.42$
Observed | Total $=38.05\pm 0.04$ | $-25.80\pm 0.03$
Table 1: The observed lunar recession rate is partly caused by tides raised on
the Earth. It might contain an additional term predicted by SID, but if so,
the effect of tides must be reduced since the total is tightly constrained by
lunar laser ranging (Folkner et al., 2014). Therefore, the scenarios predict
different angular deceleration rates for the lunar orbit (Equation 4). The
observed rate comes from their section 3c.

The conventional expectation is that $\dot{r}_{{}_{SID}}=0$ and
$\dot{r}_{tide}=38$ mm/yr, implying that
$\dot{\Omega}=-25.74\arcsec$/century2. If instead the precisely observed
$\dot{r}_{obs}$ consists of $\dot{r}_{tide}=10$ mm/yr and
$\dot{r}_{{}_{SID}}=28$ mm/yr, we should observe that
$\dot{\Omega}=-19.42\arcsec$/century2. These possibilities are summarized in
Table 1.

Space age observations of the Moon tell us that $\dot{\Omega}=-25.80\pm
0.03\arcsec$/century2 (section 3c of Folkner et al., 2014). This could be
affected by oscillatory perturbations from other planets, but fortunately a
much longer baseline is available if we also consider pre-telescopic records
of Solar and lunar eclipses over the past 2700 years (Stephenson & Morrison,
1995). This yields an estimated $\dot{\Omega}$ of ${-26\pm
1}\arcsec$/century2, with the error derived from their estimated timing
uncertainty of 0.9 s/century2 (their section 5a) and the fact that changing
$\dot{\Omega}$ by $0.5\arcsec$/century2 affects the eclipse timings by 0.46
s/century2 (their section 2b). Impressively, this section of their work found
that $\dot{\Omega}=-26.0\arcsec$/century2 fits the ancient records much better
than $\dot{\Omega}=-26.2\arcsec$/century2, correctly anticipating subsequent
refinements to the modern value of $\dot{\Omega}$.

Reconciling this estimate with SID requires changing $\dot{\Omega}$ by
$\Delta\dot{\Omega}=6\arcsec$/century2. The SID-predicted angular deceleration
of the Earth is only $0.94\arcsec$/century2, so SID-induced changes in its
orbital period have little effect on this discrepancy. Over ${T=2000}$ years,
the discrepancy amounts to a shift in the angular position of the Moon by
$T^{2}\Delta\dot{\Omega}/2=1200\arcsec$. Since the Moon takes 27.3 days to
orbit the Earth, it would take 36 minutes to rotate through ${1200\arcsec}$.
Ancient astronomers were well capable of noticing such a large time
difference, especially since some eclipses occurred close to sunrise or
sunset. Although there is some degeneracy between $\dot{\Omega}$ and changes
in the Earth’s rotation rate, the time of day provides a tight constraint on
the latter. The very fact that an eclipse occurred tightly constrains the
lunar orbit while the positions of background stars constrain Earth’s orbit
around the Sun. In this way, careful observations and record-keeping can break
the various degeneracies involved. Therefore, the SID theory is falsified at
extremely high significance by both modern and ancient observations.

##  3 Time-varying $G$ and stellar evolution

We have seen that SID predicts the gravitational parameter $GM\propto t$. If
this is caused by changes in $M$, then the masses of fundamental particles
would need to grow. In particular, changing the electron mass $\propto t$
would cause a similar change in the energy of hydrogen spectral lines, i.e.
the Rydberg constant would change. For much heavier elements, the expected
shift is not exactly the same due to special relativistic corrections that
become more important for atoms orbiting a more charged nucleus. This allowed
Karshenboim & Peik (2008) to place stringent constraints on possible time
evolution of the Rydberg constant, completely ruling out the possibility that
it grows $\propto t$.

Thus, the SID model requires $G\propto t$, implying a present value of
$\dot{G}/G=7.2\times 10^{-11}$/yr. This very slow change is in significant
tension with the precision pulsar timing constraint from PSR J0437-4715
($\dot{G}/G=\left(-0.5\pm 1.8\right)\times 10^{-11}$/yr, Verbiest et al.,
2008).

Over a Hubble time, the much larger expected change in $G$ would significantly
alter stellar evolution since the core pressure and temperature depend
sensitively on $G$. As a result, its expected variation is in strong tension
with astroseismic observations of the ancient star KIC 7970740
($\dot{G}/G=\left(1.2\pm 2.6\right)\times 10^{-12}$/yr, Bellinger &
Christensen-Dalsgaard, 2019).

More generally, since SID does not much change the age of the universe, it
would require significant revision to our understanding of main sequence
turnoffs in the oldest globular clusters (VandenBerg et al., 2013; Correnti et
al., 2018). We consider this unlikely because the nuclear processes in stars
are rather well understood and the turnoff stars have a mass of ${\approx
0.8M_{\odot}}$, making them only slightly less massive than the Sun.

Any time variation of $G$ would also change the Chandrasekhar mass $\propto
G^{-3/2}$ (Gaztañaga et al., 2002). This would affect the luminosities of Type
Ia supernovae, allowing observations of them to set limits on $\dot{G}$. Based
on this idea, the analysis of Mould & Uddin (2014) constrained $\dot{G}/G$ to
the range $\left(-3,7.3\right)\times 10^{-11}$/yr, marginally consistent with
the SID expectation of $7.2\times 10^{-11}$/yr. Future improvements to this
technique could yield more stringent constraints.

##  4 Conclusions

The SID cosmological model claims to successfully reproduce the acceleration-
dependent pattern of dynamical discrepancies in galaxies (the RAR, Maeder &
Gueorguiev, 2020). This is due to an extra term in the equations, which in the
Solar System implies that the Moon’s orbit around the Earth should expand at a
rate of $\dot{r}_{{}_{SID}}=28$ mm/yr (Section 2.1). Such an expansion of two-
body orbits has always been part of SID (Maeder, 1978; Maeder & Bouvier,
1979), with the latter work explicitly suggesting that this prediction should
be tested with laser ranging. The observed lunar recession rate of
$\dot{r}_{obs}={38.05\pm 0.04}$ mm/yr (Folkner et al., 2014) can be accounted
for by tidal dissipation in Earth’s oceans. For the SID theory to be correct,
there must be an extremely large error in our understanding of how this
process works, even though the uncertainty should be $\la 0.2$ mm/yr.
Supposing nonetheless that our understanding of terrestrial tides is
significantly in error, $\dot{r}_{obs}$ would have to consist of 28 mm/yr from
SID and 10 mm/yr from tides. This would cause the angular frequency of the
lunar orbit to decrease by $19.42\arcsec$/century2, contradicting the observed
rate of $25.80\pm 0.03\arcsec$/century2. This rate has remained stable for
several thousand years and is well explained using standard mechanics (Table
1).

Looking beyond the Earth-Moon system, SID predicts that the orbit of Saturn
expands at 104 m/yr (Section 2). This probably violates constraints from the
Cassini mission since the ranging accuracy is ${\approx 32}$ m and the orbiter
functioned for over a decade (table 11 of Viswanathan et al., 2017). Unlike
oscillatory perturbations from other planets, SID predicts a constant outwards
expansion that is difficult to mask by e.g. changing the mass of Jupiter.

In addition to predicting that orbits expand, SID also implies significant
time variation of any object’s gravitational parameter $GM$ (Section 3). Since
changing the masses of fundamental particles would violate precise laboratory
constraints, the most plausible interpretation is that the change occurs
solely in the gravitational constant $G$. This would substantially affect
stellar evolution because any ancient star would have fused hydrogen at a much
lower rate in the early universe, implying low-mass red giants must be
significantly older than in standard theory. However, the SID model also
implies the age of the universe is similar to that in $\Lambda$CDM (section
4.1 of Maeder & Gueorguiev, 2020).

In addition to these observational difficulties, it is not completely clear
that SID is in fact distinct from classical General Relativity $-$ the extra
degree of freedom may be a mathematical artefact without any physical effects
(Tsamis & Woodard, 1986; Jackiw & Pi, 2015). In any case, SID is unable to
satisfy Solar System constraints and contradicts stellar astrophysics to a
substantial extent. Therefore, some other way should be found to explain the
impressive successes of Milgromian dynamics in a covariant framework. The
recently proposed relativistic MOND theory of Skordis & Złośnik (2019) may be
a significant step in this direction.

## Acknowledgements

IB is supported by an Alexander von Humboldt postdoctoral research fellowship.
The authors are grateful to Andre Maeder for carefully presenting his work at
the Bonn Gravity2019 conference. They also thank Richard Woodard and Moti
Milgrom for helpful discussions, and the referee for comments which helped to
clarify this letter.

## Data availability

The data underlying this article are available in the article.

## References

  * Adelberger et al. (2017) Adelberger E. G., et al., 2017, Classical and Quantum Gravity, 34, 245008
  * Babcock (1939) Babcock H. W., 1939, Lick Observatory Bulletin, 19, 41
  * Banik (2019) Banik I., 2019, MNRAS, 487, 5291
  * Banik & Zhao (2018) Banik I., Zhao H., 2018, MNRAS, 480, 2660
  * Baumgardt et al. (2005) Baumgardt H., Grebel E. K., Kroupa P., 2005, MNRAS, 359, L1
  * Begeman et al. (1991) Begeman K. G., Broeils A. H., Sanders R. H., 1991, MNRAS, 249, 523
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Bellinger & Christensen-Dalsgaard (2019) Bellinger E. P., Christensen-Dalsgaard J., 2019, ApJL, 887, L1
  * Bílek et al. (2019) Bílek M., Müller O., Famaey B., 2019, A&A, 627, L1
  * Canuto et al. (1977) Canuto V., Adams P. J., Hsieh S. H., Tsiang E., 1977, Physical Review D, 16, 1643
  * Chowdhury (2019) Chowdhury A., 2019, MNRAS, 482, L99
  * Correnti et al. (2018) Correnti M., Gennaro M., Kalirai J. S., Cohen R. E., Brown T. M., 2018, ApJ, 864, 147
  * Derakhshani & Haghi (2014) Derakhshani K., Haghi H., 2014, ApJ, 783, 48
  * Desmond (2017a) Desmond H., 2017a, MNRAS, 464, 4160
  * Desmond (2017b) Desmond H., 2017b, MNRAS, 472, L35
  * Efstathiou et al. (1990) Efstathiou G., Sutherland W. J., Maddox S. J., 1990, Nature, 348, 705
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Folkner et al. (2014) Folkner W. M., Williams J. G., Boggs D. H., Park R. S., Kuchynka P., 2014, Interplanetary Network Progress Report, 42-196, 1
  * Gaztañaga et al. (2002) Gaztañaga E., García-Berro E., Isern J., Bravo E., Domínguez I., 2002, Physical Review D, 65, 023506
  * Ghari et al. (2019) Ghari A., Famaey B., Laporte C., Haghi H., 2019, A&A, 623, A123
  * Haghi et al. (2011) Haghi H., Baumgardt H., Kroupa P., 2011, A&A, 527, A33
  * Haghi et al. (2019a) Haghi H., et al., 2019a, MNRAS, 487, 2441
  * Haghi et al. (2019b) Haghi H., Amiri V., Hasani Zonoozi A., Banik I., Kroupa P., Haslbauer M., 2019b, ApJL, 884, L25
  * Hoof et al. (2020) Hoof S., Geringer-Sameth A., Trotta R., 2020, JCAP, 2, 012
  * Ibata et al. (2011a) Ibata R., Sollima A., Nipoti C., Bellazzini M., Chapman S. C., Dalessandro E., 2011a, ApJ, 738, 186
  * Ibata et al. (2011b) Ibata R., Sollima A., Nipoti C., Bellazzini M., Chapman S. C., Dalessandro E., 2011b, ApJ, 743, 43
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jackiw & Pi (2015) Jackiw R., Pi S.-Y., 2015, Physical Review D, 91, 067501
  * Karshenboim & Peik (2008) Karshenboim S. G., Peik E., 2008, European Physical Journal Special Topics, 163, 1
  * Kroupa (2012) Kroupa P., 2012, PASA, 29, 395
  * Kroupa (2015) Kroupa P., 2015, Canadian Journal of Physics, 93, 169
  * Kroupa et al. (2018a) Kroupa P., et al., 2018a, Nature Astronomy, 2, 925
  * Kroupa et al. (2018b) Kroupa P., et al., 2018b, Nature, 561, E4
  * Lelli et al. (2017) Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152
  * Li et al. (2018) Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3
  * Liu et al. (2017) Liu J., Chen X., Ji X., 2017, Nature Physics, 13, 212
  * Maeder (1978) Maeder A., 1978, A&A, 65, 337
  * Maeder & Bouvier (1979) Maeder A., Bouvier P., 1979, A&A, 73, 82
  * Maeder & Gueorguiev (2020) Maeder A., Gueorguiev V. G., 2020, MNRAS, 492, 2698
  * McGaugh (2011) McGaugh S. S., 2011, Physical Review Letters, 106, 121303
  * McGaugh et al. (2016) McGaugh S., Lelli F., Schombert J., 2016, Phys. Rev. Lett., 117, 201101
  * Milgrom (1983) Milgrom M., 1983, ApJ, 270, 365
  * Milgrom (1986) Milgrom M., 1986, ApJ, 302, 617
  * Milgrom (1999) Milgrom M., 1999, Phys. Lett. A, 253, 273
  * Milgrom (2009) Milgrom M., 2009, ApJ, 698, 1630
  * Mould & Uddin (2014) Mould J., Uddin S. A., 2014, PASA, 31, e015
  * Oman et al. (2015) Oman K. A., et al., 2015, MNRAS, 452, 3650
  * Ostriker & Peebles (1973) Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467
  * Ostriker & Steinhardt (1995) Ostriker J. P., Steinhardt P. J., 1995, Nature, 377, 600
  * Pazy (2013) Pazy E., 2013, Phys. Rev. D, 87, 084063
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, ApJ, 517, 565
  * Peters (1981) Peters P. C., 1981, American Journal of Physics, 49, 564
  * Pittordis & Sutherland (2019) Pittordis C., Sutherland W., 2019, MNRAS, 488, 4740
  * Planck Collaboration XIII (2016) Planck Collaboration XIII 2016, A&A, 594, A13
  * Renaud et al. (2016) Renaud F., Famaey B., Kroupa P., 2016, MNRAS, 463, 3637
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Rodrigues et al. (2018) Rodrigues D. C., Marra V., del Popolo A., Davari Z., 2018, Nature Astronomy, 2, 668
  * Rogstad & Shostak (1972) Rogstad D. H., Shostak G. S., 1972, ApJ, 176, 315
  * Rubin & Ford (1970) Rubin V. C., Ford Jr. W. K., 1970, ApJ, 159, 379
  * Sanders (2012a) Sanders R. H., 2012a, MNRAS, 419, L6
  * Sanders (2012b) Sanders R. H., 2012b, MNRAS, 422, L21
  * Sardone et al. (2019) Sardone A., Pisano D. J., Burke-Spolaor S., Mascoop J. L., Pol N., 2019, ApJL, 871, L31
  * Skordis & Złośnik (2019) Skordis C., Złośnik T., 2019, Physical Review D, 100, 104013
  * Smolin (2017) Smolin L., 2017, Physical Review D, 96, 083523
  * Stephenson & Morrison (1995) Stephenson F. R., Morrison L. V., 1995, Philosophical Transactions of the Royal Society of London Series A, 351, 165
  * Tsamis & Woodard (1986) Tsamis N. C., Woodard R. P., 1986, Annals of Physics, 168, 457
  * VandenBerg et al. (2013) VandenBerg D. A., Brogaard K., Leaman R., Casagrande L., 2013, ApJ, 775, 134
  * Verbiest et al. (2008) Verbiest J. P. W., et al., 2008, ApJ, 679, 675
  * Verlinde (2016) Verlinde E. P., 2016, SciPost Physics, 2, 16
  * Virgo & LIGO Collaborations (2017) Virgo & LIGO Collaborations 2017, Physical Review Letters, 119, 161101
  * Viswanathan et al. (2017) Viswanathan V., Fienga A., Gastineau M., Laskar J., 2017, NSTIM, 108
  * Wittenburg et al. (2020) Wittenburg N., Kroupa P., Famaey B., 2020, ApJ, 890, 173