# ETSformer: Exponential Smoothing Transformers for Time-series Forecasting

Gerald Woo1 2, Chenghao Liu1, Doyen Sahoo1, Akshat Kumar2 & Steven Hoi1
1Salesforce Research Asia, 2Singapore Management University
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Transformers have been actively studied for time-series forecasting in recent
years. While often showing promising results in various scenarios, traditional
Transformers are not designed to fully exploit the characteristics of time-
series data and thus suffer some fundamental limitations, e.g., they are
generally not decomposable or interpretable, and are neither effective nor
efficient for long-term forecasting. In this paper, we propose ETSformer, a
novel time-series Transformer architecture, which exploits the principle of
exponential smoothing in improving Transformers for time-series forecasting.
In particular, inspired by the classical exponential smoothing methods in
time-series forecasting, we propose the novel exponential smoothing attention
(ESA) and frequency attention (FA) to replace the self-attention mechanism in
vanilla Transformers, thus improving both accuracy and efficiency. Based on
these, we redesign the Transformer architecture with modular decomposition
blocks such that it can learn to decompose the time-series data into
interpretable time-series components such as level, growth and seasonality.
Extensive experiments on various time-series benchmarks validate the efficacy
and advantages of the proposed method. Code is available at
https://github.com/salesforce/ETSformer.

##  1 Introduction

Transformer models have achieved great success in the fields of NLP [33, 8]
and CV [4, 9] in recent times. The success is widely attributed to its self-
attention mechanism which is able to explicitly model both short and long
range dependencies adaptively via the pairwise query-key interaction. Owing to
their powerful capability to model sequential data, Transformer-based
architectures [20, 38, 41, 37, 40] have been actively explored for the time-
series forecasting, especially for the more challenging Long Sequence Time-
series Forecasting (LSTF) task. While showing promising results, it is still
quite challenging to extract salient temporal patterns and thus make accurate
long-term forecasts for large-scale data. This is because time-series data is
usually noisy and non-stationary. Without incorporating appropriate knowledge
about time-series structures [1, 31, 13], it is prone to learning the spurious
dependencies and lacks interpretability.

Moreover, the use of content-based, dot-product attention in Transformers is
not effective in detecting essential temporal dependencies for two reasons.
(1) Firstly, time-series data is usually assumed to be generated by a
conditional distribution over past observations, with the dependence between
observations weakening over time [23, 17]. Therefore, neighboring data points
have similar values, and recent tokens should be given a higher weight111An
assumption further supported by the success of classical exponential smoothing
methods and ARIMA model selection methods tending to select small lags. when
measuring their similarity [13, 14]. This indicates that attention measured by
a relative time lag is more effective than that measured by the similarity of
the content when modeling time-series. (2) Secondly, many real world time-
series display strong seasonality – patterns in time-series which repeat with
a fixed period. Automatically extracting seasonal patterns has been proved to
be critical for the success of forecasting [6, 5, 36]. However, the vanilla
attention mechanism is unlikely able to learn these required periodic
dependencies without any in-built prior structure.

Figure 1: Illustration demonstrating how ETSformer generates forecasts via a
decomposition (of intermediate representations) into seasonal and trend
components. The seasonal component extracts salient periodic patterns and
extrapolates them. The trend component which is a combination of the level and
growth terms, first estimates the current level of the time-series, and
subsequently adds a damped growth term to generate trend forecasts.

To address these limitations, we propose ETSformer, an effective and efficient
Transformer architecture for time-series forecasting, inspired by exponential
smoothing methods [13] and illustrated in Figure 1. First of all, ETSformer
incorporates inductive biases of time-series structures by performing a layer-
wise level, growth, and seasonal decomposition. By leveraging the high
capacities of deep architectures and an effective residual learning scheme,
ETSformer is able to extract a series of latent growth and seasonal patterns
and model their complex dependencies. Secondly, ETSformer introduces a novel
Exponential Smoothing Attention (ESA) and Frequency Attention (FA) to replace
vanilla attention. In particular, ESA constructs attention scores based on the
relative time lag to the query, and achieves $\mathcal{O}(L\log L)$ complexity
for the length-$L$ lookback window and demonstrates powerful capability in
modeling the growth component. FA leverages the Fourier transformation to
extract the dominating seasonal patterns by selecting the Fourier bases with
the $K$ largest amplitudes in frequency domain, and also achieves
$\mathcal{O}(L\log L)$ complexity. Finally, the predicted forecast is a
composition of level, trend, and seasonal components, which makes it human
interpretable. We conduct extensive empirical analysis and show that ETSformer
achieves state-of-the-art performance by outperforming competing approaches
over 6 real world datasets on both the multivariate and univariate settings,
and also visualize the time-series components to verify its interpretability.

##  2 Related Work

Transformer based deep forecasting. Inspired by the success of Transformers in
CV and NLP, Transformer-based time-series forecasting models have been
actively studied recently. LogTrans [20] introduces local context to
Transformer models via causal convolutions in the query-key projection layer,
and propose the LogSparse attention to reduce complexity to $\mathcal{O}(L\log
L)$. Informer [41] extends the Transformer by proposing the ProbSparse
attention and distillation operation to achieve $\mathcal{O}(L\log L)$
complexity. AST [38] leverages a sparse normalization transform,
$\alpha\mathchar 45\relax\mathrm{entmax}$, to implement a sparse attention
layer. It further incorporates an adversarial loss to mitigate the adverse
effect of error accumulation in inference. Similar to our work that
incorporates prior knowledge of time-series structure, Autoformer [37]
introduces the Auto-Correlation attention mechanism which focuses on sub-
series based similarity and is able to extract periodic patterns. Yet, their
implementation of series decomposition which performs de-trending via a simple
moving average over the input signal without any learnable parameters is
arguably a simplified assumption, insufficient to appropriately model complex
trend patterns. ETSformer on the other hand, decomposes the series by de-
seasonalization as seasonal patterns are more identifiable and easier to
detect [7]. Furthermore, the Auto-Correlation mechanism fails to attend to
information from the local context (i.e. forecast at $t+1$ is not dependent on
$t,t-1$, etc.) and does not separate the trend component into level and growth
components, which are both crucial for modeling trend patterns. Lastly,
similar to previous work, their approach is highly reliant on manually
designed dynamic time-dependent covariates (e.g. month-of-year, day-of-week),
while ETSformer is able to automatically learn and extract seasonal patterns
from the time-series signal directly.

Attention Mechanism. The self-attention mechanism in Transformer models has
recently received much attention, its necessity has been greatly investigated
in attempts to introduce more flexibility and reduce computational cost.
Synthesizer [29] empirically studies the importance of dot-product
interactions, and show that a randomly initialized, learnable attention
mechanisms with or without token-token dependencies can achieve competitive
performance with vanilla self-attention on various NLP tasks. [39] utilizes an
unparameterized Gaussian distribution to replace the original attention
scores, concluding that the attention distribution should focus on a certain
local window and can achieve comparable performance. [25] replaces attention
with fixed, non-learnable positional patterns, obtaining competitive
performance on NMT tasks. [19] replaces self-attention with a non-learnable
Fourier Transform and verifies it to be an effective mixing mechanism. While
our proposed ESA shares the spirit of designing attention mechanisms that are
not dependent on pair-wise query-key interactions, our work is inspired by
exploiting the characteristics of time-series and is an early attempt to
utilize prior knowledge of time-series for tackling the time-series
forecasting tasks.

##  3 Preliminaries and Background

Problem Formulation Let ${\bm{x}}_{t}\in\mathbb{R}^{m}$ denote an observation
of a multivariate time-series at time step $t$. Given a lookback window
${\bm{X}}_{t-L:t}=[{\bm{x}}_{t-L},\ldots,{\bm{x}}_{t-1}]$, we consider the
task of predicting future values over a horizon,
${\bm{X}}_{t:t+H}=[{\bm{x}}_{t},\ldots,{\bm{x}}_{t+H-1}]$. We denote
$\hat{{\bm{X}}}_{t:t+H}$ as the point forecast of ${\bm{X}}_{t:t+H}$. Thus,
the goal is to learn a forecasting function
$\hat{{\bm{X}}}_{t:t+H}=f({\bm{X}}_{t-L:t})$ by minimizing some loss function
$\mathcal{L}:\mathbb{R}^{H\times m}\times\mathbb{R}^{H\times m}\to\mathbb{R}$.

Exponential Smoothing We instantiate exponential smoothing methods [13] in the
univariate forecasting setting. They assume that time-series can be decomposed
into seasonal and trend components, and trend can be further decomposed into
level and growth components. Specifically, a commonly used model is the
additive Holt-Winters’ method [12, 35], which can be formulated as:

Level $\displaystyle:e_{t}=\alpha(x_{t}-s_{t-p})+(1-\alpha)(e_{t-1}+b_{t-1})$
Growth $\displaystyle:b_{t}=\beta(e_{t}-e_{t-1})+(1-\beta)b_{t-1}$ Seasonal
$\displaystyle:s_{t}=\gamma(x_{t}-e_{t})+(1-\gamma)s_{t-p}$ Forecasting
$\displaystyle:\hat{x}_{t+h|t}=e_{t}+hb_{t}+s_{t+h-p}$ (1)

where $p$ is the period of seasonality, and $\hat{x}_{t+h|t}$ is the $h$-steps
ahead forecast. The above equations state that the $h$-steps ahead forecast is
composed of the last estimated level $e_{t}$, incrementing it by $h$ times the
last growth factor, $b_{t}$, and adding the last available seasonal factor
$s_{t+h-p}$. Specifically, the level smoothing equation is formulated as a
weighted average of the seasonally adjusted observation $(x_{t}-s_{t-p})$ and
the non-seasonal forecast, obtained by summing the previous level and growth
$(e_{t-1}+b_{t-1})$. The growth smoothing equation is implemented by a
weighted average between the successive difference of the (de-seasonalized)
level, $(e_{t}-e_{t-1})$, and the previous growth, $b_{t-1}$. Finally, the
seasonal smoothing equation is a weighted average between the difference of
observation and (de-seasonalized) level, $(x_{t}-e_{t})$, and the previous
seasonal index $s_{t-p}$. The weighted average of these three equations are
controlled by the smoothing parameters $\alpha$, $\beta$ and $\gamma$,
respectively.

A widely used modification of the additive Holt-Winters’ method is to allow
the damping of trends, which has been proved to produce robust multi-step
forecasts [28, 22]. The forecast with damping trend can be rewritten as:

$\displaystyle\hat{x}_{t+h|t}=e_{t}+(\phi+\phi^{2}+\dots+\phi^{h})b_{t}+s_{t+h-p},$
(2)

where the growth is damped by a factor of $\phi$. If $\phi=1$, it degenerates
to the vanilla forecast. For $0<\phi<1$, as $h\rightarrow\infty$ this growth
component approaches an asymptote given by $\phi b_{t}/(1-\phi)$.

##  4 ETSformer

Figure 2: ETSformer model architecture.

In this section, we redesign the classical Transformer architecture into an
exponential smoothing inspired encoder-decoder architecture specialized for
tackling the time-series forecasting problem. Our architecture design
methodology relies on three key principles: (1) the architecture leverages the
stacking of multiple layers to progressively extract a series of level,
growth, and seasonal representations from the intermediate latent residual;
(2) following the spirit of exponential smoothing, we extract the salient
seasonal patterns while modeling level and growth components by assigning
higher weight to recent observations; (3) the final forecast is a composition
of level, growth, and seasonal components making it human interpretable. We
now expound how our ETSformer architecture encompasses these principles.

###  4.1 Overall Architecture

Figure 2 illustrates the overall encoder-decoder architecture of ETSformer. At
each layer, the encoder is designed to iteratively extract growth and seasonal
latent components from the lookback window. The level is then extracted in a
similar fashion to classical level smoothing in Section 3. These extracted
components are then fed to the decoder to further generate the final $H$-step
ahead forecast via a composition of level, growth, and seasonal forecasts,
which is defined:

$\displaystyle\hat{{\bm{X}}}_{t:t+H}={\bm{E}}_{t:t+H}+\mathrm{Linear}\Big{(}\sum^{N}_{n=1}({\bm{B}}_{t:t+H}^{(n)}+{\bm{S}}_{t:t+H}^{(n)})\Big{)},$
(3)

where ${\bm{E}}_{t:t+H}\in\mathbb{R}^{H\times m}$, and
${\bm{B}}_{t:t+H}^{(n)},{\bm{S}}_{t:t+H}^{(n)}\in\mathbb{R}^{H\times d}$
represent the level forecasts, and the growth and seasonal latent
representations of each time step in the forecast horizon, respectively. The
superscript represents the stack index, for a total of $N$ encoder stacks.
Note that $\mathrm{Linear}(\cdot):\mathbb{R}^{d}\to\mathbb{R}^{m}$ operates
element-wise along each time step, projecting the extracted growth and
seasonal representations from latent to observation space.

####  4.1.1 Input Embedding

Raw signals from the lookback window are mapped to latent space via the input
embedding module, defined by
${\bm{Z}}_{t-L:t}^{(0)}={\bm{E}}_{t-L:t}^{(0)}=\mathrm{Conv}({\bm{X}}_{t-L:t}),$
where $\mathrm{Conv}$ is a temporal convolutional filter with kernel size 3,
input channel $m$ and output channel $d$. In contrast to prior work [20, 38,
37, 41], the inputs of ETSformer do not rely on any other manually designed
dynamic time-dependent covariates (e.g. month-of-year, day-of-week) for both
the lookback window and forecast horizon. This is because the proposed
Frequency Attention module (details in Section 4.2.2) is able to automatically
uncover these seasonal patterns, which renders it more applicable for
challenging scenarios without these discriminative covariates and reduces the
need for feature engineering.

####  4.1.2 Encoder

The encoder focuses on extracting a series of latent growth and seasonality
representations in a cascaded manner from the lookback window. To achieve this
goal, traditional methods rely on the assumption of additive or multiplicative
seasonality which has limited capability to express complex patterns beyond
these assumptions. Inspired by [24, 10], we leverage residual learning to
build an expressive, deep architecture to characterize the complex intrinsic
patterns. Each layer can be interpreted as sequentially analyzing the input
signals. The extracted growth and seasonal signals are then removed from the
residual and undergo a nonlinear transformation before moving to the next
layer. Each encoder layer takes as input the residual from the previous
encoder layer ${\bm{Z}}_{t-L:t}^{(n-1)}$ and emits
${\bm{Z}}_{t-L:t}^{(n)},{\bm{B}}_{t-L:t}^{(n)},{\bm{S}}_{t-L:t}^{(n)}$, the
residual, latent growth, and seasonal representations for the lookback window
via the Multi-Head Exponential Smoothing Attention (MH-ESA) and Frequency
Attention (FA) modules (detailed description in Section 4.2). The following
equations formalizes the overall pipeline in each encoder layer, and for ease
of exposition, we use the notation $\vcentcolon=$ for a variable update.
$\displaystyle\textrm{Seasonal:}\quad{\bm{S}}_{t-L:t}^{(n)}$
$\displaystyle=\mathrm{FA}_{t-L:t}({\bm{Z}}_{t-L:t}^{(n-1)})$
$\displaystyle{\bm{Z}}_{t-L:t}^{(n-1)}$
$\displaystyle\vcentcolon={\bm{Z}}_{t-L:t}^{(n-1)}-{\bm{S}}_{t-L:t}^{(n)}$
$\displaystyle\textrm{Growth:}\quad{\bm{B}}_{t-L:t}^{(n)}$
$\displaystyle=\mathrm{MH}\mathchar
45\relax\mathrm{ESA}({\bm{Z}}_{t-L:t}^{(n-1)})$
$\displaystyle{\bm{Z}}_{t-L:t}^{(n-1)}$
$\displaystyle\vcentcolon=\mathrm{LN}({\bm{Z}}_{t-L:t}^{(n-1)}-{\bm{B}}_{t-L:t}^{(n)})$
$\displaystyle{\bm{Z}}_{t-L:t}^{(n)}$
$\displaystyle=\mathrm{LN}({\bm{Z}}_{t-L:t}^{(n-1)}+\mathrm{FF}({\bm{Z}}_{t-L:t}^{(n-1)}))$

$\mathrm{LN}$ is layer normalization [2],
$\mathrm{FF}(x)=\mathrm{Linear}(\sigma(\mathrm{Linear}(x)))$ is a position-
wise feedforward network [33] and $\sigma(\cdot)$ is the sigmoid function.

Level Module Given the latent growth and seasonal representations from each
layer, we extract the level at each time step $t$ in the lookback window in a
similar way as the level smoothing equation in Section 3. Formally, the
adjusted level is a weighted average of the current (de-seasonalized) level
and the level-growth forecast from the previous time step $t-1$. It can be
formulated as:

$\displaystyle{\bm{E}}_{t}^{(n)}=\bm{\alpha}*\Big{(}{\bm{E}}_{t}^{(n-1)}-\mathrm{Linear}({\bm{S}}_{t}^{(n)})\Big{)}+(1-\bm{\alpha})*\Big{(}{\bm{E}}_{t-1}^{(n)}+\mathrm{Linear}({\bm{B}}_{t-1}^{(n)})\Big{)},$

where $\bm{\alpha}\in\mathbb{R}^{m}$ is a learnable smoothing parameter, $*$
is an element-wise multiplication term, and
$\mathrm{Linear}(\cdot):\mathbb{R}^{d}\to\mathbb{R}^{m}$ maps representations
to observation space. Finally, the extracted level in the last layer
${\bm{E}}_{t-L:t}^{(N)}$ can be regarded as the corresponding level for the
lookback window. We show in Section A.3 that this recurrent exponential
smoothing equation can also be efficiently evaluated using the efficient
$\mathcal{A}_{\mathrm{ES}}$ algorithm (Algorithm 1) with an auxiliary term.

####  4.1.3 Decoder

The decoder is tasked with generating the $H$-step ahead forecasts. As shown
in Equation 3, the final forecast is a composition of level forecasts
${\bm{E}}_{t:t+H}$, growth representations ${\bm{B}}_{t:t+H}^{(n)}$ and
seasonal representations ${\bm{S}}_{t:t+H}^{(n)}$ in the forecast horizon. It
comprises $N$ Growth + Seasonal (G+S) Stacks, and a Level Stack. The G+S Stack
consists of the Growth Damping (GD) and FA blocks, which leverage
${\bm{B}}_{t}^{(n)}$, ${\bm{S}}_{t-L:t}^{(n)}$ to predict
${\bm{B}}_{t:t+H}^{(n)}$, ${\bm{S}}_{t:t+H}^{(n)}$, respectively.
$\displaystyle\text{Growth:}\quad{\bm{B}}_{t:t+H}^{(n)}=\mathrm{TD}({\bm{B}}_{t}^{(n)})$
$\displaystyle\text{Seasonal:}\quad{\bm{S}}_{t:t+H}^{(n)}=\mathrm{FA}_{t:t+H}({\bm{S}}_{t-L:t}^{(n)})$

To obtain the level in the forecast horizon, the Level Stack repeats the level
in the last time step $t$ along the forecast horizon. It can be defined as
${\bm{E}}_{t:t+H}=\mathrm{Repeat}_{H}({\bm{E}}_{t}^{(N)})=[{\bm{E}}_{t}^{(N)},\ldots,{\bm{E}}_{t}^{(N)}]$,
with $\mathrm{Repeat}_{H}(\cdot):\mathbb{R}^{1\times m}\to\mathbb{R}^{H\times
m}$.

Growth Damping To obtain the growth representation in the forecast horizon, we
follow the idea of trend damping in Equation 2 to make robust multi-step
forecast. Thus, the trend representations can be formulated as:

$\displaystyle\mathrm{TD}({\bm{B}}_{t}^{(n)})_{j}$
$\displaystyle=\sum_{i=1}^{j}\gamma^{i}{\bm{B}}_{t}^{(n)},$
$\displaystyle\mathrm{TD}({\bm{B}}_{t-L:t}^{(n)})$
$\displaystyle=[\mathrm{TD}({\bm{B}}_{t}^{(n)})_{t},\ldots,\mathrm{TD}({\bm{B}}_{t}^{(n)})_{t+H-1}],$

where $0<\gamma<1$ is the damping parameter which is learnable, and in
practice, we apply a multi-head version of trend damping by making use of
$n_{h}$ damping parameters. Similar to the implementation for level forecast
in the Level Stack, we only use the last trend representation in the lookback
window ${\bm{B}}_{t}^{(n)}$ to forecast the trend representation in the
forecast horizon.

###  4.2 Exponential Smoothing Attention and Frequency Attention Mechanism

(a) Full Attention (2017)

(b) Sparse Attention (2020, 2021)

(c) Log-sparse Attention (2019)

(d) Autocorrelation Attention (2021)

(e) Exponential Smoothing Attention (Ours)

(f) Frequency Attention (Ours)

Figure 3: Comparison between different attention mechanisms. (a) Full, (b)
Sparse, and (c) Log-sparse Attentions are adaptive mechanisms, where the green
circles represent the attention weights adaptively calculated by a point-wise
dot-product query, and depends on various factors including the time-series
value, additional covariates (e.g. positional encodings, time features, etc.).
(d) Autocorrelation attention considers sliding dot-product queries to
construct attention weights for each rolled input series. We introduce (e)
Exponential Smoothing Attention (ESA) and (f) Frequency Attention (FA). ESA
directly computes attention weights based on the relative time lag, without
considering the input content, while FA attends to patterns which dominate
with large magnitudes in the frequency domain.

Considering the ineffectiveness of existing attention mechanisms in handling
time-series data, we develop the Exponential Smoothing Attention (ESA) and
Frequency Attention (FA) mechanisms to extract latent growth and seasonal
representations. ESA is a non-adaptive, learnable attention scheme with an
inductive bias to attend more strongly to recent observations by following an
exponential decay, while FA is a non-learnable attention scheme, that
leverages Fourier transformation to select dominating seasonal patterns. A
comparison between existing work and our proposed ESA and FA is illustrated in
Figure 3.

####  4.2.1 Exponential Smoothing Attention

Vanilla self-attention can be regarded as a weighted combination of an input
sequence, where the weights are normalized alignment scores measuring the
similarity between input contents [32]. Inspired by the exponential smoothing
in Section 3, we aim to assign a higher weight to recent observations. It can
be regarded as a novel form of attention whose weights are computed by the
relative time lag, rather than input content. Thus, the ESA mechanism can be
defined as $\mathcal{A}_{\mathrm{ES}}:\mathbb{R}^{L\times
d}\to\mathbb{R}^{L\times d}$, where
$\mathcal{A}_{\mathrm{ES}}({\bm{V}})_{t}\in\mathbb{R}^{d}$ denotes the $t$-th
row of the output matrix, representing the token corresponding to the $t$-th
time step. Its exponential smoothing formula can be further written as:

$\displaystyle\mathcal{A}_{\mathrm{ES}}({\bm{V}})_{t}$
$\displaystyle=\alpha{\bm{V}}_{t}+(1-\alpha)\mathcal{A}_{\mathrm{ES}}({\bm{V}})_{t-1}=\sum_{j=0}^{t-1}\alpha(1-\alpha)^{j}{\bm{V}}_{t-j}+(1-\alpha)^{t}{\bm{v}}_{0},$

where $0<\alpha<1$ and ${\bm{v}}_{0}$ are learnable parameters known as the
smoothing parameter and initial state respectively.

Efficient $\mathcal{A}_{\mathrm{ES}}$ algorithm The straightforward
implementation of the ESA mechanism by constructing the attention matrix,
${\bm{A}}_{\mathrm{ES}}$ and performing a matrix multiplication with the input
sequence (detailed algorithm in Section A.4) results in an
$\mathcal{O}(L^{2})$ computational complexity.

$\displaystyle\mathcal{A}_{\mathrm{ES}}({\bm{V}})=\begin{bmatrix}\mathcal{A}_{\mathrm{ES}}({\bm{V}})_{1}\\\
\vdots\\\
\mathcal{A}_{\mathrm{ES}}({\bm{V}})_{L}\end{bmatrix}={\bm{A}}_{\mathrm{ES}}\cdot\begin{bmatrix}{\bm{v}}_{0}^{T}\\\
{\bm{V}}\end{bmatrix},$

Yet, we are able to achieve an efficient algorithm by exploiting the unique
structure of the exponential smoothing attention matrix,
${\bm{A}}_{\mathrm{ES}}$, which is illustrated in Section A.1. Each row of the
attention matrix can be regarded as iteratively right shifting with padding
(ignoring the first column). Thus, a matrix-vector multiplication can be
computed with a cross-correlation operation, which in turn has an efficient
fast Fourier transform implementation [21]. The full algorithm is described in
Algorithm 1, Section A.2, achieving an $\mathcal{O}(L\log L)$ complexity.

Multi-Head Exponential Smoothing Attention (MH-ESA) We use
$\mathcal{A}_{\mathrm{ES}}$ as a basic building block, and develop the Multi-
Head Exponential Smoothing Attention to extract latent growth representations.
Formally, we obtain the growth representations by taking the successive
difference of the residuals.

$\displaystyle\tilde{{\bm{Z}}}_{t-L:t}^{(n)}$
$\displaystyle=\mathrm{Linear}({\bm{Z}}_{t-L:t}^{(n-1)}),$
$\displaystyle{\bm{B}}_{t-L:t}^{(n)}$ $\displaystyle=\mathrm{MH}\mathchar
45\relax\mathcal{A}_{\mathrm{ES}}(\tilde{{\bm{Z}}}_{t-L:t}^{(n)}-[\tilde{{\bm{Z}}}_{t-L:t-1}^{(n)},{\bm{v}}_{0}^{(n)}]),$
$\displaystyle{\bm{B}}_{t-L:t}^{(n)}$
$\displaystyle\coloneqq\mathrm{Linear}({\bm{B}}_{t-L:t}^{(n)}),$

where $\mathrm{MH}\mathchar 45\relax\mathcal{A}_{\mathrm{ES}}$ is a multi-head
version of $\mathcal{A}_{\mathrm{ES}}$ and ${\bm{v}}_{0}^{(n)}$ is the initial
state from the ESA mechanism.

####  4.2.2 Frequency Attention

The goal of identifying and extracting seasonal patterns from the lookback
window is twofold. Firstly, it can be used to perform de-seasonalization on
the input signals such that downstream components are able to focus on
modeling the level and growth information. Secondly, we are able to
extrapolate the seasonal patterns to build representations for the forecast
horizon. The main challenge is to automatically identify seasonal patterns.
Fortunately, the use of power spectral density estimation for periodicity
detection has been well studied [34]. Inspired by these methods, we leverage
the discrete Fourier transform (DFT, details in Appendix B) to develop the FA
mechanism to extract dominant seasonal patterns.

Specifically, FA first decomposes input signals into their Fourier bases via a
DFT along the temporal dimension,
$\mathcal{F}({\bm{Z}}_{t-L:t}^{(n-1)})\in\mathbb{C}^{F\times d}$ where
$F=\lfloor L/2\rfloor+1$, and selects bases with the $K$ largest amplitudes.
An inverse DFT is then applied to obtain the seasonality pattern in time
domain. Formally, this is given by the following equations:

$\displaystyle\bm{\Phi}_{k,i}=\phi\Big{(}\mathcal{F}({\bm{Z}}_{t-L:t}^{(n-1)})_{k,i}\Big{)},\quad{\bm{A}}_{k,i}=\Big{|}\mathcal{F}({\bm{Z}}_{t-L:t}^{(n-1)})_{k,i}\Big{|},$
$\displaystyle\kappa_{i}^{(1)},\ldots,\kappa_{i}^{(K)}=\mathop{\arg\mathrm{Top\mathchar
45\relax K}}_{k\in\\{2,\ldots,F\\}}~\Big{\\{}{\bm{A}}_{k,i}\Big{\\}},$
$\displaystyle{\bm{S}}_{j,i}^{(n)}=\sum_{k=1}^{K}{\bm{A}}_{\kappa_{i}^{(k)},i}\Big{[}\cos(2\pi
f_{\kappa_{i}^{(k)}}j+\bm{\Phi}_{\kappa_{i}^{(k)},i})+\cos(2\pi\bar{f}_{\kappa_{i}^{(k)}}j+\bar{\bm{\Phi}}_{\kappa_{i}^{(k)},i})\Big{]},$
(4)

where $\bm{\Phi}_{k,i},{\bm{A}}_{k,i}$ are the phase/amplitude of the $k$-th
frequency for the $i$-th dimension, $\mathop{\arg\mathrm{Top\mathchar 45\relax
K}}$ returns the arguments of the top $K$ amplitudes, $K$ is a hyperparameter,
$f_{k}$ is the Fourier frequency of the corresponding index, and
$\bar{f}_{k},\bar{\bm{\Phi}}_{k,i}$ are the Fourier frequency/amplitude of the
corresponding conjugates.

Finally, the latent seasonal representation of the $i$-th dimension for the
lookback window is formulated as
${\bm{S}}_{t-L:t,i}^{(n)}=[{\bm{S}}_{t-L,i}^{(n)},\ldots,{\bm{S}}_{t-1,i}^{(n)}]$.
For the forecast horizon, the FA module extrapolates beyond the lookback
window via,
${\bm{S}}_{t:t+H,i}^{(n)}=[{\bm{S}}_{t,i}^{(n)},\ldots,{\bm{S}}_{t+H-1,i}^{(n)}]$.
Since $K$ is a hyperparameter typically chosen for small values, the
complexity for the FA mechanism is similarly $\mathcal{O}(L\log L)$.

##  5 Experiments

This section presents extensive empirical evaluations on the LSTF task over 6
real world datasets, ETT, ECL, Exchange, Traffic, Weather, and ILI, coming
from a variety of application areas (details in Appendix D) for both
multivariate and univariate settings. This is followed by an ablation study of
the various contributing components, and interpretability experiments of our
proposed model. An additional analysis on computational efficiency can be
found in Appendix H for space. For the main benchmark, datasets are split into
train, validation, and test sets chronologically, following a 60/20/20 split
for the ETT datasets and 70/10/20 split for other datasets. Inputs are zero-
mean normalized and we use MSE and MAE as evaluation metrics. Further details
on implementation and hyperparameters can be found in Appendix C.

###  5.1 Results

Table 1: Multivariate forecasting results over various forecast horizons. Best
results are bolded, and second best results are underlined.

Methods | ETSformer | Autoformer | Informer | LogTrans | Reformer | LSTnet | LSTM
---|---|---|---|---|---|---|---
Metrics | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE
ETTm2 | 96 | 0.189 | 0.280 | 0.255 | 0.339 | 0.365 | 0.453 | 0.768 | 0.642 | 0.658 | 0.619 | 3.142 | 1.365 | 2.041 | 1.073
192 | 0.253 | 0.319 | 0.281 | 0.340 | 0.533 | 0.563 | 0.989 | 0.757 | 1.078 | 0.827 | 3.154 | 1.369 | 2.249 | 1.112
336 | 0.314 | 0.357 | 0.339 | 0.372 | 1.363 | 0.887 | 1.334 | 0.872 | 1.549 | 0.972 | 3.160 | 1.369 | 2.568 | 1.238
720 | 0.414 | 0.413 | 0.422 | 0.419 | 3.379 | 1.388 | 3.048 | 1.328 | 2.631 | 1.242 | 3.171 | 1.368 | 2.720 | 1.287
ECL | 96 | 0.187 | 0.304 | 0.201 | 0.317 | 0.274 | 0.368 | 0.258 | 0.357 | 0.312 | 0.402 | 0.680 | 0.645 | 0.375 | 0.437
192 | 0.199 | 0.315 | 0.222 | 0.334 | 0.296 | 0.386 | 0.266 | 0.368 | 0.348 | 0.433 | 0.725 | 0.676 | 0.442 | 0.473
336 | 0.212 | 0.329 | 0.231 | 0.338 | 0.300 | 0.394 | 0.280 | 0.380 | 0.350 | 0.433 | 0.828 | 0.727 | 0.439 | 0.473
720 | 0.233 | 0.345 | 0.254 | 0.361 | 0.373 | 0.439 | 0.283 | 0.376 | 0.340 | 0.420 | 0.957 | 0.811 | 0.980 | 0.814
Exchange | 96 | 0.085 | 0.204 | 0.197 | 0.323 | 0.847 | 0.752 | 0.968 | 0.812 | 1.065 | 0.829 | 1.551 | 1.058 | 1.453 | 1.049
192 | 0.182 | 0.303 | 0.300 | 0.369 | 1.204 | 0.895 | 1.040 | 0.851 | 1.188 | 0.906 | 1.477 | 1.028 | 1.846 | 1.179
336 | 0.348 | 0.428 | 0.509 | 0.524 | 1.672 | 1.036 | 1.659 | 1.081 | 1.357 | 0.976 | 1.507 | 1.031 | 2.136 | 1.231
720 | 1.025 | 0.774 | 1.447 | 0.941 | 2.478 | 1.310 | 1.941 | 1.127 | 1.510 | 1.016 | 2.285 | 1.243 | 2.984 | 1.427
Traffic | 96 | 0.607 | 0.392 | 0.613 | 0.388 | 0.719 | 0.391 | 0.684 | 0.384 | 0.732 | 0.423 | 1.107 | 0.685 | 0.843 | 0.453
192 | 0.621 | 0.399 | 0.616 | 0.382 | 0.696 | 0.379 | 0.685 | 0.390 | 0.733 | 0.420 | 1.157 | 0.706 | 0.847 | 0.453
336 | 0.622 | 0.396 | 0.622 | 0.337 | 0.777 | 0.420 | 0.733 | 0.408 | 0.742 | 0.420 | 1.216 | 0.730 | 0.853 | 0.455
720 | 0.632 | 0.396 | 0.660 | 0.408 | 0.864 | 0.472 | 0.717 | 0.396 | 0.755 | 0.423 | 1.481 | 0.805 | 1.500 | 0.805
Weather | 96 | 0.197 | 0.281 | 0.266 | 0.336 | 0.300 | 0.384 | 0.458 | 0.490 | 0.689 | 0.596 | 0.594 | 0.587 | 0.369 | 0.406
192 | 0.237 | 0.312 | 0.307 | 0.367 | 0.598 | 0.544 | 0.658 | 0.589 | 0.752 | 0.638 | 0.560 | 0.565 | 0.416 | 0.435
336 | 0.298 | 0.353 | 0.359 | 0.359 | 0.578 | 0.523 | 0.797 | 0.652 | 0.639 | 0.596 | 0.597 | 0.587 | 0.455 | 0.454
720 | 0.352 | 0.388 | 0.419 | 0.419 | 1.059 | 0.741 | 0.869 | 0.675 | 1.130 | 0.792 | 0.618 | 0.599 | 0.535 | 0.520
ILI | 24 | 2.527 | 1.020 | 3.483 | 1.287 | 5.764 | 1.677 | 4.480 | 1.444 | 4.400 | 1.382 | 6.026 | 1.770 | 5.914 | 1.734
36 | 2.615 | 1.007 | 3.103 | 1.148 | 4.755 | 1.467 | 4.799 | 1.467 | 4.783 | 1.448 | 5.340 | 1.668 | 6.631 | 1.845
48 | 2.359 | 0.972 | 2.669 | 1.085 | 4.763 | 1.469 | 4.800 | 1.468 | 4.832 | 1.465 | 6.080 | 1.787 | 6.736 | 1.857
60 | 2.487 | 1.016 | 2.770 | 1.125 | 5.264 | 1.564 | 5.278 | 1.560 | 4.882 | 1.483 | 5.548 | 1.720 | 6.870 | 1.879

For the multivariate benchmark, baselines include recently proposed time-
series/efficient Transformers – Autoformer, Informer, LogTrans, and Reformer
[16], and RNN variants – LSTnet [18], and LSTM [11]. Univariate baselines
further include N-BEATS [24], DeepAR [26], ARIMA, Prophet [30], and AutoETS
[3]. We obtain baseline results from the following papers: [37, 41], and
further run AutoETS from the Merlion library [3]. Table 1 summarize the
results of ETSformer against top performing baselines on a selection of
datasets, for the multivariate setting, and Table 4 in Appendix F for space.
Results for ETSformer are averaged over three runs (standard deviation in
Appendix G).

Overall, ETSformer achieves state-of-the-art performance, achieving the best
performance (across all datasets/settings, based on MSE) on 35 out of 40
settings for the multivariate case, and 17 out of 23 for the univariate case.
Notably, on Exchange, a dataset with no obvious periodic patterns, ETSformer
demonstrates an average (over forecast horizons) improvemnt of 39.8% over the
best performing baseline, evidencing its strong trend forecasting
capabilities. We highlight that for cases where ETSformer does not achieve the
best performance, it is still highly competitive, and is always within the top
2 performing methods, based on MSE, for 40 out of 40 settings in the
multivariate benchmark , and 21 out of 23 settings of the univariate case.

###  5.2 Ablation Study

Table 2: Ablation study on the various components of ETSformer, on the
horizon$=24$ setting.

Datasets | ETTh2 | ETTm2 | ECL | Traffic
---|---|---|---|---
ETSformer | MSE | 0.262 | 0.110 | 0.163 | 0.571
MAE | 0.337 | 0.222 | 0.287 | 0.373
w/o Level | MSE | 0.434 | 0.464 | 0.275 | 0.649
MAE | 0.466 | 0.518 | 0.373 | 0.393
w/o Season | MSE | 0.521 | 0.131 | 0.696 | 1.334
MAE | 0.450 | 0.236 | 0.677 | 0.779
w/o Growth | MSE | 0.290 | 0.115 | 0.167 | 0.583
MAE | 0.359 | 0.226 | 0.288 | 0.383
MH-ESA $\to$ MHA | MSE | 0.656 | 0.343 | 0.205 | 0.586
MAE | 0.639 | 0.451 | 0.323 | 0.380

We study the contribution of each major component which the final forecast is
composed of level, growth, and seasonality. Table 2 first presents the
performance of the full model, and subsequently, the performance of the
resulting model by removing each component. We observe that the composition of
level, growth, and season provides the most accurate forecasts across a
variety of application areas, and removing any one component results in a
deterioration. In particular, estimation of the level of the time-series is
critical. We also analyse the case where MH-ESA is replaced with a vanilla
multi-head attention, and observe that our trend attention formulation indeed
is more effective.

Figure 4:  Left: Visualization of decomposed forecasts from ETSformer and
Autoformer on a synthetic datset. (i) Ground truth and non-decomposed
forecasts of ETSformer and Autoformer on synthetic data. (ii) Trend component.
(iii) Seasonal component. The data sample on which Autoformer obtained lowest
MSE was selected for visualization. Right: Visualization of decomposed
forecasts from ETSformer on real world datasets, ETTh1, ECL, and Weather. Note
that season is zero-centered, and trend successfully tracks the level of the
time-series. Due to the long sequence forecasting setting and with a damping,
the growth component is not visually obvious, but notice for the Weather
dataset, the trend pattern is has a strong downward slope initially (near time
step 0), and is quickly damped.

###  5.3 Interpretability

ETSformer generates forecasts based on a composition of interpretable time-
series components. This means we can visualize each component individually,
and understand how seasonality and trend affects the forecasts. We showcase
this ability in Figure 4 on both synthetic and real world data. Experiments
with synthetic data are crucial in this case, since we are not able to obtain
the ground truth decomposition from real world data. ETSformer is first
trained on the synthetic dataset (details in Appendix E) with clear
(nonlinear) trend and seasonality patterns which we can control. Given a
lookback window (without noise), we visualize the forecast, as well as
decomposed trend and seasonal forecasts. ETSformer successfully forecasts
interpretable level, trend (level + growth), and seasonal components, as
observed in the trend and seasonality components closely tracking the ground
truth patterns. Despite obtaining a low MSE, the competing decomposition based
approach, Autoformer, struggles to disambiguate between trend and seasonality.

##  6 Conclusion

Inspired by the classical exponential smoothing methods and emerging
Transformer approaches for time-series forecasting, we proposed ETSformer, a
novel Transformer-based architecture for time-series forecasting which learns
level, growth, and seasonal latent representations and their complex
dependencies. ETSformer leverages the novel Exponential Smoothing Attention
and Frequency Attention mechanisms which are more effective at modeling time-
series than vanilla self-attention mechanism, and at the same time achieves
$\mathcal{O}(L\log L)$ complexity, where $L$ is the length of lookback window.
Our extensive empirical evaluation shows that ETSformer achieves state-of-the-
art performance, beating competing baselines in 35 out of 40 and 17 out of 23
settings for multivariate and univariate forecasting respectively. Future
directions include including additional covariates such as holiday indicators
and other dummy variables to consider holiday effects which cannot be captured
by the FA mechanism.

## References

  * Assimakopoulos and Nikolopoulos [2000] Vassilis Assimakopoulos and Konstantinos Nikolopoulos.  The theta model: a decomposition approach to forecasting.  _International journal of forecasting_ , 16(4):521–530, 2000.
  * Ba et al. [2016] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton.  Layer normalization, 2016.
  * Bhatnagar et al. [2021] Aadyot Bhatnagar, Paul Kassianik, Chenghao Liu, Tian Lan, Wenzhuo Yang, Rowan Cassius, Doyen Sahoo, Devansh Arpit, Sri Subramanian, Gerald Woo, et al.  Merlion: A machine learning library for time series.  _arXiv preprint arXiv:2109.09265_ , 2021.
  * Carion et al. [2020] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with transformers.  In _European Conference on Computer Vision_ , pages 213–229. Springer, 2020.
  * Cleveland et al. [1990] Robert B Cleveland, William S Cleveland, Jean E McRae, and Irma Terpenning.  Stl: A seasonal-trend decomposition.  _J. Off. Stat_ , 6(1):3–73, 1990.
  * Cleveland and Tiao [1976] William P Cleveland and George C Tiao.  Decomposition of seasonal time series: a model for the census x-11 program.  _Journal of the American statistical Association_ , 71(355):581–587, 1976.
  * De Livera et al. [2011] Alysha M De Livera, Rob J Hyndman, and Ralph D Snyder.  Forecasting time series with complex seasonal patterns using exponential smoothing.  _Journal of the American statistical association_ , 106(496):1513–1527, 2011.
  * Devlin et al. [2019] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _ArXiv_ , abs/1810.04805, 2019.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=YicbFdNTTy.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Holt [2004] Charles C Holt.  Forecasting seasonals and trends by exponentially weighted moving averages.  _International journal of forecasting_ , 20(1):5–10, 2004.
  * Hyndman et al. [2008] Rob Hyndman, Anne B Koehler, J Keith Ord, and Ralph D Snyder.  _Forecasting with exponential smoothing: the state space approach_.  Springer Science & Business Media, 2008.
  * Hyndman and Khandakar [2008] Rob J Hyndman and Yeasmin Khandakar.  Automatic time series forecasting: the forecast package for r.  _Journal of statistical software_ , 27(1):1–22, 2008.
  * Kingma and Ba [2015] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1412.6980.
  * Kitaev et al. [2020] Nikita Kitaev, Lukasz Kaiser, and Anselm Levskaya.  Reformer: The efficient transformer.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=rkgNKkHtvB.
  * Kuznetsov and Mohri [2015] Vitaly Kuznetsov and Mehryar Mohri.  Learning theory and algorithms for forecasting non-stationary time series.  In _NIPS_ , pages 541–549. Citeseer, 2015.
  * Lai et al. [2018] Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu.  Modeling long-and short-term temporal patterns with deep neural networks.  In _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_, pages 95–104, 2018.
  * Lee-Thorp et al. [2021] James Lee-Thorp, Joshua Ainslie, Ilya Eckstein, and Santiago Ontanon.  Fnet: Mixing tokens with fourier transforms.  _arXiv preprint arXiv:2105.03824_ , 2021.
  * Li et al. [2019] Shiyang Li, Xiaoyong Jin, Yao Xuan, Xiyou Zhou, Wenhu Chen, Yu-Xiang Wang, and Xifeng Yan.  Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting.  _ArXiv_ , abs/1907.00235, 2019.
  * Mathieu et al. [2014] Michaël Mathieu, Mikael Henaff, and Yann LeCun.  Fast training of convolutional networks through ffts.  _CoRR_ , abs/1312.5851, 2014.
  * McKenzie and Gardner Jr [2010] Eddie McKenzie and Everette S Gardner Jr.  Damped trend exponential smoothing: a modelling viewpoint.  _International Journal of Forecasting_ , 26(4):661–665, 2010.
  * Mohri and Rostamizadeh [2010] Mehryar Mohri and Afshin Rostamizadeh.  Stability bounds for stationary $\varphi$-mixing and $\beta$-mixing processes.  _Journal of Machine Learning Research_ , 11(2), 2010.
  * Oreshkin et al. [2019] Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio.  N-beats: Neural basis expansion analysis for interpretable time series forecasting.  In _International Conference on Learning Representations_ , 2019.
  * Raganato et al. [2020] Alessandro Raganato, Yves Scherrer, and Jörg Tiedemann.  Fixed encoder self-attention patterns in transformer-based machine translation.  In Trevor Cohn, Yulan He, and Yang Liu, editors, _Findings of the Association for Computational Linguistics: EMNLP 2020, Online Event, 16-20 November 2020_ , volume EMNLP 2020 of _Findings of ACL_ , pages 556–568. Association for Computational Linguistics, 2020.  doi: 10.18653/v1/2020.findings-emnlp.49.  URL https://doi.org/10.18653/v1/2020.findings-emnlp.49.
  * Salinas et al. [2020] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 36(3):1181–1191, 2020.  ISSN 0169-2070.  doi: https://doi.org/10.1016/j.ijforecast.2019.07.001.  URL https://www.sciencedirect.com/science/article/pii/S0169207019301888.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(56):1929–1958, 2014.  URL http://jmlr.org/papers/v15/srivastava14a.html.
  * Svetunkov [2016] Ivan Svetunkov.  _Complex exponential smoothing_.  Lancaster University (United Kingdom), 2016.
  * Tay et al. [2021] Yi Tay, Dara Bahri, Donald Metzler, Da-Cheng Juan, Zhe Zhao, and Che Zheng.  Synthesizer: Rethinking self-attention for transformer models.  In _International Conference on Machine Learning_ , pages 10183–10192. PMLR, 2021.
  * Taylor and Letham [2018] Sean J Taylor and Benjamin Letham.  Forecasting at scale.  _The American Statistician_ , 72(1):37–45, 2018\.
  * Theodosiou [2011] Marina Theodosiou.  Forecasting monthly and quarterly time series using stl decomposition.  _International Journal of Forecasting_ , 27(4):1178–1195, 2011.
  * Tsai et al. [2019] Yao-Hung Hubert Tsai, Shaojie Bai, Makoto Yamada, Louis-Philippe Morency, and Ruslan Salakhutdinov.  Transformer dissection: An unified understanding for transformer’s attention via the lens of kernel.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 4344–4353, Hong Kong, China, November 2019. Association for Computational Linguistics.  doi: 10.18653/v1/D19-1443.  URL https://aclanthology.org/D19-1443.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In _Advances in neural information processing systems_ , pages 5998–6008, 2017.
  * Vlachos et al. [2005] Michail Vlachos, Philip Yu, and Vittorio Castelli.  On periodicity detection and structural periodic similarity.  In _Proceedings of the 2005 SIAM international conference on data mining_ , pages 449–460. SIAM, 2005.
  * Winters [1960] Peter R Winters.  Forecasting sales by exponentially weighted moving averages.  _Management science_ , 6(3):324–342, 1960.
  * Woo et al. [2022] Gerald Woo, Chenghao Liu, Doyen Sahoo, Akshat Kumar, and Steven Hoi.  CoST: Contrastive learning of disentangled seasonal-trend representations for time series forecasting.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=PilZY3omXV2.
  * Wu et al. [2021] Haixu Wu, Jiehui Xu, Jianmin Wang, and Mingsheng Long.  Autoformer: Decomposition transformers with Auto-Correlation for long-term series forecasting.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Wu et al. [2020] Sifan Wu, Xi Xiao, Qianggang Ding, Peilin Zhao, Ying Wei, and Junzhou Huang.  Adversarial sparse transformer for time series forecasting.  In _NeurIPS_ , 2020.
  * You et al. [2020] Weiqiu You, Simeng Sun, and Mohit Iyyer.  Hard-coded gaussian attention for neural machine translation.  _arXiv preprint arXiv:2005.00742_ , 2020.
  * Zerveas et al. [2021] George Zerveas, Srideepika Jayaraman, Dhaval Patel, Anuradha Bhamidipaty, and Carsten Eickhoff.  A transformer-based framework for multivariate time series representation learning.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pages 2114–2124, 2021.
  * Zhou et al. [2021] Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang.  Informer: Beyond efficient transformer for long sequence time-series forecasting.  In _Proceedings of AAAI_ , 2021.

##  Appendix A Exponential Smoothing Attention

###  A.1 Exponential Smoothing Attention Matrix

$\displaystyle{\bm{A}}_{\mathrm{ES}}=\begin{bmatrix}(1-\alpha)^{1}&\alpha&0&0&\ldots&0\\\
(1-\alpha)^{2}&\alpha(1-\alpha)&\alpha&0&\ldots&0\\\
(1-\alpha)^{3}&\alpha(1-\alpha)^{2}&\alpha(1-\alpha)&\alpha&\ldots&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
(1-\alpha)^{L}&\alpha(1-\alpha)^{L-1}&\ldots&\alpha(1-\alpha)^{j}&\ldots&\alpha\end{bmatrix}$

###  A.2 Efficient Exponential Smoothing Attention Algorithm

Algorithm 1 PyTorch-style pseudocode of efficient $\mathcal{A}_{\mathrm{ES}}$

conv1d_fft: efficient convolution operation implemented with fast Fourier
transform (Appendix A, Algorithm 3),
outer: outer product

⬇

# V: value matrix, shape: L x d

# v0: initial state, shape: d

# alpha: smoothing parameter, shape: 1

# obtain exponentially decaying weights

# and compute weighted combination

powers = arange(L) # L

weight = alpha * (1 - alpha) ** flip(powers) # L

output = conv1d_fft(V, weight, dim=0) # L x d

# compute contribution from initial state

init_weight = (1 - alpha) ** (powers + 1) # L

init_output = outer(init_weight, v0) # L x d

return init_output + output

###  A.3 Level Smoothing via Exponential Smoothing Attention

$\displaystyle{\bm{E}}_{t}^{(n)}$
$\displaystyle=\bm{\alpha}*({\bm{E}}_{t}^{(n-1)}-{\bm{S}}_{t}^{(n)})+(1-\bm{\alpha})*({\bm{E}}_{t-1}^{(n)}+{\bm{B}}_{t-1}^{(n)})$
$\displaystyle=\bm{\alpha}*({\bm{E}}_{t}^{(n-1)}-{\bm{S}}_{t}^{(n)})+(1-\bm{\alpha})*{\bm{B}}_{t-1}^{(n)}$
$\displaystyle\hskip
30.00005pt+(1-\bm{\alpha})*[\bm{\alpha}*({\bm{E}}_{t-1}^{(n-1)}-{\bm{S}}_{t-1}^{(n)})+(1-\bm{\alpha})*({\bm{E}}_{t-2}^{(n)}+{\bm{B}}_{t-2}^{(n)})]$
$\displaystyle=\bm{\alpha}*({\bm{E}}_{t}^{(n-1)}-{\bm{S}}_{t}^{(n)})+\bm{\alpha}*(1-\bm{\alpha})*({\bm{E}}_{t-1}^{(n-1)}-{\bm{S}}_{t-1}^{(n)})$
$\displaystyle\hskip
30.00005pt+(1-\bm{\alpha})*{\bm{B}}_{t-1}^{(n)}+(1-\bm{\alpha})^{2}*{\bm{B}}_{t-2}^{(n)}$
$\displaystyle\hskip
30.00005pt+(1-\bm{\alpha})^{2}[\bm{\alpha}*({\bm{E}}_{t-2}^{(n-1)}-{\bm{S}}_{t-2}^{(n)})+(1-\bm{\alpha})*({\bm{E}}_{t-3}^{(n)}+{\bm{B}}_{t-3}^{(n)})]$
$\displaystyle\mathmakebox[\widthof{{}={}}][c]{\vdots}$
$\displaystyle=(1-\bm{\alpha})^{t}({\bm{E}}_{0}^{(n)}-{\bm{S}}_{0}^{(n)})+\sum_{j=0}^{t-1}\bm{\alpha}*(1-\bm{\alpha})^{j}*({\bm{E}}_{t-j}^{(n-1)}-{\bm{S}}_{t-j}^{(n)})+\sum_{k=1}^{t}(1-\bm{\alpha})^{k}*{\bm{B}}_{t-k}^{(n)}$
$\displaystyle=\mathcal{A}_{\mathrm{ES}}({\bm{E}}_{t-L:t}^{(n-1)}-{\bm{S}}_{t-L:t}^{(n)})+\sum_{k=1}^{t}(1-\bm{\alpha})^{k}*{\bm{B}}_{t-k}^{(n)}$

Based on the above expansion of the level equation, we observe that
${\bm{E}}_{n}^{(t)}$ can be computed by a sum of two terms, the first of which
is given by an $\mathcal{A}_{\mathrm{ES}}$ term, and we finally, we note that
the second term can also be calculated using the conv1d_fft algorithm,
resulting in a fast implementation of level smoothing.

###  A.4 Further Details on ESA Implementation

Algorithm 2 PyTorch-style pseudocode of naive $\mathcal{A}_{\mathrm{ES}}$

mm: matrix multiplication, outer: outer product
repeat: einops style tensor operations,
gather: gathers values along an axis specified by dim

⬇

# V: value matrix, shape: L x d

# v0: initial state, shape: d

# alpha: smoothing parameter, shape: 1

L, d = V.shape

# obtain exponentially decaying weights

powers = arange(L) # L

weight = alpha * (1 - alpha).pow(flip(powers)) # L

# perform a strided roll operation

# rolls a matrix along the columns in a strided manner

# i.e. first row is shifted right by L-1 positions,

# second row is shifted L-2, …, last row is shifted by 0.

weight = repeat(weight, 'L -> T L', T=L) # L x L

indices = repeat(arange(L), 'L -> T L', T=L)

indices = (indices - (arange(L) + 1).unsqueeze(1)) % L

weight = gather(weight, dim=-1, index=indices)

# triangle masking to achieve the exponential smoothing attention matrix

weight = triangle_causal_mask(weight)

output = mm(weight, V)

init_weight = (1 - alpha) ** (powers + 1)

init_output = outer(init_weight, v0)

return init_output + output

Algorithm 3 PyTorch-style pseudocode of conv1d_fft

next_fast_len: find the next fast size of input data to fft, for zero-padding,
etc.
rfft: compute the one-dimensional discrete Fourier Transform for real input
x.conj(): return the complex conjugate, element-wise
irfft: computes the inverse of rfft
roll: roll array elements along a given axis
index_select: returns a new tensor which index es the input tensor along
dimension dim using the entries in index

⬇

# V: value matrix, shape: L x d

# weight: exponential smoothing attention vector, shape: L

# dim: dimension to perform convolution on

# obtain lengths of sequence to perform convolution on

N = V.size(dim)

M = weight.size(dim)

# Fourier transform on inputs

fast_len = next_fast_len(N + M - 1)

F_V = rfft(V, fast_len, dim=dim)

F_weight = rfft(weight, fast_len, dim=dim)

# multiplication and inverse

F_V_weight = F_V * F_weight.conj()

out = irfft(F_V_weight, fast_len, dim=dim)

out = out.roll(-1, dim=dim)

# select the correct indices

idx = range(fast_len - N, fast_len)

out = out.index_select(dim, idx)

return out

Algorithm 2 describes the naive implementation for ESA by first constructing
the exponential smoothing attention matrix, ${\bm{A}}_{\mathrm{ES}}$, and
performing the full matrix-vector multiplication. Efficient
$\mathcal{A}_{\mathrm{ES}}$ relies on Algorithm 3, to achieve an
$\mathcal{O}(L\log L)$ complexity, by speeding up the matrix-vector
multiplication. Due to the structure lower triangular structure of
${\bm{A}}_{\mathrm{ES}}$ (ignoring the first column), we note that performing
a matrix-vector multiplication with it is equivalent to performing a
convolution with the last row. Algorithm 3 describes the pseudocode for fast
convolutions using fast Fourier transforms.

##  Appendix B Discrete Fourier Transform

The DFT of a sequence with regular intervals,
${\bm{x}}=({x}_{0},{x}_{1},\ldots,{x}_{N-1})$ is a sequence of complex
numbers,

${c}_{k}=\sum_{n=0}^{N-1}{x}_{n}\cdot\exp(-i2\pi kn/N),$

for $k=0,1,\ldots,N-1$, where $c_{k}$ are known as the Fourier coefficients of
their respective Fourier frequencies. Due to the conjugate symmetry of DFT for
real-valued signals, we simply consider the first $\lfloor N/2\rfloor+1$
Fourier coefficients and thus we denote the DFT as
$\mathcal{F}:\mathbb{R}^{N}\to{\mathbb{C}}^{\lfloor N/2\rfloor+1}$. The DFT
maps a signal to the frequency domain, where each Fourier coefficient can be
uniquely represented by the amplitude, $|{c}_{k}|$, and the phase,
$\phi({c}_{k})$,

$\displaystyle|{c}_{k}|$
$\displaystyle=\sqrt{\mathfrak{R}\\{{c}_{k}\\}^{2}+\mathfrak{I}\\{{c}_{k}\\}^{2}}$
$\displaystyle\phi({c}_{k})$
$\displaystyle=\tan^{-1}\bigg{(}\frac{\mathfrak{I}\\{{c}_{k}\\}}{\mathfrak{R}\\{{c}_{k}\\}}\bigg{)}$

where $\mathfrak{R}\\{{c}_{k}\\}$ and $\mathfrak{I}\\{{c}_{k}\\}$ are the real
and imaginary components of ${c}_{k}$ respectively. Finally, the inverse DFT
maps the frequency domain representation back to the time domain,

${x}_{n}=\mathcal{F}^{-1}({\bm{c}})_{n}=\frac{1}{N}\sum_{k=0}^{N-1}c_{k}\cdot\exp(i2\pi
kn/N),$

##  Appendix C Implementation Details

###  C.1 Hyperparameters

For all experiments, we use the same hyperparameters for the encoder layers,
decoder stacks, model dimensions, feedforward layer dimensions, number of
heads in multi-head exponential smoothing attention, and kernel size for input
embedding as listed in Table 3. We perform hyperparameter tuning via a grid
search over the number of frequencies $K$, lookback window size, and learning
rate, selecting the settings which perform the best on the validation set
based on MSE (on results averaged over three runs). The search range is
reported in Table 3, where the lookback window size search range was decided
to be set as the values for the horizon sizes for the respective datasets.

Table 3: Hyperparameters used in ETSformer. Hyperparameter | Value
---|---
Encoder layers | 2
Decoder stacks | 2
Model dimension | 512
Feedforward dimension | 2048
Multi-head ESA heads | 8
Input embedding kernel size | 3
K | $K\in\\{0,1,2,3\\}$
Lookback window size | $L\in\\{96,192,336,720\\}$
Lookback window size (ILI) | $L\in\\{24,36,48,60\\}$
Learning rate | $lr\in\\{1\mathrm{e}{-3},3\mathrm{e}{-4},1\mathrm{e}{-4},3\mathrm{e}{-5},1\mathrm{e}{-5}\\}$

###  C.2 Optimization

We use the Adam optimizer [15] with $\beta_{1}=0.9$, $\beta_{2}=0.999$, and
$\epsilon=1e-08$, and a batch size of 32. We schedule the learning rate with
linear warmup over 3 epochs, and cosine annealing thereafter for a total of 15
training epochs for all datasets. The minimum learning rate is set to 1e-30.
For smoothing and damping parameters, we set the learning rate to be 100 times
larger and do not use learning rate scheduling. Training was done on an Nvidia
A100 GPU.

###  C.3 Regularization

We apply two forms of regularization during the training phase.

##### Data Augmentations

We utilize a composition of three data augmentations, applied in the following
order - scale, shift, and jitter, activating with a probability of 0.5.

  1. 1.

Scale – The time-series is scaled by a single random scalar value, obtained by
sampling $\epsilon\sim\mathcal{N}(0,0.2)$, and each time step is
$\tilde{x}_{t}=\epsilon x_{t}$.

  2. 2.

Shift – The time-series is shifted by a single random scalar value, obtained
by sampling $\epsilon\sim\mathcal{N}(0,0.2)$ and each time step is
$\tilde{x}_{t}=x_{t}+\epsilon$.

  3. 3.

Jitter – I.I.D. Gaussian noise is added to each time step, from a distribution
$\epsilon_{t}\sim\mathcal{N}(0,0.2)$, where each time step is now
$\tilde{x}_{t}=x_{t}+\epsilon_{t}$.

##### Dropout

We apply dropout [27] with a rate of $p=0.2$ across the model. Dropout is
applied on the outputs of the Input Embedding, Frequency Self-Attention and
Multi-Head ES Attention blocks, in the Feedforward block (after activation and
before normalization), on the attention weights, as well as damping weights.

##  Appendix D Datasets

ETT222https://github.com/zhouhaoyi/ETDataset Electricity Transformer
Temperature [41] is a multivariate time-series dataset, comprising of load and
oil temperature data recorded every 15 minutes from electricity transformers.
ETT consists of two variants, ETTm and ETTh, whereby ETTh is the hourly-
aggregated version of ETTm, the original 15 minute level dataset.

ECL333lhttps://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014
Electricity Consuming Load measures the electricity consumption of 321
households clients over two years, the original dataset was collected at the
15 minute level, but is pre-processed into an hourly level dataset.

Exchange444https://github.com/laiguokun/multivariate-time-series-data Exchange
[18] tracks the daily exchange rates of eight countries (Australia, United
Kingdom, Canada, Switzerland, China, Japan, New Zealand, and Singapore) from
1990 to 2016.

Traffic555https://pems.dot.ca.gov/ Traffic is an hourly dataset from the
California Department of Transportation describing road occupancy rates in San
Francisco Bay area freeways.

Weather666https://www.bgc-jena.mpg.de/wetter/ Weather measures 21
meteorological indicators like air temperature, humidity, etc., every 10
minutes for the year of 2020.

ILI777https://gis.cdc.gov/grasp/fluview/fluportaldashboard.html Influenza-like
Illness records the ratio of patients seen with ILI and the total number of
patients on a weekly basis, obtained by the Centers for Disease Control and
Prevention of the United States between 2002 and 2021.

##  Appendix E Synthetic Dataset

The synthetic dataset is constructed by a combination of trend and seasonal
component. Each instance in the dataset has a lookack window length of 192 and
forecast horizon length of 48. The trend pattern follows a nonlinear,
saturating pattern, $b(t)=\frac{1}{1+\exp{\beta_{0}(t-\beta_{1})}}$, where
$\beta_{0}=-0.2,\beta_{1}=192$. The seasonal pattern follows a complex
periodic pattern formed by a sum of sinusoids. Concretely,
$s(t)=A_{1}\cos(2\pi f_{1}t)+A_{2}\cos(2\pi f_{2}t$, where
$f_{1}=1/10,f_{2}=1/13$ are the frequencies, $A_{1}=A_{2}=0.15$ are the
amplitudes. During training phase, we use an additional noise component by
adding i.i.d. gaussian noise with 0.05 standard deviation. Finally, the $i$-th
instance of the dataset is $x_{i}=[x_{i}(1),x_{i}(2),\ldots,x_{i}(192+48)]$,
where $x_{i}(t)=b(t)+s(t+i)+\epsilon$.

##  Appendix F Univariate Forecasting Benchmark

Table 4: Univariate forecasting results over various forecast horizons. Best
results are bolded, and second best results are underlined.

Methods | ETSformer | Autoformer | Informer | N-BEATS | DeepAR | Prophet | ARIMA | AutoETS
---|---|---|---|---|---|---|---|---
Metrics | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE
ETTm2 | 96 | 0.080 | 0.212 | 0.065 | 0.189 | 0.088 | 0.225 | 0.082 | 0.219 | 0.099 | 0.237 | 0.287 | 0.456 | 0.211 | 0.362 | 0.794 | 0.617
192 | 0.150 | 0.302 | 0.118 | 0.256 | 0.132 | 0.283 | 0.120 | 0.268 | 0.154 | 0.310 | 0.312 | 0.483 | 0.261 | 0.406 | 1.078 | 0.740
336 | 0.175 | 0.334 | 0.154 | 0.305 | 0.180 | 0.336 | 0.226 | 0.370 | 0.277 | 0.428 | 0.331 | 0.474 | 0.317 | 0.448 | 1.279 | 0.822
720 | 0.224 | 0.379 | 0.182 | 0.335 | 0.300 | 0.435 | 0.188 | 0.338 | 0.332 | 0.468 | 0.534 | 0.593 | 0.366 | 0.487 | 1.541 | 0.924
Exchange | 96 | 0.099 | 0.230 | 0.241 | 0.299 | 0.591 | 0.615 | 0.156 | 0.299 | 0.417 | 0.515 | 0.828 | 0.762 | 0.112 | 0.245 | 0.192 | 0.316
192 | 0.223 | 0.353 | 0.273 | 0.665 | 1.183 | 0.912 | 0.669 | 0.665 | 0.813 | 0.735 | 0.909 | 0.974 | 0.304 | 0.404 | 0.355 | 0.442
336 | 0.421 | 0.497 | 0.508 | 0.605 | 1.367 | 0.984 | 0.611 | 0.605 | 1.331 | 0.962 | 1.304 | 0.988 | 0.736 | 0.598 | 0.577 | 0.578
720 | 1.114 | 0.807 | 0.991 | 0.860 | 1.872 | 1.072 | 1.111 | 0.860 | 1.890 | 1.181 | 3.238 | 1.566 | 1.871 | 0.935 | 1.242 | 0.865

##  Appendix G ETSformer Standard Deviation

Table 5: ETSformer main benchmark results with standard deviation. Experiments
are performed over three runs.

(a) Multivariate benchmark. Metrics | MSE (SD) | MAE (SD)
---|---|---
ETTm2 | 96 | 0.189 (0.002) | 0.280 (0.001)
192 | 0.253 (0.002) | 0.319 (0.001)
336 | 0.314 (0.001) | 0.357 (0.001)
720 | 0.414 (0.000) | 0.413 (0.001)
ECL | 96 | 0.187 (0.001) | 0.304 (0.001)
192 | 0.199 (0.001) | 0.315 (0.002)
336 | 0.212 (0.001) | 0.329 (0.002)
720 | 0.233 (0.006) | 0.345 (0.006)
Exchange | 96 | 0.085 (0.000) | 0.204 (0.001)
192 | 0.182 (0.003) | 0.303 (0.002)
336 | 0.348 (0.004) | 0.428 (0.003)
720 | 1.025 (0.031) | 0.774 (0.014)
Traffic | 96 | 0.607 (0.005) | 0.392 (0.005)
192 | 0.621 (0.015) | 0.399 (0.013)
336 | 0.622 (0.003) | 0.396 (0.003)
720 | 0.632 (0.004) | 0.396 (0.004)
Weather | 96 | 0.197 (0.007) | 0.281 (0.008)
192 | 0.237 (0.005) | 0.312 (0.004)
336 | 0.298 (0.003) | 0.353 (0.003)
720 | 0.352 (0.007) | 0.388 (0.002)
ILI | 24 | 2.527 (0.061) | 1.020 (0.021)
36 | 2.615 (0.103) | 1.007 (0.013)
48 | 2.359 (0.056) | 0.972 (0.011)
60 | 2.487 (0.006) | 1.016 (0.007)

(b) Univariate benchmark. Metrics | MSE (SD) | MAE (SD)
---|---|---
ETTm2 | 96 | 0.080 (0.001) | 0.212 (0.001)
192 | 0.150 (0.024) | 0.302 (0.026)
336 | 0.175 (0.012) | 0.334 (0.014)
720 | 0.224 (0.008) | 0.379 (0.006)
Exchange | 96 | 0.099 (0.003) | 0.230 (0.003)
192 | 0.223 (0.015) | 0.353 (0.009)
336 | 0.421 (0.002) | 0.497 (0.000)
720 | 1.114 (0.049) | 0.807 (0.016)

##  Appendix H Computational Efficiency

(a) Runtime Efficiency Analysis

(b) Memory Efficiency Analysis

Figure 5: Computational Efficiency Analysis. Values reported are based on the
training phase of ETTm2 multivariate setting. Horizon is fixed to 48 for
lookback window plots, and lookback is fixed to 48 for forecast horizon plots.
For runtime efficiency, values refer to the time for one iteration. The “ "
marker indicates an out-of-memory error for those settings.

In this section, our goal is to compare ETSformer’s computational efficiency
with that of competing Transformer-based approaches. Visualized in Figure 5,
ETSformer maintains competitive efficiency with compting quasilinear
complexity Transformers, while obtaining state-of-the-art performance.
Furthermore, due to ETSformer’s unique decoder architecture which relies on
its Trend Damping and Frequency Attention modules rather than output
embeddings, ETSformer maintains superior efficiency as forecast horizon
increases.


††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

# Logarithmic Quantum Time Crystal

Haipeng Xue (薛海鹏)  International Center for Quantum Materials, School of
Physics, Peking University, Beijing 100871, China  Lingchii Kong (孔令琦)
International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China  Biao Wu (吴飙<EMAIL_ADDRESS>International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China  Wilczek Quantum Center, School of Physics
and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Collaborative Innovation Center of Quantum Matter, Beijing 100871, China

###### Abstract

We investigate a time-independent many-boson system, whose ground states are
quasi-degenerate and become infinitely degenerate in the thermodynamic limit.
Out of these quasi-degenerate ground states we construct a quantum state that
evolves in time with a period that is logarithmically proportional to the
number of particles, that is, $T\sim\log N$. This boson system in such a state
is a quantum time crystal as it approaches the ground state in the
thermodynamic limit. The logarithmic dependence of its period on the total
particle number $N$ makes it observable experimentally even for systems with
very large number of particles. Possible experimental proposals are discussed.

##  I INTRODUCTION

Spontaneous symmetry breaking is a well known physical phenomenon, where the
observed ground state of a many-particle system does not possess the
symmetries of its Hamiltonian[1], e.g., ferromagnetic materials break the
rotational symmetry and crystals break the spatial translational symmetry.
Wilczek suggested the possibility of spontaneous breaking of time
translational symmetry[2, 3]: the observed ground state of a closed quantum
system may oscillate periodically in time. This suggestion had drawn some
quick criticism[4, 5]. In 2015, a no-go theorem was proved to exclude the
possibility of spontaneous continuous time translational symmetry breaking in
the ground state for a wide class of Hamiltonians with short-range
interactions[6]. However, it was realized later that quantum time crystals can
exist in a periodically driven system[7]. This is now known as discrete time
crystal; its equilibrium state can vary with a period that is multiple of the
driving period[8, 9, 10]. Two experiments with trapped ions and nitrogen-
vacancy centers were performed, confirming the existence of discrete time
crystals[11, 12]. Recently, the discussion about time crystals has expanded to
the systems with long-range interaction[13].

When spontaneous symmetry breaking occurs in a time independent system, the
observed ground state is not the true ground state but a superposition of many
(quasi-)degenerate ground states[14]. Two energy states are quasi-degenerate
if the energy gap between them approaches zero in the thermodynamic limit.
There have some efforts to construct a time crystal with these
(quasi-)degenerate ground states[15, 16, 17]. For such a time crystal, it
typically oscillates with a period that grows polynomially with particle
number, i.e. $T\sim O(N^{n})$. It will be impossible in experiments to keep a
many-particle state from decoherence for such a long time.

In this work we study a time-independent many-boson system. This boson system
has many quasi-degenerate ground states, which become infinitely degenerate in
the thermodynamic limit. With these ground states, we are able to construct
states that oscillate with time with a period $T$ that is logarithmically
proportional to the number of particles $N$, that is, $T\sim\log N$. Due to
this logarithmic dependence, the period is short enough for possible
experimental observation even for very large $N$. We find that the observables
that are usually used in time crystal experiments are not suitable for our
systems. We suggest to observe the time crystal by measuring the square of
overlap between the states of such a quantum time crystal at different times,
i.e., $\left|\langle\Psi(0)|\Psi(t)\rangle\right|^{2}$ with a technique based
on the Hong–Ou–Mandel interference[18].

##  II  Two-mode interacting boson systems

The interacting system of $N$ bosons that we are going to consider has only
two modes. In a certain parameter range, this system has many quasi-ground
states, which approach the true ground state in the thermodynamics limit
$N\rightarrow\infty$. This is the critical feature of this system and a
necessary condition for spontaneous symmetry breaking to occur.

###  II.1 Theoretical model

The two-mode interacting boson system is described by the following
Hamiltonian[19, 20, 14]

$\displaystyle\hat{H}$
$\displaystyle=-(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}_{1})+\frac{\gamma}{N}[\hat{n}_{1}(\hat{n}_{1}-1)+\hat{n}_{2}(\hat{n}_{2}-1)]$
$\displaystyle\quad+\frac{4\gamma}{N}\hat{n}_{1}\hat{n}_{2}+\frac{\gamma}{N}(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{1}\hat{a}_{2}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{2}\hat{a}_{1}\hat{a}_{1})\,,$
(1)

where $\hat{a}^{\dagger}_{1}(\hat{a}_{1})$ and
$\hat{a}^{\dagger}_{2}(\hat{a}_{2})$ are the creation (annihilation) operators
of two modes and $N$ is the total number of bosons. This simple theoretical
model can now be realized in experiments[21]. In the above, we have set the
strength of single particle hopping as the unit of energy and therefore
$\gamma$ is a dimensionless parameter characterizing the interaction strength
and the pair hopping. In addition, for simplicity, we set $\hbar\equiv 1$.

In the thermodynamic limit $N\rightarrow\infty$, this quantum many-body system
can be described by a mean-field model. Mathematically, one can simply replace
the operators $\hat{a}^{\dagger}_{1,2}$ and $\hat{a}_{1,2}$ with their
corresponding complex variables $a^{*}_{1,2}$ and $a_{1,2}$ in Eq. (II.1)[14].
We choose a different set of variables $P=(|a_{2}|^{2}-|a_{1}|^{2})/2$ and
$Q=\arg(a_{2})-\arg(a_{1})$. Physically, $P$ and $Q$ are the particle number
difference and the phase difference between the two modes, respectively. In
terms of $P$ and $Q$, the mean-field model of the system is

$H_{m}=-\sqrt{1-4P^{2}}\cos Q+\gamma\left(1-4P^{2}\right)\cos^{2}Q\,.$ (2)

Note that $P$ and $Q$ are a pair of dynamical variables that are canonically
conjugate to each other. Fig. 1 is the energy landscape in the phase space of
this mean-field Hamiltonian. When $\gamma<1/2$, the center of the phase space
$P=Q=0$ has the lowest energy. When $\gamma>1/2$, all the points on the dashed
line in Fig. 1(b) have the same lowest energy. That means that the system is
infinitely degenerate when $\gamma>1/2$. As we shall see in the next
subsection, this infinite degeneracy corresponds to a set of quasi-degenerate
ground states in the quantum model (II.1).

Figure 1: The energy landscape in the phase space of the mean-field model
$H_{m}$. (a) When $\gamma<1/2$, the system has only one ground state at
$(Q,P)=(0,0)$; (b) when $\gamma>1/2$, the system has infinitely degenerate
ground states, which are located on the dashed line.

###  II.2 Quantum energy levels

For mathematical simplicity, we transform this boson system into a
mathematically equivalent spin system. With
$\hat{n}_{1,2}=\hat{a}^{\dagger}_{1,2}\hat{a}_{1,2}$, we introduce
$\hat{S}_{x}=(\hat{a}^{\dagger}_{1}\hat{a}_{2}+\hat{a}^{\dagger}_{2}\hat{a}_{1})/2,\quad\hat{S}_{y}=(\hat{a}^{\dagger}_{1}\hat{a}_{2}-\hat{a}^{\dagger}_{2}\hat{a}_{1})/2i,\quad\hat{S}_{z}=(\hat{n}_{1}-\hat{n}_{2})/2$.
The quantum model (II.1) becomes

$\hat{H}^{\prime}=\hat{H}/2=-\hat{S}_{x}+\frac{2\gamma}{N}\hat{S}^{2}_{x}\,.$
(3)

This Hamiltonian is also called Lipkin-Meshkov-Glick model and describes a
spin system with $S=N/2$ [22]. In this spin formalism the energy eigenstates
and energy levels are rather obvious. As $\hat{S}_{x}$ commutes with
$\hat{H}$, the eigenstates of the system are the eigenstates of $\hat{S}_{x}$,

$\hat{S}_{x}|m\rangle=m|m\rangle\,,~{}~{}m=-S,-S+1,\cdots,0,\cdots,S\,.$ (4)

For eigenstate $\ket{m}$, its corresponding energy level is

$E_{m}=\frac{2\gamma}{N}m^{2}-m\,.$ (5)

When $\gamma<1/2$, the lowest energy level is at $m=S=N/2$ and the energy
level increases monotonically as $m$ decreases. And the energy gap $\delta
E_{m}$ between two neighboring levels $\ket{m-1}$ and $\ket{m}$ is

$\delta E_{m}=1-\frac{2\gamma}{N}(2m-1)>1-2\gamma\,.$ (6)

This means that the energy gap remains finite even in the thermodynamic limit
$N\rightarrow\infty$.

Figure 2: The energy levels for $\gamma>1/2$ which are labeled by quantum
number $m$ and denoted by black bars. The ground state is at $m_{0}$. The
energy gap $\delta E$ between neighboring levels around $m_{0}$ is $\propto
1/N$. All the energy levels marked by the gray area are quasi-ground states,
whose energies become degenerate when $N\rightarrow\infty$. The width of the
gray area is about $\sim\sqrt{N}$. The energy gaps outside the gray area
remain finite even at the large $N$ limit.

When $\gamma>1/2$, the energy level is the lowest at $m_{0}$, which is the
largest integer smaller than or equal to $N/(4\gamma)$. The highest energy
level is at $m=-S$. Without loss of essential physics, we choose
$m_{0}=N/(4\gamma)$ by choosing an appropriate value of $\gamma$. The energy
gap between the ground state and the first excited state is $2\gamma/N$; the
gap between the highest energy level and the second highest is
$1+2\gamma-2\gamma/N$. The former approaches zero and the latter remains
finite at the limit $N\rightarrow\infty$. This means that for $\gamma>1/2$ the
energy gap $\delta E_{m}$ between two neighboring levels has different
dependence on $N$ at different levels $m$. We find that for a set of energy
levels near the ground state the gap approaches zero and for others it remains
finite at $N\rightarrow\infty$. This is shown schematically in Fig. 2. To
accurately to describe this behavior, we define the energy gap $\Delta E_{m}$
between energy level $\ket{m}$ and the ground state $\ket{m_{0}}$ and we find
that

$\Delta E_{m}=\frac{2\gamma}{N}(m-m_{0})^{2}\,.$ (7)

This means that any energy level $\ket{m}$ will approach the ground state
energy in the thermodynamic limit $N\rightarrow\infty$ when $m$ satisfies

$|m-m_{0}|\sim O\left(N^{1/2-\delta}\right)\,,$ (8)

where $0<\delta<1/2$. These energy levels are located in the shadow area in
Fig. 2. We call them quasi-ground states, which form a sub-Hilbert space $C$.
At the limit $N\rightarrow\infty$, all these quasi-ground states become the
mean-field ground states on the black dashed line in Fig. 1(b).

Figure 3: Energy eigenstates $\ket{m}$ in the quantum phase space for $N=440$
and $\gamma=3/4$. (a) The true ground state; (b) the $29$th exited state. The
black lines represent the mean-field infinite degenerate ground states.

###  II.3 Eigenstates in quantum phase space

As discussed above, in the case of $\gamma>1/2$ there are many quasi-ground
states in the quantum model $\hat{H}$, which correspond to the infinitely many
ground states in the mean-field model. Such a correspondence becomes more
apparent and insightful when we plot these quasi-ground states in the quantum
phase space (see Fig. 3).

The quantum phase space is obtained by dividing the mean-field phase space in
Fig. 1 into Planck cells. The size of a Planck cell, which is denoted by
$(Q_{i},P_{j})$, is $1/N$. A Wannier function
$\left\\{\ket{Q_{i},P_{j}}\right\\}$, which is localized along both $P$ and
$Q$ directions, is assigned for each Planck cell and all the Wannier functions
form a complete set of orthonormal basis[23, 24, 25, 26, 27]. An energy
eigenstate $\ket{m}$ is then projected to the quantum phase space as

$\ket{m}=\sum_{i,j}\ket{Q_{i},P_{j}}\braket{Q_{i},P_{j}}{m}\,.$ (9)

What is plotted in Fig. 3 is $|\langle Q,P|m\rangle|^{2}$. Two examples are
shown in Fig. 3: one is the true ground state $\ket{m_{0}}$ and the other is a
quasi-ground state. It is clear from the figure that both eigenstates are
concentrated around the mean-field ground states, which are marked by black
solid lines in Fig. 3.

##  III  LOGARITHMIC TIME CRYSTAL

The quantum time crystal can be constructed as follows. We choose an
eigenstate $\ket{m_{1}}$ with

$m_{1}=m_{0}-\Big{\lfloor}\sqrt{\frac{N}{\log N}}\Big{\rfloor}\,,$ (10)

where $\lfloor x\rfloor$ is the largest integer $\leq x$. This state
$\ket{m_{1}}$ certainly belongs to the sub-Hilbert space $C$ as $\Delta
E_{m_{1}}=2\gamma/\log N$ approaches zero when $N\rightarrow\infty$. We
construct a superposition state
$\ket{\Phi_{0}}=(\ket{m_{0}}+\ket{m_{1}})/\sqrt{2}$. It will evolve with time
as

$|\Phi(t)\rangle=\frac{e^{-iE_{m_{0}}t}}{\sqrt{2}}\left[|m_{0}\rangle+\exp\left(-\frac{i2\pi
t}{T_{0}\log{N}}\right)|m_{1}\rangle\right]\,,$ (11)

where $T_{0}=\pi/\gamma$ is a constant with order $O(1)$ and we have ignored
the minor difference between $x$ and $\lfloor x\rfloor$. It is clear that this
quantum state will oscillate with a period of $T_{0}\log(N)$.

In real experiments, it is hard to prepare a state that is a superposition of
only two eigenstates. We consider a double-Gaussian superposition state
$|\Psi(0)\rangle=\sum_{m}c_{m}|m\rangle$, where

$c_{m}=\frac{1}{\sqrt{2\pi\sigma}}\Big{[}e^{-\frac{(m-m_{0})^{2}}{4\sigma^{2}}}+e^{-\frac{(m-m_{1})^{2}}{4\sigma^{2}}}\Big{]}\,,$
(12)

where $\sigma$ is the width of the Gaussian distribution. Since we do not want
the Gaussian distributions to spread out the whole sub-Hilbert space $C$ of
quasi-ground states, the width $\sigma$ should be much smaller than
$|m_{1}-m_{0}|$. For simplicity, we choose $\sigma\sim O(1)$.

The overlap between $|\Psi(t)\rangle$ and $|\Psi(0)\rangle$ can be evaluated
as

$\displaystyle\langle\Psi(0)|\Psi(t)\rangle=\sum_{m=-N/2}^{N/2}|c_{m}|^{2}\exp\left(-i\Delta
E_{m}t\right)\,.$ (13)

When $N$ is very large, we can replace the summation by integral,

$\displaystyle\langle\Psi(0)|\Psi(t)\rangle\approx\int_{-\infty}^{\infty}dm\;|c_{m}|^{2}\exp\left(-i\Delta
E_{m}t\right)$ $\displaystyle=$ $\displaystyle
A(t)\left\\{1+\exp\left[-\frac{N^{2}}{4\Sigma^{2}(t)\log^{2}N}\right]\exp\left(-i\frac{2\pi
t}{T(t)}\right)\right\\}\,,$ (14)

where

$\displaystyle T(t)$ $\displaystyle=T_{0}\log
N\left(1+\frac{64\gamma^{2}\sigma^{4}t^{2}}{N^{2}}\right),$
$\displaystyle\Sigma(t)$ $\displaystyle=\sqrt{\left(\frac{N}{8\gamma\sigma
t}\right)^{2}+{\sigma^{2}}},$ $\displaystyle A(t)$
$\displaystyle=\frac{1}{2}\sqrt{\frac{N}{{N}+8i{\gamma\sigma^{2}t}}}.$ (15)

It is clear from $T(t)$ that when $t\ll O(N)$ the system oscillates with a
period of $T_{0}\log N$. $A(t)$ shows that the amplitude of the oscillation
will gradually decrease with time. However, within the time interval $0<t\ll
O(N)$ we should be able to observe many periods of oscillations since the
period of oscillation is of order $O(\log(N))$. In addition, the oscillations
are enveloped by a Gaussian function with width $\Sigma(t)$ which increases
with time. Instead of doing integral, one can also do the summation in Eq.(13)
numerically. Both results are shown in Fig. 4 for comparison; they agree very
well.

Figure 4: The oscillations of logarithmic time crystal (12). The solid line is
the analytical result in Eq.(III) and the squares are numerical results
obtained by directly summing Eq.(13). $N=29240$, $m_{1}=m_{0}+91$, $\sigma=1$.

##  IV EXPERIMENTAL Scheme

In literature, the time crystal dynamics is usually illustrated with the
correlation function of the polarization along the $z$ direction,
$\langle\hat{S}_{z}(0)\hat{S}_{z}(t)\rangle$ [28]. For the quantum time
crystal Eq.(11), we have

$\displaystyle\langle\hat{S}_{z}(0)\hat{S}_{z}(t)\rangle$
$\displaystyle=A_{0,+}e^{-i\delta E_{m_{0}+1}t}+A_{0,-}e^{-i\delta
E_{m_{0}}t}$ $\displaystyle\quad+A_{1,+}e^{-i\delta
E_{m_{1}+1}t}+A_{1,-}e^{-i\delta E_{m_{1}}t}\,,$ (16)

where $A_{j,\pm}=(N/2+1)N/2-m_{j}(m_{j}\pm 1)$, $j=0,1$. The four
characteristic frequencies $\delta E_{m_{j}\pm 1}$ are just the energy gaps
between the neighboring levels at the two states $|m_{0}\rangle$ and
$|m_{1}\rangle$. Therefore, we have $\delta E_{m_{j}\pm 1}\ll O(1/\log N)$ and
the oscillation period is much longer than $\log N$. So the correlation
function Eq.(IV) can not be used to observe the dynamic behavior of our
logarithmic quantum time crystal. We propose to use the experimental scheme
shown in Fig. 5, where $\left|\langle\Psi(0)|\Psi(t)\rangle\right|^{2}$ is
measured.

Figure 5: The experiment scheme. (i) Two identical boson systems $a$ and $b$;
each has two modes, denoted by $1,2$. Initially, the two systems are
independent and they are prepared in states, $|\Psi(0)\rangle$ and
$|\Psi(t)\rangle$, respectively. (ii) We allow the atom tunneling between the
modes $a_{i}$ and $b_{i}$ ($i=1,2$) to realize Hong-Ou-Mandel interference.
(iii) We count particle number parity in system $b$. Repeat the above process
for many times. This scheme effectively measures
$\left|\langle\Psi(0)|\Psi(t)\rangle\right|^{2}$.

In Fig. 5, there are two identical two-mode boson systems, which are denoted
as $a$ and $b$, created with atoms and optical devices. The tunneling between
system $a$ and $b$ can be tuned: initially the optical potential barrier
between them is high so that there is no tunneling; the barrier is lowered in
the later stage of the experiment to allow Hong-Ou-Mandel interference.
Initially when the tunneling is off, both system $a$ and $b$ have the same
number of bosons and evolve with time under the Hamiltonian Eq. II.1 .

A swap operator $\hat{V}$ between these two systems is defined as

$\hat{V}|\psi\rangle_{a}\otimes|\phi\rangle_{b}=|\phi\rangle_{a}\otimes|\psi\rangle_{b}\,,$
(17)

where $|\cdot\rangle_{a}$ and $|\cdot\rangle_{b}$ are states of systems $a$
and $b$. It is easy to see that
$\langle\hat{V}\rangle=\left|\langle\psi|\phi\rangle\right|^{2}$, i.e., the
expectation value of the swap operator $\hat{V}$ is the overlap between two
quantum states $\ket{\psi}$ and $\ket{\phi}$. As $\hat{V}$ has only two
eigenvalues $\pm 1$, measuring the expectation of operator $\hat{V}$ requires
statistically counting its eigenvalues within a large ensemble. In our system,
$\hat{V}=\hat{V}_{1}\otimes\hat{V}_{2}$, where $\hat{V}_{i}$ is the swapping
on individual modes $i=1,2$. So measuring $\hat{V}$ can be achieved by
measuring $\hat{V}_{1}$ and $\hat{V}_{2}$.

The symmetric eigenstate $|S\rangle$ and asymmetric eigenstate $|AS\rangle$ of
$\hat{V}_{i}$ can be expressed in terms of Fock states as follows

$\displaystyle|S\rangle$
$\displaystyle=\text{superpose}\left\\{(\hat{a}_{i}^{\dagger}-\hat{b}_{i}^{\dagger})^{2\alpha}(\hat{a}_{i}^{\dagger}+\hat{b}_{i}^{\dagger})^{\beta}|\text{vac}\rangle\right\\},$
(18) $\displaystyle|AS\rangle$
$\displaystyle=\text{superpose}\left\\{(\hat{a}_{i}^{\dagger}-\hat{b}_{i}^{\dagger})^{2\alpha-1}(\hat{a}_{i}^{\dagger}+\hat{b}_{i}^{\dagger})^{\beta}|\text{vac}\rangle\right\\}\,,$
(19)

where $\hat{a}_{i}^{\dagger}(\hat{b}_{i}^{\dagger})$ refers to the creation
operator on mode $a_{i}(b_{i})$ and $\alpha$ and $\beta$ are positive
integers. We perform the following transformation,

$\left(\hat{a}_{i}^{\dagger}+\hat{b}_{i}^{\dagger}\right)/\sqrt{2}\rightarrow\hat{a}_{i}^{\dagger}\,,~{}~{}~{}\left(\hat{a}_{i}^{\dagger}-\hat{b}_{i}^{\dagger}\right)/\sqrt{2}\rightarrow\hat{b}_{i}^{\dagger}\,.$
(20)

This transformation can be realized experimentally with atoms tunneling
between $a_{i}$ and $b_{i}$ under a Bose-Hubbard Hamiltonian [29]. This is
essentially Hong-Ou-Mandel interference. After the transformation, the
eigenstates of $\hat{V}_{i}$ become

$\displaystyle|S^{\prime}\rangle$
$\displaystyle=\text{superpose}\left\\{(\hat{b}_{i}^{\dagger})^{2\alpha}(\hat{a}_{i}^{\dagger})^{\beta}|\text{vac}\rangle\right\\},$
$\displaystyle|AS^{\prime}\rangle$
$\displaystyle=\text{superpose}\left\\{(\hat{b}_{i}^{\dagger})^{2\alpha-1}(\hat{a}_{i}^{\dagger})^{\beta}|\text{vac}\rangle\right\\}.$
(21)

This means that observing even(odd) particle numbers in mode $b_{i}$
corresponds to eigenvalues $+1(-1)$ of $\hat{V}_{i}$. Hence if the particle
numbers in mode $b_{1}$ and $b_{2}$ have the same parity, we know the measured
quantity of $\hat{V}$ is $+1$; otherwise, it is $-1$.

The basic experimental steps are shown in Fig. 5. Initially, we prepare system
$b$ in state $|\Psi(0)\rangle$ and let it evolve for a period of time $t$ so
the state becomes $|\Psi(t)\rangle$. At this same moment, we prepare system
$a$ in state $|\Psi(0)\rangle$. Then we turn off the tunneling between
different modes and turn on the the tunneling between system $a$ and $b$,
realizing Hong-Ou-Mandel interference. Finally, we observe the number of
particles in both modes of system $b$. By repeat this experiment many times
for different $t$,
$\langle\hat{V}\rangle=\left|\langle\Psi(0)|\Psi(t)\rangle\right|^{2}$ is
obtained.

##  V  Discussion And Conclusion

Before concluding, we need to mention that performing perturbed external
fields [16] actually fails to produce a logarithmic time crystal state.
Because the field, like $g\hat{S}_{z}$, can only superpose the state
$|m_{0}\rangle$ with its nearest levels. To be honest, how to realize
logarithmic time crystal is still an open question for us. But we do have some
speculations. In solid physics, electrons can be excited from the band bottom
to the band top by the phonons with specific momentum and energy, where the
momentum of phonons is tuned by the structure of Brillouin-zone. In our model,
the system can be excited by a laser from $|m_{0}\rangle$ to $|m_{1}\rangle$
when the frequency of laser is equal to $\Delta
E_{m_{1}}=E_{m_{1}}-E_{m_{0}}$. But only controlling the frequency cannot
assure that the intermediate states are not superposed. So, we need to
speculate that there may be some extra characteristics (like Brillouin-zone
with respect to electrons) about our logarithmic time crystal state to rule
the momentum of laser so that it can be achieved.

In conclusion, we have studied a two-mode boson system of $N$ particles. In a
certain parameter regime, this system has many quasi-ground states, which
forms a sub-Hilbert space $C$ and become degenerate with the true ground state
in the thermodynamic limit $N\rightarrow\infty$. We have been able to
superpose these quasi-ground states and construct a quantum time crystal that
oscillates with a period that is logarithmically proportional to the number of
particles, $T\sim O(\log(N))$. Such a logarithmic dependence makes the period
short enough for possible experimental observation even when the number of
particle is large. A experimental scheme based on Hong-Ou-Mandel interference
has been proposed.

###### Acknowledgements.

We thank Yijie Wang for helpful discussion. This work is supported by the
National Key R&D Program of China (Grants No. 2017YFA0303302, No.
2018YFA0305602), National Natural Science Foundation of China (Grant No.
11921005), and Shanghai Municipal Science and Technology Major Project (Grant
No.2019SHZDZX01).

## References

  * Beekman  _et al._ [2019] A. J. Beekman, L. Rademaker, and J. van Wezel, An introduction to spontaneous symmetry breaking, SciPost Physics Lecture Notes , 11 (2019).
  * Shapere and Wilczek [2012] A. Shapere and F. Wilczek, Classical time crystals, Physical Review Letters 109, 160402 (2012).
  * Wilczek [2012] F. Wilczek, Quantum time crystals, Physical Review Letters 109, 160401 (2012).
  * Nozières [2013] P. Nozières, Time crystals: Can diamagnetic currents drive a charge density wave into rotation?, EPL (Europhysics Letters) 103, 57008 (2013).
  * Bruno [2013] P. Bruno, Impossibility of spontaneously rotating time crystals: a no-go theorem, Physical Review Letters 111, 070402 (2013).
  * Watanabe and Oshikawa [2015] H. Watanabe and M. Oshikawa, Absence of quantum time crystals, Physical Review Letters 114, 251603 (2015).
  * Sacha [2015] K. Sacha, Modeling spontaneous breaking of time-translation symmetry, Physical Review A 91, 033617 (2015).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Physical Review Letters 117, 090402 (2016).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Physical Review Letters 116, 250401 (2016).
  * Yao  _et al._ [2017] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations, Physical Review Letters 118, 030401 (2017).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, _et al._ , Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017).
  * Zhang  _et al._ [2017] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, _et al._ , Observation of a discrete time crystal, Nature 543, 217 (2017).
  * Kozin and Kyriienko [2019] V. K. Kozin and O. Kyriienko, Quantum time crystals from hamiltonians with long-range interactions, Physical Review Letters 123, 210602 (2019).
  * Zhu  _et al._ [2015] Q. Zhu, Q. Zhang, and B. Wu, Extended two-site bose–hubbard model with pair tunneling: spontaneous symmetry breaking, effective ground state and fragmentation, Journal of Physics B: Atomic, Molecular and Optical Physics 48, 045301 (2015).
  * Shapere and Wilczek [2019] A. D. Shapere and F. Wilczek, Regularizations of time-crystal dynamics, Proceedings of the National Academy of Sciences 116, 18772 (2019).
  * Huang  _et al._ [2018] Y. Huang, T. Li, and Z. Yin, Symmetry-breaking dynamics of the finite-size lipkin-meshkov-glick model near ground state, Physical Review A 97, 012115 (2018).
  * Syrwid  _et al._ [2021] A. Syrwid, A. Kosior, and K. Sacha, Can a bright soliton model reveal a genuine time crystal for a finite number of bosons?, EPL (Europhysics Letters) 134, 66001 (2021).
  * Hong  _et al._ [1987] C.-K. Hong, Z.-Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Physical Review Letters 59, 2044 (1987).
  * Liang  _et al._ [2009] J.-Q. Liang, J.-L. Liu, W.-D. Li, and Z.-J. Li, Atom-pair tunneling and quantum phase transition in the strong-interaction regime, Physical Review A 79, 033617 (2009).
  * Cao and Fu [2012] H. Cao and L. Fu, Quantum phase transition and dynamics induced by atom-pair tunnelling of bose-einstein condensates in a double-well potential, The European Physical Journal D 66, 1 (2012).
  * Hemmerich [2019] A. Hemmerich, Bosons condensed in two modes with flavor-changing interaction, Physical Review A 99, 013623 (2019).
  * Lipkin  _et al._ [1965] H. J. Lipkin, N. Meshkov, and A. Glick, Validity of many-body approximation methods for a solvable model:(i). exact solutions and perturbation theory, Nuclear Physics 62, 188 (1965).
  * Neumann [1929] J. v. Neumann, Beweis des ergodensatzes und des h-theorems in der neuen mechanik, Zeitschrift für Physik 57, 30 (1929).
  * Han and Wu [2015] X. Han and B. Wu, Entropy for quantum pure states and quantum h theorem, Physical Review E 91, 062106 (2015).
  * Jiang  _et al._ [2017] J. Jiang, Y. Chen, and B. Wu, Quantum ergodicity and mixing and their classical limits with quantum kicked rotor, arXiv preprint arXiv:1712.04533 (2017).
  * Fang  _et al._ [2018] Y. Fang, F. Wu, and B. Wu, Quantum phase space with a basis of wannier functions, Journal of Statistical Mechanics: Theory and Experiment 2018, 023113 (2018).
  * Wang  _et al._ [2021] Z. Wang, J. Feng, and B. Wu, Microscope for quantum dynamics with planck cell resolution, Physical Review Research 3, 033239 (2021).
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).