# Automated Knee X-ray Report Generation

Aydan Gasimova
Biomedical Image Analysis Group
Imperial College London
&Giovanni Montana
Imaging and Biomedical Engineering
King’s College London
&Daniel Rueckert
Biomedical Image Analysis Group
Imperial College London

###### Abstract

Gathering manually annotated images for the purpose of training a predictive
model is far more challenging in the medical domain than for natural images as
it requires the expertise of qualified radiologists. We therefore propose to
take advantage of past radiological exams (specifically, knee X-ray
examinations) and formulate a framework capable of learning the correspondence
between the images and reports, and hence be capable of generating diagnostic
reports for a given X-ray examination consisting of an arbitrary number of
image views. We demonstrate how aggregating the image features of individual
exams and using them as conditional inputs when training a language generation
model results in auto-generated exam reports that correlate well with
radiologist-generated reports.

##  1 Introduction

The burden on diagnostic radiologists has been increasing quite rapidly with
the advancement of a variety of modern imaging modalities, which, while
allowing for higher resolution images in 3D and even 4D, dramatically increase
the complexity of the diagnostic process. It has become common for
radiologists to rely on various image analysis and automated decision support
systems to facilitate the interpretation process. These computer aided
diagnostic systems, or CAD, have been able to make great advances with the
help of machine learning algorithms and a large amount of clinical imaging
data, and have been successfully tested in clinical settings [1, 2, 3, 4].

Many of these machine learning algorithms are supervised learning approaches
that require large amounts of annotated image data for training. Gathering
suitable data is especially challenging in the medical domain as it is
incredibly time consuming for radiologists to generate ground-truths of the
standard and volume required for training a predictive model. An alternative
approach is to use past clinical images and corresponding radiological reports
available through a hospital’s picture archiving and communication system
(PACS); the advantage being that, although this data is largely unstructured
(free text), it is available to us in high volumes and removes the need for
manual annotation.

Hence, it is necessary to develop an appropriate learning framework: one that
can take advantage of past radiological examinations and their corresponding
textual reports, and formulate a method for prediction that can best expedite
the diagnostic process. An additional motivation to learning from clinical
reports is that reports provide a much richer context to the pathology than a
single class label, such as severity and location.

Here we present a model for captioning knee X-ray exams that is capable of
generating radiological reports summarising present pathologies and the
contexts within which they are presented. A key contribution is that the
proposed model is able to handle an arbitrary number of input image views
(which is typical in the context of knee X-ray exams) and is able to learn
from free-text reports. The model is based on the neural image caption (NIC)
model of [5], where a word-level language generation model is conditioned on
image features. In our implementation, the X-ray image features for each input
image are derived from a pre-trained classification network, and aggregated
into a fixed-sized input to the language model.

##  2 Related Work

In the field of computer vision, the use of human generated visually
descriptive text to infer the contents of an image has primarily been applied
to image caption generation. Earlier models relied on linking template-based
language models to objects and spatial contexts in the image [6, 7, 8]. More
recently, interest has moved to the combined potential of convolutional neural
networks (CNNs) and recurrent neural networks (RNNs) for describing images
using natural language [9, 5, 10, 11, 12]. The advantage of using neural
networks for caption generation is that the model is not constricted by hard-
coded language templates and is able to learn more freely from the training
data.

By contrast, using clinical reports to perform automated pathology detection
in medical images has only been explored in a few studies [13, 14, 15, 16]. In
addition, these models are limited to learning from carefully curated reports,
either medical subject heading annotations (MeSH® [17]) [13, 14], or templated
reports made specifically for the purpose of training [15, 16], which are both
very time consuming for radiologists/pathologists to create.

To our knowledge, this is the first study exploring the use of past medical
examinations and their corresponding raw reports in order to develop a
predictive model for an arbitrary number of input radiological images.

##  3 Dataset

The knee X-ray dataset has been extracted from the PAC system of St Thomas
Hospital (part of Guys and St Thomas NHS Foundation Trust) and has been fully
anonymised to remove sensitive patient information. It comprises a total of
330 knee X-ray exams collected over the years 2015 and 2016. Each exam
consists of a textual report and one or more X-ray images (left/right knees or
both, taken from different views: anteroposterior (AP), lateral (L) and
skyline (S), and different positions: weight-bearing (WB) and non-weight-
bearing(nonWB)). The most common exam consists of both AP and L views of left
and right knees separately, making up 42% of total exams. The reports vary in
length between 2 and 145 words, with an average of 30 and standard deviation
of 18.7; and between 1 and 16 sentences with an average of 2.7 per report. The
X-ray images vary in sizes between $420\times 650\times 3$ and $3056\times
3056\times 3$.

Prior to using the raw reports for training an image captioning model, they
were cleaned to remove phrases not directly describing the images, for
instance, instructions to compare to past exams (which are not available to us
due to the anonymisation process). These phrases do not allow any inferences
about current images, so are not useful for finding image-text correspondence.
In addition, common ‘stopwords’ (such as ‘and’, ‘the’) were removed, and all
punctuation was replaced with either commas or full stops where appropriate.
All report cleaning was done using regular expressions. Exams were randomly
split into 80%-10%-10% for training, testing and validation respectively.

##  4 Medical Image Report Generation Model

Figure 1: Medical image report generation network.

###  4.1 Image Modelling

We adopt the GoogLeNet [18] CNN architecture, pre-trained on the ImageNet
dataset [19], and extract image features of each X-ray image view ($V_{1}$ \-
$V_{K}$) from the last spatial average pooling layer ($\mathbb{R}^{1024}$).
The maximum value of each feature is aggregated across the exam images to
create a fixed-size input to the RNN of dimension $\mathbb{R}^{1024}$, which
is then passed through a fully connected layer in order to reduce the
dimension to $\mathbb{R}^{256}$, equal to the RNN input size.

###  4.2 Language Modelling

We use a recurrent neural network, specifically the Long Short-Term Memory
(LSTM) implementation proposed in [20], to model the report sequence. LSTM
models are widely used in machine translation [21, 22, 23] and natural image
and video captioning [9, 24, 12] due to their ability to capture long-term
dependencies, and to reduce the problem of vanishing gradients in vanilla
RNNs. Each LSTM unit has three sigmoid gates to control the internal state:
‘input’, ‘output’ and ‘forget’. At each time step, the gates control how much
of the previous time steps is propagated through to determine the output. For
an input word sequence $\\{x_{1},\dots,x_{n}\\}$, the internal hidden state
$h_{t}\in\mathbb{R}^{256}$ and memory state $m_{t}\in\mathbb{R}^{256}$ are
updated as follows:

$h_{t}=f_{t}\odot h_{t-1}+i_{t}\odot\tanh(W^{(hx)}x_{t}+W^{(hm)}m_{t-1})$
$m_{t}=o_{t}\odot\tanh(h_{t})$

where $x_{t}\in\mathbb{R}^{D}$ is the word embedding, $W^{(hx)}$ and
$W^{(hm)}$ are the trainable weight parameters, and $i_{t}$, $o_{t}$ and
$f_{t}$ are the input, output and forget gates respectively.

In order for the sequence generation model to be conditioned on the set of
input images, the combined images features for each exam are input at time
step $t=0$, and the words in the report input at consequent time steps. The
complete network, including a bounding box (BBox) regressor to extract the
knee joint, is illustrated in Figure 1.

##  5 Training

For training the report generation model, the LSTM is unrolled up to 33 time
steps (1 for image features, 1 per start and end token, and 30 for the average
number of words in the reports). The images were rescaled and padded
(preserving the aspect ratio) to $224\times 224$. Words were one-hot encoded,
and ones that occurred at a frequency <5 were removed.

To improve the X-ray image features, a simple CNN BBox regressor was built and
trained on a subset of the training images and ground-truth BBoxes manually
created for this purpose (231 image-BBox pairs). In addition, the training set
was augmented eight-fold by random cropping the original images from
$256\times 256$ to $224\times 224$, (in the case of the BBox images, small,
random translations (between -5 and 5 in x and y) were applied to the BBox of
each image), flipping the images along the vertical axis, and shuffling the
sentences in the reports.

We then trained the LSTM model for report generation conditioned on the
combined image features by minimising the negative log-likelihood between the
output and true sequence:

$L(S,I)=-\sum_{t=0}^{N}\log p(P_{t}=T_{t}|\text{CNN($I$)},P_{0}\ldots
P_{t-1})$

where $p$ is the probability that the predicted word $P_{t}$ equals the true
word $T_{t}$ at time step $t$ given image features CNN($I$) and previous words
$P_{0}\ldots P_{t-1}$, and $N$ is the LSTM sequence length. At training time,
loss is minimised over the training set using stochastic gradient descent
(batch size 20, learning rate $1\times 10^{-5}$), and parameters are updated
using Adam [25] optimisation.

##  6 Evaluation

Reports are generated on the test set by max-aggregating the image features of
an exam, using it as the first word input into the LSTM, and sampling
consequent words. The quality of the generated reports was evaluated by
measuring BLUE [26] and METEOR [27] scores averaged over all the reports,
which are a modified form of n-gram precision commonly used for evaluating
image captioning and machine translation as they maintain high correlation
with human judgement. As a comparison, we trained the model on single image
inputs, duplicating the reports for each image in an exam. We also evaluated
whether the performance improved if we narrow the input to the CNN to a BBox
surrounding the knee joint.

The report generation model performed better when trained on a combination of
image features than on single images as the conditional input (see Table 1),
suggesting that all images are required, in part, in order to make the
pathology assessment. This is in line with expectations since pathologies may
only be present in one leg, or may only be seen from a particular view.
Contrary to expectation, bounding box extraction did not improve the
performance of report generation as it may have instead removed pertinent
image features. A sample exam with the original report and two auto-generated
reports of contrasting BLEU-1 scores is presented in Figure 2.

| BLEU-1 | BLEU-2 | BLEU-3 | BLEU-4 | METEOR
---|---|---|---|---|---
| tr | te | tr | te | tr | te | tr | te | tr | te
Baseline, single image input | 42.2 | 33.2 | 13.3 | 5.7 | 3.8 | 1.9 | 1.3 | 1.1 | 26.7 | 22.2
Max-aggr. of image features | 60.7 | 40.4 | 32.6 | 10.1 | 19.4 | 2.6 | 12.3 | 1.2 | 41.1 | 35.7
Max-aggr.+BBoxes | 38.9 | 37.4 | 11.3 | 7.1 | 3.4 | 1.1 | 1.2 | 0.2 | 28.3 | 28.9

Table 1: BLEU n-gram and METEOR scores evaluated on report generation model
trained on single image inputs (Baseline), max-aggregate of image features,
and max-aggregate of image features extracted from BBox detections. Evaluated
on training (tr) and test (te) data.

---

Original: Joint spaces articular surfaces appear preserved. Significant
degenerative erosive change seen. Good prediction (B1=87.5): Joint spaces
articular surfaces appear preserved bilaterally. Poor prediction (B1=28.6):
Joint space narrowing medial compartments bilaterally.
---

Figure 2: Sample knee X-ray exam with corresponding original report (with
common stopwords removed) and two sampled predictions generated by the model
trained on max-aggregates of image features.

##  7 Conclusion

We demonstrate how past knee X-ray exams consisting of sets of images and
corresponding radiological reports can be used as part of a learning framework
to determine image-text correspondence and hence automate the generation of
such reports for new X-ray exams. Preliminary results look promising as the
auto-generated reports correlate well with true reports, and we hope to train
the model on additional knee X-ray exams as these become available to us.
Further developments to the model can be made by incorporating the knowledge
of the view-type of each image, keeping them as separate inputs, and finding
correspondence between each image and parts of text (for instance, through the
use of an attention mechanism [28]). We will explore this in a future study.

## References

  * Freer and Ulissey [2001] Timothy W Freer and Michael J Ulissey.  Screening mammography with computer-aided detection: Prospective study of 12,860 patients in a community breast center 1.  _Radiology_ , 220(3):781–786, 2001.
  * Ramírez et al. [2013] Javier Ramírez, JM Górriz, Diego Salas-Gonzalez, A Romero, Míriam López, Ignacio Álvarez, and Manuel Gómez-Río.  Computer-aided diagnosis of alzheimer’s type dementia combining support vector machines and discriminant set of features.  _Information Sciences_ , 237:59–72, 2013.
  * Litjens et al. [2015] Geert J S Litjens, Jelle O. Barentsz, Nico Karssemeijer, and Henkjan J. Huisman.  Clinical evaluation of a computer-aided diagnosis system for determining cancer aggressiveness in prostate MRI.  _European Radiology_ , 25(11):3187–3199, 2015\.
  * Kobayashi et al. [2016] Hajime Kobayashi, Masaki Ohkubo, Akihiro Narita, Janaka C Marasinghe, Kohei Murao, Toru Matsumoto, Shusuke Sone, and Shinichi Wada.  A method for evaluating the performance of computer-aided detection of pulmonary nodules in lung cancer ct screening: detection limit for nodule size and density.  _The British Journal of Radiology_ , 90(0):20160313, 2016.
  * Vinyals et al. [2015] Oriol Vinyals, Alexander Toshev, Samy Bengio, and Dumitru Erhan.  Show and tell: A neural image caption generator.  In _CVPR_ , pages 3156–3164, 2015.
  * Farhadi et al. [2010] Ali Farhadi, Mohsen Hejrati, Mohammad Amin Sadeghi, Peter Young, Cyrus Rashtchian, Julia Hockenmaier, and David Forsyth.  Every picture tells a story: Generating sentences from images.  In _ECCV_ , pages 15–29. Springer, 2010.
  * Kulkarni et al. [2013] Girish Kulkarni, Visruth Premraj, Vicente Ordonez, Sagnik Dhar, Siming Li, Yejin Choi, Alexander C Berg, and Tamara L Berg.  Babytalk: Understanding and generating simple image descriptions.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 35(12):2891–2903, 2013.
  * Fidler et al. [2013] Sanja Fidler, Abhishek Sharma, and Raquel Urtasun.  A sentence is worth a thousand pixels.  In _CVPR_ , pages 1995–2002, 2013.
  * Kiros et al. [2014] Ryan Kiros, Ruslan Salakhutdinov, and Richard S Zemel.  Unifying visual-semantic embeddings with multimodal neural language models.  _arXiv preprint arXiv:1411.2539_ , 2014.
  * Karpathy and Li [2015] Andrej Karpathy and Fei-fei Li.  Deep visual-semantic alignments for generating image descriptions.  In _CVPR_ , pages 3128–3137, 2015.
  * Chen and Lawrence Zitnick [2015] Xinlei Chen and C Lawrence Zitnick.  Mind’s eye: A recurrent visual representation for image caption generation.  In _CVPR_ , pages 2422–2431, 2015.
  * Xu et al. [2015a] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron C Courville, Ruslan Salakhutdinov, Richard S Zemel, and Yoshua Bengio.  Show, attend and tell: Neural image caption generation with visual attention.  In _ICML_ , volume 14, pages 77–81, 2015a.
  * Shin et al. [2016a] Hoo-Chang Shin, Le Lu, Lauren Kim, Ari Seff, Jianhua Yao, and Ronald M Summers.  Interleaved Text / Image Deep Mining on a Large-Scale Radiology Database for Automated Image Interpretation.  _Journal of Machine Learning Research (JMLR)_ , 17:1–31, 2016a.  doi: 10.1109/CVPR.2015.7298712.
  * Shin et al. [2016b] Hoo-chang Shin, Kirk Roberts, Le Lu, Dina Demner-Fushman, Jianhua Yao, and Ronald M Summers.  Learning to read chest x-rays: Recurrent neural cascade model for automated image annotation.  In _CVPR_ , pages 2497–2506, 2016b.
  * Zhang et al. [2017a] Zizhao Zhang, Yuanpu Xie, Fuyong Xing, Mason Mcgough, and Lin Yang.  Mdnet: A semantically and visually interpretable medical image diagnosis network.  In _CVPR_ , July 2017a.
  * Zhang et al. [2017b] Zizhao Zhang, Pingjun Chen, Manish Sapkota, and Lin Yang.  Tandemnet: Distilling knowledge from medical images using diagnostic reports as optional semantic references.  In _MICCAI_ , 2017b.
  * Lipscomb [2000] Carolyn E Lipscomb.  Medical subject headings (mesh).  _Bulletin of the Medical Library Association_ , 88(3):265, 2000.
  * Szegedy et al. [2015] Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, and Andrew Rabinovich.  Going deeper with convolutions.  In _CVPR_ , pages 1–9, 2015.
  * Deng et al. [2009] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _CVPR_ , pages 248–255. IEEE, 2009.
  * Hochreiter and Schmidhuber [1997] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Sutskever et al. [2014] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  In _NIPS_ , pages 3104–3112, 2014.
  * Cho et al. [2014] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  _arXiv preprint arXiv:1406.1078_ , 2014.
  * Luong et al. [2015] Minh-Thang Luong, Hieu Pham, and Christopher D Manning.  Effective approaches to attention-based neural machine translation.  _arXiv preprint arXiv:1508.04025_ , 2015.
  * Donahue et al. [2015] Jeffrey Donahue, Lisa Anne Hendricks, Sergio Guadarrama, Marcus Rohrbach, Subhashini Venugopalan, Kate Saenko, and Trevor Darrell.  Long-term recurrent convolutional networks for visual recognition and description.  In _CVPR_ , pages 2625–2634, 2015.
  * Kingma and Ba [2014] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Papineni et al. [2002] Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th annual meeting on association for computational linguistics_ , pages 311–318. Association for Computational Linguistics, 2002.
  * Banerjee and Lavie [2005] Satanjeev Banerjee and Alon Lavie.  Meteor: An automatic metric for mt evaluation with improved correlation with human judgments.  In _Proceedings of the acl workshop on intrinsic and extrinsic evaluation measures for machine translation and/or summarization_ , volume 29, pages 65–72, 2005.
  * Xu et al. [2015b] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhudinov, Rich Zemel, and Yoshua Bengio.  Show, attend and tell: Neural image caption generation with visual attention.  In _ICML_ , pages 2048–2057, 2015b.


# A Decade of Time Crystals: Quo Vadis?

Peter Hannaford  Optical Sciences Centre, Swinburne University of Technology,
Hawthorn, Victoria 3122, Australia  Krzysztof Sacha  Instytut Fizyki
Teoretycznej, Uniwersytet Jagielloński, ulica Profesora Stanisława
Łojasiewicza 11, PL-30-348 Kraków, Poland

###### Abstract

Ten years ago, the new era of time crystals began. Time crystals are systems
that behave in the time dimension like ordinary space crystals do in space
dimensions. We present a brief history of a decade of research on time
crystals, describe current research directions, indicate challenges, and
discuss some future perspectives for condensed matter physics in the time
domain.

Ten years of history of time crystals. – The new era of time crystals began in
2012 when Frank Wilczek alone and together with Alfred Shapere published two
seminal papers on quantum and classical time crystals Wilczek (2012); Shapere
and Wilczek (2012). It is very well known that interacting quantum particles
can self-organize their distribution in space and spontaneously form space
crystals. However, prior to 2012 no one raised the question of whether similar
self-organization could also happen in time. Wilczek’s idea of quantum time
crystals was ambitious because he proposed that such a crystal could
spontaneously form even if a quantum many-body system were in the ground
state. Soon after, it turned out that for many-body systems with quite general
two-body interactions this was not feasible Bruno (2013); Watanabe and
Oshikawa (2015); Syrwid et al. (2017).

The classical version of time crystals proposed by Shapere and Wilczek relied
on a single-particle system with kinetic energy that possessed a minimum at
non-zero velocity Shapere and Wilczek (2012). At the lowest possible energy,
the particle performed periodic motion, and thus revealed crystalline
structure in time. In the Hamilton description, this seemed impossible because
a minimum of the classical Hamiltonian implied vanishing velocity. Shapere and
Wilczek solved this puzzle by showing that at a minimum, the Hamiltonian had a
cusp where the Hamilton equations were not well defined Shapere and Wilczek
(2012).

Let us return to quantum time crystals. The original Wilczek idea did not work
but three years later there was a resurrection of quantum time crystals in
closed but periodically driven systems (i.e. systems that are described by
unitary evolution, and thus there is no coupling to the environment, but the
energy is not a constant of motion because the Hamiltonian depends explicitly
on time). In 2015 it was shown that ultra-cold atoms bouncing on an
oscillating atom mirror could spontaneously self-reorganize their motion and
start evolving with a period twice longer than the period dictated by the
drive Sacha (2015a) — a new periodic motion appeared, i.e., a novel
crystalline structure emerged in the time domain. The following year, the same
phenomena were proposed in spin-based solid state systems Khemani et al.
(2016); Else et al. (2016). The spins were being periodically flipped and if
there were interactions between them, the systems evolved with a period twice
as long as the driving period even when the spin flipping process was
imperfect. Periodically driven systems possess discrete time translation
symmetry which corresponds to a translation in time by the period of the
external drive. The symmetry is spontaneously broken when time crystals form
and therefore they have been dubbed discrete (or Floquet) time crystals Yao et
al. (2017).

The first experiments on discrete time crystals were performed in 2017. Two
groups implemented periodically driven spin systems and demonstrated that the
presence of the spin-spin interactions resulted in periodic evolution of the
systems with a period twice or three times as long as the period of the drive
Zhang et al. (2017); Choi et al. (2017). Period-doubling time crystals were
later realized in other spin-based experiments Pal et al. (2018); Rovny et al.
(2018); Kyprianidis et al. (2021); Randall et al. (2021); Mi et al. (2022);
Frey and Rachel (2022) and in ultra-cold atoms Smits et al. (2018).

Discrete time crystals that evolve with a period twice as long as the period
of the drive are small crystals in the time dimension. However, it was shown
that big discrete time crystals evolving with a period even 100 times longer
than the driving period are also attainable in the laboratory Giergiel et al.
(2018a); Surace et al. (2019); Pizzi et al. (2019); Giergiel et al. (2020).
These systems consist of many elementary cells in time and allow for
realization of condensed matter phenomena in the time domain like in the case
of ordinary space crystals, which consist of many elementary cells in space
Hannaford and Sacha (2020). It was also shown that spontaneous breaking of
discrete time translation symmetry can lead to the formation of not only
crystal structures in time, but also quasi-crystal structures in time Li et
al. (2012); Giergiel et al. (2019a).

Apart from the described discrete time crystals in closed systems, other
research directions of the field of time crystals have been developed, which
we list in the following.

In 2019 it turned out that in time-independent systems if, instead of two-body
interactions, multi-body interactions were considered, it was possible to
construct a system where time crystal behavior could emerge in the ground
state Kozin and Kyriienko (2019). How to realize such a system experimentally,
which would be a realization of Wilczek’s original idea of a quantum time
crystal, is not known yet.

All time crystals described so far correspond to closed systems. However, it
is known from the systematic studies of Prigogine and others that open systems
that are out of thermodynamic equilibrium can possess steady states that
spontaneously break the time translation symmetry and reveal periodic
evolution in time Nicolis and Prigogine (1977); Prigogine (1978). In the
context of time crystals, dissipative time-independent and periodically driven
systems have recently been investigated theoretically and experimentally Gong
et al. (2018); Iemini et al. (2018); Buča et al. (2019); Dogra et al. (2019);
Cosme et al. (2019); Keßler et al. (2021); Taheri et al. (2022). They are
important from the point of view of practical applications because in the real
world, i.e., outside advanced laboratories, it can be difficult to avoid
contact with the environment.

If bosons form a Bose-Einstein condensate and the system is described in the
framework of the grand canonical ensemble, then the number of particles
fluctuates and only the average particle number can be well defined. Then,
Bose–Einstein condensation is identified with breaking of the U(1) symmetry
and the phase of the order parameter oscillates periodically in time,
indicating that the time translation symmetry is spontaneously broken. The
oscillations can be observed because the system is not isolated and can
exchange particles with another system that plays the role of a reference
frame with respect to which the oscillations can be measured Volovik (2013).
Such spontaneous emergence of time periodic oscillations can be investigated
in a condensate formed by bosonic magnetic quasi-particles (magnons) excited
in superfluid 3He-B Volovik (2013) (see also Homann et al. (2020); Ojeda
Collado et al. (2021)). Application of a radio-frequency pulse deflects the
magnetization of the superfluid 3He-B and the resulting induction signal first
decays but afterwards spontaneously reappears with its own phase completely
forgetting the phase of the radio-frequency pulse. Spontaneous precession of
the magnetization in magnon condensates in the context of time crystals has
been experimentally investigated Autti et al. (2018); Kreil et al. (2019);
Autti et al. (2021). Note that when the total number of bosons is well
defined, a Bose-Einstein condensate is identified with a single non-vanishing
eigenvalue of the reduced single-particle density matrix and the U(1) symmetry
broken approach cannot be used to define a time crystal.

Classical time crystals proposed by Shapere and Wilczek are time-independent
systems that can perform periodic motion even if their energy is the lowest
possible Shapere and Wilczek (2012, 2019); Das et al. (2018). However, the
research on quantum discrete time crystals also inspired investigation of
classical discrete time crystals. The latter are periodically driven classical
many-body systems which break ergodicity and perform motion with a period
longer than the period of the drive Kim et al. (2006); Heo et al. (2010); Yao
et al. (2018); Simula (2022).

We have sketched a brief history of the new era of time crystals. More details
can be found in review articles Sacha and Zakrzewski (2018); Else et al.
(2020); Khemani et al. (2019); Guo and Liang (2020) and in recent books Sacha
(2020); Guo (2021).

Current research directions. – The existence of discrete time crystals in
closed systems is surprising because a generic periodically driven quantum
many-body system is normally expected to heat up to an infinite temperature
state where no regular temporal structure can be observed. There is currently
intensive research concerning the lack of heating and hence the ergodicity
breaking by discrete time crystals in closed systems, which we will describe
shortly.

Discrete time crystals in spin-based or lattice systems can be roughly divided
into two groups: (i) prethermal time crystals and (ii) time crystals protected
by many-body localization. In group (i) there are systems which are driven
with high frequency $\omega$ so that to absorb a single quantum portion of
energy from the time-periodic drive, a large number of local rearrangements in
the system is needed Abanin et al. (2017); Else et al. (2017); Mizuta et al.
(2019); Machado et al. (2020); Khemani et al. (2019). Consequently, they do
not heat up for a long time and only thermalize at a time that increases
exponentially with $\omega/J_{local}$, where the frequency $J_{local}$ is
associated with the local interaction energy in the system. In the group (ii),
there are systems whose static counterparts reveal many-body localization.
Time-independent quantum many-body systems in the presence of disorder can be
many-body localized. That is, they do not thermalize (in the sense of the
eigenstate thermalization hypothesis Deutsch (1991); Srednicki (1994)) because
there are local integrals of motion and they exhibit emergent integrability
Huse et al. (2014); Nandkishore and Huse (2015); Mierzejewski et al. (2018).
When a time-periodic driving is turned on and a discrete time crystal forms,
the crystal will, in principle, live forever because the ergodicity breaking
is protected by the many-body localization Khemani et al. (2016); Else et al.
(2016). While such a mechanism is plausible, one should bear in mind that
there is no rigorous proof that it is true in the thermodynamic limit and for
an infinite evolution time. Even in the case of time-independent systems,
many-body localization is a highly non-trivial phenomenon and there is only a
class of one-dimensional time-independent systems where its existence is
mathematically proven Imbrie (2016). It is very difficult to verify ergodicity
breaking in the laboratory because in all experiments performed so far there
were other processes, not included in the time crystal Hamiltonians, which
made the lifetime of the crystals finite (typically a hundred of the driving
periods Zhang et al. (2017); Choi et al. (2017); Pal et al. (2018); Rovny et
al. (2018); Smits et al. (2018); Kyprianidis et al. (2021); Randall et al.
(2021); Mi et al. (2022); Frey and Rachel (2022)). However, prethermal
properties of discrete time crystals not protected by many-body localization
have been demonstrated experimentally Kyprianidis et al. (2021).

The first discrete time crystal was proposed in ultra-cold bosonic atoms
bouncing on an oscillating atom mirror, i.e., in a system which is very
different from the spin-based or lattice systems Sacha (2015a). First, the
single-particle version of this system, contrary to a single spin, possesses a
well-defined classical limit. Second, this is a bosonic many-body system. The
lifetime of an ordinary space crystal is infinite if first the thermodynamic
limit is taken, and only then the infinite evolution time is considered.
Otherwise, a sufficiently long quantum evolution of a space crystal results in
spreading of the center of mass position of the system and consequently
blurring of the crystalline structure in space Sacha (2020). In the case of
ultra-cold bosonic atoms bouncing on a mirror, if the limit of an infinite
number of the particles is taken, the description of a discrete time crystal
reduces to the mean-field Gross-Pitaevskii equation Sacha (2015a); Wang et al.
(2021a, b). Period-doubling solutions of the Gross-Pitaevskii equation are
stable and consequently the discrete time crystal does not decay Sacha
(2015a).

The periodic evolution of discrete time crystals in closed systems implies
that the average transfer of energy from the drive to the systems is zero. For
these systems, the environment is an enemy which can destroy the time crystals
Lazarides et al. (2015). In the case of dissipative discrete time crystals the
environment participates in the formation of the crystals Gong et al. (2018).
The balance between energy pumped by an external drive and energy released to
the environment is responsible for the formation of periodically evolving
steady states. Treating the environment as an ally is thus another strategy to
realize time crystals, which can be used in future practical applications. Two
experiments demonstrating dissipative discrete time crystals have been
performed recently, in ultra-cold atoms Keßler et al. (2021) and in an optical
system Taheri et al. (2022). In the latter case, one can realize big discrete
time crystals that can live forever, i.e., as long as the optical resonator is
pumped by laser beams. Moreover, the experiment was carried out at room
temperature and does not require advanced experimental equipment, and
consequently, transfer from the laboratory to the real world seems
straightforward.

One may ask the question why do we need time crystals? They may be interesting
from the point of view of basic research, but what are their practical
applications? Discrete time crystals are many-body systems that can maintain
coherent evolution for a long time. This is a very important property which
can be used in future applications. Whether this will allow one to build
quantum computers, we shall see. Will time crystals allow us to build precise
clocks? At least in the case of discrete time crystals, the precision of their
periodic evolution is limited by the precision of the external periodic drive.
However, the fact that they evolve with subharmonic frequencies can be used as
a frequency division, which in the optical frequency regime, like in the
experiment Taheri et al. (2022), is highly non-trivial. In order to think
about practical applications of time crystals, they must be easy to implement
and well controlled experimentally. Therefore, we await the development of
experimental techniques and new groundbreaking realizations of time crystals.

All discrete time crystals realized in spin-based systems and ultra-cold atoms
so far have been small in the time dimension Zhang et al. (2017); Choi et al.
(2017); Pal et al. (2018); Rovny et al. (2018); Smits et al. (2018);
Kyprianidis et al. (2021); Randall et al. (2021); Keßler et al. (2021); Mi et
al. (2022); Frey and Rachel (2022). Typically, they evolved with a period
twice as long as the driving period. In other words, they consisted of only
two elementary cells in time. With such time crystals it is not possible to do
condensed matter physics in the time domain. However, discrete time crystals
that are big in the time dimension can be realized in ultra-cold atoms
Giergiel et al. (2018a); Giergiel et al. (2020); Pizzi et al. (2019) and also
in optical systems Taheri et al. (2022). Especially the former system can host
a variety of solid-state phases which we are going to describe in the
following Sacha (2020); Guo (2021). The development of solid-state physics in
the time dimension is a natural next stage of time crystal research with
powerful potential for future practical applications. This research which is
in the stage of experimental implementation Giergiel et al. (2020) is
described below.

Condensed Matter Physics in Big Time Crystals. – Discrete time crystals
created by a Bose-Einstein condensate (BEC) of ultracold atoms bouncing on an
oscillating atom mirror can exhibit dramatic breaking of discrete time
translational symmetry, where the bouncing atoms can evolve with periods up to
about 100 times longer than the period of the driving mirror Giergiel et al.
(2018a); Giergiel et al. (2020). Such a system allows the creation of big time
crystals possessing a large number of temporal lattice sites and are suitable
for investigating condensed matter physics in the time dimension Sacha
(2015b); Sacha and Delande (2016); Delande et al. (2017); Matus and Sacha
(2019); Mierzejewski et al. (2017); Giergiel et al. (2019b, 2018b, a); Kosior
and Sacha (2018); Giergiel et al. (2021); Kuroś et al. (2021). These condensed
matter phenomena may or may not involve spontaneous breaking of the discrete
time translation symmetry. In the latter case, the formation of time lattices
resembles the situation of photonic time crystals that do not arise
spontaneously because they are externally imposed by periodic fabrication of
the refractive index in space Joannopoulos et al. (2008).

For ultra-cold atoms bouncing on an oscillating mirror, the time-periodic
driving function of the mirror can be expanded as a Fourier series
$f(t)=f(t+2\pi/\omega)=\sum_{k}f_{k}e^{ik\omega t}$, which by a proper choice
of the Fourier components $f_{k}$ allows one to engineer effective temporal
lattice potentials with nearly any shape. This provides a highly flexible
platform to realize a broad range of condensed matter phenomena in the time
dimension. We first consider the case of a single particle bouncing on a
mirror oscillating with period $T$. The classical effective Hamiltonian
describing the motion in the vicinity of a resonant trajectory evolving with
period $sT$ (where $s$ is integer), in the moving frame of the resonant
trajectory, has the form Sacha (2020)

$H_{eff}=\frac{P^{2}}{2m_{eff}}+V_{eff}\left(\Theta\right),$ (1)

where we use action-angle coordinates $(I,\theta)$ and all quantities are in
gravitational units Sacha (2015a). $P=I-I_{s}$, $I_{s}$ is the action for the
$s$ resonant trajectory, $\Theta=\theta-\omega t/s$, $m_{eff}$ is an effective
mass and $V_{eff}(\Theta)$ is the effective potential which depends on the
Fourier components $f_{k}$ of the time-periodic driving function $f(t)$ of the
mirror. For pure harmonic driving in time, $f(t)\propto\cos\omega t$, the
effective potential $V_{eff}\propto\cos(s\Theta)$. Then, in the quantum
description of the particle Sacha (2020); Giergiel et al. (2020), the eigen-
energies form energy bands and the corresponding eigenstates have the form of
Bloch waves in time when they are observed in the laboratory frame. If we are
interested in just the first energy band, the description reduces to the
tight-binding model Dutta, O. and Gajda, M. and Hauke, P. and Lewenstein, M.
and Lühmann, D.-S. and Malomed, B. A. and Sowiński, T. and Zakrzewski, J.
(2015) in which there are $s$ Wannier-like states $w_{j}(z,t)$ corresponding
to localized wave-packets evolving along the resonant trajectory with period
$sT$, in a similar way as for the problem of an electron in a periodic space
crystal.

For the many-body case of $N$ indistinguishable bosonic particles, if the
interaction energy per particle is much smaller than the energy gap between
the first and second energy bands of the corresponding single-particle
problem, we may restrict attention to the first energy band. The many-body
Floquet Hamiltonian then takes the form of an effective many-body Bose-Hubbard
Hamiltonian which describes $N$ bosons in a time lattice Sacha (2020)

$\hat{H}_{F}\approx-\frac{J}{2}\sum_{j=1}^{s}\left(\hat{a}_{j+1}^{\dagger}\hat{a}_{j}+H.c.\right)+\frac{1}{2}\sum_{i,j=1}^{s}U_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{j}\hat{a}_{i},$
(2)

where the bosonic annihilation operators $\hat{a}_{i}$ correspond to
periodically evolving Wannier-like wave-packets $w_{j}(z,t)$. $J$ is the
tunnelling amplitude of the particle between neighbouring wave-packets, and
$U_{ii}$ and $U_{ij}$ are the intra-site and inter-site interaction
coefficients, which are proportional to the interparticle contact interaction
$g_{0}(t)$, and hence to the s-wave scattering length, which can be modulated
in time by modulating the magnetic field of a Feshbach resonance.

In space crystals we are interested in the regular distribution of particles
in space at a fixed moment of time (i.e., at the detection moment). In time
crystals, the roles of space and time are exchanged, i.e., we fix a position
in space and ask if the probability of clicking of a detector located at this
point behaves periodically in time Sacha (2020). The Hamiltonian (2) describes
a many-body system loaded into a time lattice. The detector located close to
the resonant trajectory of this system will reveal a crystalline structure in
time, and the underlying temporal behavior is identical to the spatial
behavior of a solid state system described by the static version of the
Hamiltonian (2), i.e., when the operators $\hat{a}_{i}$ correspond to time-
independent Wannier modes of a space lattice Dutta, O. and Gajda, M. and
Hauke, P. and Lewenstein, M. and Lühmann, D.-S. and Malomed, B. A. and
Sowiński, T. and Zakrzewski, J. (2015). In the following we present a few
examples of non-trivial condensed matter phenomena which can be realized in
the time domain with the help of ultra-cold atoms bouncing on a mirror
Hannaford and Sacha (2020); Sacha (2020).

Condensed matter phenomena can also be investigated in photonic time crystals
which are temporal analogues of photonic space crystals but the refractive
index of the dielectric material changes periodically in time and not in space
Lustig et al. (2018); Sharabi et al. (2021); Galiffi et al. (2022). Another
direction of condensed matter research in periodically driven systems is
phase-space crystals Guo et al. (2013); Guo and Marthaler (2016); Guo et al.
(2016); Pengfei et al. (2018); Guo and Liang (2020); Guo et al. (2022); Guo
(2021). They reveal crystal structures in phase space where single-particle
and many-body solid state phenomena can be investigated.

Figure 1: Comparison of Anderson localization in a 1D space crystal with
periodic boundary conditions (left) with Anderson localization in the time
domain (right). A particle moving in the space crystal with disorder can
Anderson localize, i.e., if we go around the ring, we observe an exponential
localization of a particle around a certain point in space. When switching
from Anderson localization in space to Anderson localization in time, we have
to exchange the role of space and time. We fix the position in space and ask
if the probability of clicking of a detector is exponentially localized around
a certain moment of time — such behavior is repeated periodically due to the
periodic boundary conditions in time. This circulating ring system is
equivalent to an atom bouncing resonantly on an oscillating mirror. Reprinted
from Sacha and Zakrzewski (2018).

Anderson localization in the time dimension. Anderson localization normally
involves the exponential localization of eigenstates of a particle in space
due to the presence of a spatially disordered potential Anderson (1958).
Anderson localization can also be observed in the time domain due to the
presence of temporal disorder Sacha (2015b).

If a _single_ _particle_ is bouncing resonantly on a mirror which is driven in
time as $f(t)=\lambda\cos(\omega t)+\sum_{k}f_{k}e^{ik\omega t/s}$, where the
$f_{k}$ are random complex numbers with $|f_{k}|\ll\lambda$, there are
additional temporal disorder terms, $\sum_{j=1}^{s}\epsilon_{j}|a_{j}|^{2}$
(where the $\epsilon_{j}$ are real random numbers), in the tight-binding
Hamiltonian [i.e., the single-particle version of the Hamiltonian (2)], and we
observe an exponential localization in time, i.e., the probability of clicking
of a detector located close to the resonant trajectory is exponentially
localized around a certain moment of time Sacha (2015b); Giergiel et al.
(2018a), see Fig. 1. It is also possible to observe the localized-delocalized
Anderson transition in the time domain if a particle moving in 3D space is
properly driven Delande et al. (2017). Recently, signatures of temporal
Anderson localization have also been studied in photonic time crystals Sharabi
et al. (2021).

If a resonantly driven _many-body_ system experiences temporal disorder, there
are additional temporal disorder terms,
$\sum_{j=1}^{s}\epsilon_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}$ in the effective
Bose-Hubbard model (2), and the system can reveal many-body localization in
time Mierzejewski et al. (2017).

Topological time crystals. Topological insulators are insulators in the bulk
but possess topologically protected conducting edge states or surface states.

If a _single particle_ is resonantly bouncing on a mirror which is driven with
two harmonics, $f(t)=\lambda\cos(s\omega t)+\lambda^{\prime}\cos(s\omega
t/2)$, then the effective potential in Eq. (1) reveals a crystalline structure
with a two-point basis
$V_{eff}(\Theta)=V_{0}\cos(s\Theta)+V_{0}^{\prime}\cos(s\Theta/2)$. In the
quantum description, the first two energy bands of $H_{eff}$ can be described
by the Su-Schriefer-Heeger (SSH) model which is a tight-binding model with
staggered hopping amplitudes Su et al. (1979). By introducing a narrow barrier
in $V_{eff}(\Theta)$, by applying a properly chosen additional modulation
$\sum_{k}f_{k}e^{ik\omega t}$, we can create a topological system with an edge
and two eigenstates which are localized close to the edge and appear localized
in time when the particle is observed in the laboratory frame, see Fig. 2
Giergiel et al. (2019b).

Lusting et al. have considered topological photonic time crystals Lustig et
al. (2018). If an electromagnetic wave propagates in a time-varying medium so
that two photonic time crystals are realized one by one, and if the wavenumber
of the wave is located in a common gap of the crystals with different
topological properties, then a temporal analogue of a topologically protected
edge state can be observed.

In a _many-body_ system, if the interaction $g_{0}$ is modulated not only in
time but also in space, via a Feshbach resonance, there are additional
possibilities to engineer different effective long-range interactions. For
example, if $g_{0}\left(z,t\right)=\sum_{m=2}^{2}{\alpha_{m}(t)z^{m}}$, the
coefficients $\alpha_{m}(t)$ can be chosen so that the Bose-Hubbard model (2)
reduces to a form which can reveal the bosonic analogue of the topological
Haldane insulator in a spin-one chain Giergiel et al. (2019b).

Figure 2: A particle driven resonantly in time so that the effective
Hamiltonian (1) reproduces the Su-Schrieffer-Heeger model Su et al. (1979)
with an edge for the case $s=42$. In the laboratory frame, the bulk states are
delocalized along the entire resonant orbit while the topologically protected
edge states are localized close to the edge created in time. Reprinted from
Sacha (2020).

Mott insulator phases in the time dimension. If in a resonantly driven
$N$-atom system, Eq. (2), the effective repulsive on-site interactions
$U_{ii}$ dominate over the long-range repulsion $U_{ij\neq i}$ and are
sufficiently strong compared with the tunneling rate $J$, the system reveals a
Mott insulator phase in the time domain Sacha (2015b). Then, the gap between
the ground and excited states of the Bose-Hubbard Hamiltonian (2) is opened,
fluctuations of the number of atoms in each site are suppressed, and the
system is not compressible. When we look at a fixed position close to the
resonant trajectory, well-defined bunches of atoms (with number $N/s$) arrive
periodically at the observation point like on a conveyer belt and there is no
coherence between them Sacha (2015b). This is in contrast to the superfluid
phase where the repulsive interactions are weak and the ground state of the
Bose-Hubbard model (2) is a BEC. The superfluid-Mott insulator transition can
be realized in the time lattice described here by changing the tunneling rate
$J$ Sacha (2015b).

Exotic long-range interactions in the time dimension. If we periodically
modulate the contact interaction $g_{0}(t)$ between ultra-cold atoms by
modulating the magnetic field around a Feshbach resonance, we can make the
(normally negligible) long-range interactions $U_{ij\neq i}$ in the effective
Hamiltonian (2) significant and engineer how they change with distance $|i-j|$
between temporal lattice sites Giergiel et al. (2018b). If $g_{0}(t)$ changes
in time with the period $sT$, one can engineer exotic long-range interactions
including interactions not existing in nature, such as interactions that
change the repulsive and attractive character in an oscillatory way with
increasing $|i-j|$.

Figure 3: Trajectories of atomic wave-packets bouncing between two orthogonal
oscillating mirrors, located at $x=0$ and $y=0$, on a 2D $5\times 5$ time
lattice. Reprinted from Giergiel et al. (2018b).

Higher-dimensional time lattices. Time crystals can be extended to 2D (or 3D)
by allowing the BEC to bounce resonantly between two (or three) oscillating
atom mirrors Giergiel et al. (2018b).

We first consider the case of two _orthogonal_ oscillating mirrors, located at
$x=0$ and $y=0$, oscillating at frequency $\omega$ and inclined at 45o to the
gravitational force (Fig. 3). With this configuration, the localized wave-
packets are products of 1D wave-packets. If the $s$ resonance condition is
fulfilled for both mirrors, then a square $s\times s$ time lattice is created
(Fig. 3). Use of this 2D time lattice allows gradual breaking of the discrete
time translation symmetry when the motion along one of the two orthogonal
directions breaks the symmetry and with a range of system parameters the
motion along the other direction also reveals symmetry breaking Kuroś et al.
(2021).

For atoms bouncing between a pair of orthogonal atom mirrors, it is possible
to realize spontaneous formation of _time quasi-crystals_ when the ratio
$s_{x}/s_{y}$ (where $s_{x,y}$ are the resonance numbers corresponding to the
particle’s motion along the $x$ and $y$ directions) approximates the golden
number $(1+\sqrt{5})/2\approx 1.618$ Giergiel et al. (2019a). Then, when the
discrete time translation symmetry is spontaneously broken, the sequence of
atoms reflected from one mirror and the other mirror forms a Fibonacci quasi-
crystal in time.

Figure 4: (a) Honeycomb and (b) Lieb time lattices presented in the frame
moving along a resonant trajectory. Atoms loaded to lattices fulfill the
boundary conditions of the Möbius strip geometry, i.e., points
$\\{\Theta_{+}=\pi,\Theta_{-}\\}$ are identified with
$\\{\Theta_{+}=0,\pi-\Theta_{-}\\}$. Reprinted from Giergiel et al. (2021).

We now consider the case when the oscillating mirror at $y=0$ is moved to the
vertical position to form a 45o wedge. When a particle strikes the vertical
mirror, its momenta are exchanged $p_{x}\leftrightarrow p_{y}$, whereas when
it strikes the 45o mirror, the momentum $p_{y}$ remains the same while the
momentum $p_{x}$ is reversed $p_{x}\rightarrow-p_{x}$ Giergiel et al. (2021).
With this configuration, the system is inseparable and forms a time lattice
with a Möbius strip geometry. Using different relative amplitudes and phases
of the oscillations of the two mirrors we can realize different crystalline
structures of the effective potential, such as a honeycomb time lattice or a
Lieb square time lattice (Fig. 4) with a ‘flat band’, in which the dynamics of
the atoms is governed solely by interactions Giergiel et al. (2021).

Summary. – Seminal papers published by Wilczek and Shapere in 2012 initiated
the new era of time crystals. The idea of time crystals is inspiring
scientists and attracting public attention. We have briefly described the
history of time crystals and pointed out current research directions. Among
them, condensed matter physics in the time domain seems to be a promising
horizon for novel theoretical and experimental research. It is also a natural
next stage of development of time crystals toward practical applications.

###### Acknowledgements.

This work was supported by the Australian Research Council Discovery Grant No.
DP190100815 and the National Science Centre, Poland via Project No.
218/31/B/ST2/00319.

## References

  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.160401.
  * Shapere and Wilczek (2012) A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.160402.
  * Bruno (2013) P. Bruno, Phys. Rev. Lett. 111, 070402 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.111.070402.
  * Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.114.251603.
  * Syrwid et al. (2017) A. Syrwid, J. Zakrzewski, and K. Sacha, Phys. Rev. Lett. 119, 250602 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.250602.
  * Sacha (2015a) K. Sacha, Phys. Rev. A 91, 033617 (2015a), URL http://link.aps.org/doi/10.1103/PhysRevA.91.033617.
  * Khemani et al. (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016), URL http://link.aps.org/doi/10.1103/PhysRevLett.116.250401.
  * Else et al. (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016), URL http://link.aps.org/doi/10.1103/PhysRevLett.117.090402.
  * Yao et al. (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017), URL http://link.aps.org/doi/10.1103/PhysRevLett.118.030401.
  * 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., Nature 543, 217 (2017), ISSN 0028-0836, URL http://dx.doi.org/10.1038/nature21413.
  * 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., Nature 543, 221 (2017), ISSN 0028-0836, URL http://dx.doi.org/10.1038/nature21426.
  * Pal et al. (2018) S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys. Rev. Lett. 120, 180602 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180602.
  * Rovny et al. (2018) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180603.
  * Kyprianidis et al. (2021) A. Kyprianidis, F. Machado, W. Morong, P. Becker, K. S. Collins, D. V. Else, L. Feng, P. W. Hess, C. Nayak, G. Pagano, et al., Science 372, 1192 (2021), eprint 2102.01695.
  * Randall et al. (2021) J. Randall, C. E. Bradley, F. V. van der Gronden, A. Galicia, M. H. Abobeih, M. Markham, D. J. Twitchen, F. Machado, N. Y. Yao, and T. H. Taminiau, Science 374, 1474 (2021), eprint 2107.00736.
  * Mi et al. (2022) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, et al., Nature 601, 531 (2022), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-021-04257-w.
  * Frey and Rachel (2022) P. Frey and S. Rachel, Science Advances 8, eabm7652 (2022), URL https://www.science.org/doi/abs/10.1126/sciadv.abm7652.
  * Smits et al. (2018) J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Phys. Rev. Lett. 121, 185301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.185301.
  * Giergiel et al. (2018a) K. Giergiel, A. Kosior, P. Hannaford, and K. Sacha, Phys. Rev. A 98, 013613 (2018a), URL https://link.aps.org/doi/10.1103/PhysRevA.98.013613.
  * Surace et al. (2019) F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva, R. Fazio, and F. Iemini, Phys. Rev. B 99, 104303 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.99.104303.
  * Pizzi et al. (2019) A. Pizzi, J. Knolle, and A. Nunnenkamp, Phys. Rev. Lett. 123, 150601 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.150601.
  * Giergiel et al. (2020) K. Giergiel, T. Tran, A. Zaheer, A. Singh, A. Sidorov, K. Sacha, and P. Hannaford, New Journal of Physics 22, 085004 (2020), URL https://doi.org/10.1088/1367-2630/aba3e6.
  * Hannaford and Sacha (2020) P. Hannaford and K. Sacha, Physics World 33, 42 (2020).
  * Li et al. (2012) T. Li, Z.-X. Gong, Z.-Q. Yin, H. T. Quan, X. Yin, P. Zhang, L.-M. Duan, and X. Zhang, Phys. Rev. Lett. 109, 163001 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.163001.
  * Giergiel et al. (2019a) K. Giergiel, A. Kuroś, and K. Sacha, Phys. Rev. B 99, 220303 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevB.99.220303.
  * Kozin and Kyriienko (2019) V. K. Kozin and O. Kyriienko, Phys. Rev. Lett. 123, 210602 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.210602.
  * Nicolis and Prigogine (1977) G. Nicolis and I. Prigogine, _Self-organization in nonequilibrium systems : from dissipative structures to order through fluctuations / G. Nicolis, I. Prigogine_ (Wiley New York, 1977), ISBN 0471024015\.
  * Prigogine (1978) I. Prigogine, Science 201, 777 (1978), ISSN 0036-8075, URL https://science.sciencemag.org/content/201/4358/777.
  * Gong et al. (2018) Z. Gong, R. Hamazaki, and M. Ueda, Phys. Rev. Lett. 120, 040404 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.040404.
  * Iemini et al. (2018) F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio, Phys. Rev. Lett. 121, 035301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.035301.
  * Buča et al. (2019) B. Buča, J. Tindall, and D. Jaksch, Nature Communications 10, 1730 (2019), URL https://doi.org/10.1038/s41467-019-09757-y.
  * Dogra et al. (2019) N. Dogra, M. Landini, K. Kroeger, L. Hruby, T. Donner, and T. Esslinger, Science 366, 1496 (2019), URL https://www.science.org/doi/abs/10.1126/science.aaw4465.
  * Cosme et al. (2019) J. G. Cosme, J. Skulte, and L. Mathey, Phys. Rev. A 100, 053615 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.100.053615.
  * Keßler et al. (2021) H. Keßler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, Phys. Rev. Lett. 127, 043602 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.127.043602.
  * Taheri et al. (2022) H. Taheri, A. B. Matsko, L. Maleki, and K. Sacha, Nature Communications 13, 848 (2022), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-022-28462-x.
  * Volovik (2013) G. E. Volovik, JETP Lett. 98, 491 (2013).
  * Homann et al. (2020) G. Homann, J. G. Cosme, and L. Mathey, Phys. Rev. Research 2, 043214 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.043214.
  * Ojeda Collado et al. (2021) H. P. Ojeda Collado, G. Usaj, C. A. Balseiro, D. H. Zanette, and J. Lorenzana, Phys. Rev. Research 3, L042023 (2021), URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.L042023.
  * Autti et al. (2018) S. Autti, V. B. Eltsov, and G. E. Volovik, Phys. Rev. Lett. 120, 215301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.215301.
  * Kreil et al. (2019) A. J. E. Kreil, H. Y. Musiienko-Shmarova, S. Eggert, A. A. Serga, B. Hillebrands, D. A. Bozhko, A. Pomyalov, and V. S. L’vov, Phys. Rev. B 100, 020406 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.020406.
  * Autti et al. (2021) S. Autti, P. J. Heikkinen, J. T. Mäkinen, G. E. Volovik, V. V. Zavjalov, and V. B. Eltsov, Nature Materials 20, 171 (2021), ISSN 1476-4660, URL https://doi.org/10.1038/s41563-020-0780-y.
  * Shapere and Wilczek (2019) A. D. Shapere and F. Wilczek, Proceedings of the National Academy of Sciences 116, 18772 (2019), URL https://www.pnas.org/doi/abs/10.1073/pnas.1908758116.
  * Das et al. (2018) P. Das, S. Pan, S. Ghosh, and P. Pal, Phys. Rev. D 98, 024004 (2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.024004.
  * Kim et al. (2006) K. Kim, M.-S. Heo, K.-H. Lee, K. Jang, H.-R. Noh, D. Kim, and W. Jhe, Phys. Rev. Lett. 96, 150601 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.150601.
  * Heo et al. (2010) M.-S. Heo, Y. Kim, K. Kim, G. Moon, J. Lee, H.-R. Noh, M. I. Dykman, and W. Jhe, Phys. Rev. E 82, 031134 (2010), URL https://link.aps.org/doi/10.1103/PhysRevE.82.031134.
  * Yao et al. (2018) N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, ArXiv e-prints (2018), eprint 1801.02628.
  * Simula (2022) T. Simula, arXiv e-prints arXiv:2202.05407 (2022), eprint 2202.05407.
  * Sacha and Zakrzewski (2018) K. Sacha and J. Zakrzewski, Rep. Prog. Phys. 81, 016401 (2018), URL https://doi.org/10.1088/1361-6633/aa8b38.
  * Else et al. (2020) D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Annual Review of Condensed Matter Physics 11, 467 (2020), URL https://doi.org/10.1146/annurev-conmatphys-031119-050658.
  * Khemani et al. (2019) V. Khemani, R. Moessner, and S. L. Sondhi, arXiv e-prints arXiv:1910.10745 (2019).
  * Guo and Liang (2020) L. Guo and P. Liang, New Journal of Physics 22, 075003 (2020), URL https://doi.org/10.1088/1367-2630/ab9d54.
  * Sacha (2020) K. Sacha, _Time Crystals_ (Springer International Publishing, Switzerland, Cham, 2020), ISBN 978-3-030-52523-1, URL https://doi.org/10.1007/978-3-030-52523-1.
  * Guo (2021) L. Guo, _Phase Space Crystals_ , 2053-2563 (IOP Publishing, 2021), ISBN 978-0-7503-3563-8, URL https://dx.doi.org/10.1088/978-0-7503-3563-8.
  * Abanin et al. (2017) D. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, Communications in Mathematical Physics 354, 809 (2017), ISSN 1432-0916, URL https://doi.org/10.1007/s00220-017-2930-x.
  * Else et al. (2017) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. X 7, 011026 (2017), URL http://link.aps.org/doi/10.1103/PhysRevX.7.011026.
  * Mizuta et al. (2019) K. Mizuta, K. Takasan, and N. Kawakami, Phys. Rev. B 100, 020301 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.020301.
  * Machado et al. (2020) F. Machado, D. V. Else, G. D. Kahanamoku-Meyer, C. Nayak, and N. Y. Yao, Phys. Rev. X 10, 011043 (2020), URL https://link.aps.org/doi/10.1103/PhysRevX.10.011043.
  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991), URL https://link.aps.org/doi/10.1103/PhysRevA.43.2046.
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994), URL http://link.aps.org/doi/10.1103/PhysRevE.50.888.
  * Huse et al. (2014) D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.90.174202.
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat. Phys. 6, 15 (2015).
  * Mierzejewski et al. (2018) M. Mierzejewski, M. Kozarzewski, and P. Prelovšek, Phys. Rev. B 97, 064204 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.064204.
  * Imbrie (2016) J. Z. Imbrie, Journal of Statistical Physics 163, 998 (2016), eprint 1403.7837.
  * Wang et al. (2021a) J. Wang, P. Hannaford, and B. J. Dalton, New Journal of Physics 23, 063012 (2021a), URL https://doi.org/10.1088/1367-2630/abea45.
  * Wang et al. (2021b) J. Wang, K. Sacha, P. Hannaford, and B. J. Dalton, Phys. Rev. A 104, 053327 (2021b), URL https://link.aps.org/doi/10.1103/PhysRevA.104.053327.
  * Lazarides et al. (2015) A. Lazarides, A. Das, and R. Moessner, Phys. Rev. Lett. 115, 030402 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.030402.
  * Sacha (2015b) K. Sacha, Sci. Rep. 5, 10787 (2015b), URL https://www.nature.com/articles/srep10787.
  * Sacha and Delande (2016) K. Sacha and D. Delande, Phys. Rev. A 94, 023633 (2016), URL http://link.aps.org/doi/10.1103/PhysRevA.94.023633.
  * Delande et al. (2017) D. Delande, L. Morales-Molina, and K. Sacha, Phys. Rev. Lett. 119, 230404 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.230404.
  * Matus and Sacha (2019) P. Matus and K. Sacha, Phys. Rev. A 99, 033626 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.033626.
  * Mierzejewski et al. (2017) M. Mierzejewski, K. Giergiel, and K. Sacha, Phys. Rev. B 96, 140201 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.140201.
  * Giergiel et al. (2019b) K. Giergiel, A. Dauphin, M. Lewenstein, J. Zakrzewski, and K. Sacha, New Journal of Physics 21, 052003 (2019b), URL https://doi.org/10.1088/1367-2630/ab1e5f.
  * Giergiel et al. (2018b) K. Giergiel, A. Miroszewski, and K. Sacha, Phys. Rev. Lett. 120, 140401 (2018b), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.140401.
  * Kosior and Sacha (2018) A. Kosior and K. Sacha, Phys. Rev. A 97, 053621 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.053621.
  * Giergiel et al. (2021) K. Giergiel, A. Kuroś, A. Kosior, and K. Sacha, Phys. Rev. Lett. 127, 263003 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.127.263003.
  * Kuroś et al. (2021) A. Kuroś, R. Mukherjee, F. Mintert, and K. Sacha, Phys. Rev. Research 3, 043203 (2021), URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.043203.
  * Joannopoulos et al. (2008) J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, _Photonic crystals: Molding the flow of light (Second Edition)_ (Princeton University Press, 2008), 2nd ed., ISBN 0691124566\.
  * Dutta, O. and Gajda, M. and Hauke, P. and Lewenstein, M. and Lühmann, D.-S. and Malomed, B. A. and Sowiński, T. and Zakrzewski, J. (2015) Dutta, O. and Gajda, M. and Hauke, P. and Lewenstein, M. and Lühmann, D.-S. and Malomed, B. A. and Sowiński, T. and Zakrzewski, J., Reports on Progress in Physics 78, 066001 (2015), ISSN 0034-4885, URL http://stacks.iop.org/0034-4885/78/i=6/a=066001.
  * Lustig et al. (2018) E. Lustig, Y. Sharabi, and M. Segev, Optica 5, 1390 (2018), URL http://www.osapublishing.org/optica/abstract.cfm?URI=optica-5-11-1390.
  * Sharabi et al. (2021) Y. Sharabi, E. Lustig, and M. Segev, Phys. Rev. Lett. 126, 163902 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.163902.
  * Galiffi et al. (2022) E. Galiffi, R. Tirole, S. Yin, H. Li, S. Vezzoli, P. A. Huidobro, M. G. Silveirinha, R. Sapienza, A. Alù, and J. B. Pendry, Advanced Photonics 4, 1 (2022), URL https://doi.org/10.1117/1.AP.4.1.014002.
  * Guo et al. (2013) L. Guo, M. Marthaler, and G. Schön, Phys. Rev. Lett. 111, 205303 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.205303.
  * Guo and Marthaler (2016) L. Guo and M. Marthaler, New Journal of Physics 18, 023006 (2016), URL http://stacks.iop.org/1367-2630/18/i=2/a=023006.
  * Guo et al. (2016) L. Guo, M. Liu, and M. Marthaler, Phys. Rev. A 93, 053616 (2016), URL https://link.aps.org/doi/10.1103/PhysRevA.93.053616.
  * Pengfei et al. (2018) L. Pengfei, M. Michael, and L. Guo, New Journal of Physics 20, 023043 (2018), ISSN 1367-2630, URL http://stacks.iop.org/1367-2630/20/i=2/a=023043.
  * Guo et al. (2022) L. Guo, V. Peano, and F. Marquardt, Phys. Rev. B 105, 094301 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.094301.
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958), URL http://link.aps.org/doi/10.1103/PhysRev.109.1492.
  * Su et al. (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979), URL https://link.aps.org/doi/10.1103/PhysRevLett.42.1698.


# Fractal Spectrum of the Aubry-André Model

Ang-Kun Wu  Department of Physics and Astronomy, Center for Materials Theory,
Rutgers University, Piscataway, NJ 08854 USA

###### Abstract

The Aubry-André model is a one-dimensional lattice model for quasicrystals
with localized and delocalized phases. At the localization transition point,
the system displays fractal spectrum, which relates to the Hofstadter
butterfly. In this work, we uncover the exact self-similarity structures in
the energy spectrum. We separate the fractal structures into two parts: the
fractal filling positions of gaps and the scaling of gap sizes. We show that
the fractal fillings emerge for a certain type of incommensurate periodicity
regardless of potential strength. However, the power-law scaling of gap sizes
is characteristic for general incommensurability at the critical point of the
model.

Physical systems with incommensurability demonstrate novel physical
properties, e.g. the bilayer systems with incommensurate Moiré superlattices
Cano  _et al._ (2021); Wilson  _et al._ (2020); Fu  _et al._
(2020),topological materials driven by quasiperiodic potential Fu  _et al._
(2021); Chou  _et al._ (2020); Pixley  _et al._ (2018) and strongly correlated
quasicrystals An  _et al._ (2021) with non-Fermi-liquid behavior Thiem and
Chalker (2015); Deguchi  _et al._ (2012). In the exploration of correlation
effects in one-dimensional quasicrystals, we found that the incommensurability
between the discrete lattice and the onsite potential gives rise to self-
similar energy spectrum, which creates singularities and changes correlation
fundamentally Wu  _et al._ (2019). It turns out that the filling position of
conduction electrons is critical for the study of correlations. A casual
choice of the electron filling will gap the system by an emerging mini-gap
from increasing system sizes. Besides, in the quest of the impurity problem on
such one-dimensional quasicrystals, the numerical renormalization group asks
for the spectrum logarithmically close to the Fermi energy, which is
challenging with all the structures emerging from higher spectrum resolutions.
Thus, the understanding of the fractal spectrum in the thermodynamic limit
becomes important for the strongly correlated quasicrystals.

The study of the fractal spectrum could be traced back to the incommensurate
magnetic fields upon two-dimensional Bloch electrons, known as the Hofstadter
butterflyHofstadter (1976). In the problem, the Schrödinger equation turns
into the Harper’s equation Hsu and Falicov (1976); Harper (1955); Azbel’
(1964); Rauh (1974), a one-dimensional difference equation. By varying the
magnetic field, the periodicity of the equation is continuously changed and
the Hofstadter butterfly could be generated by the stacking the energy
spectrum for each field. Hofstadter, as well as othersAzbel’ (1964); Xu
(1986), summarized rules based on empirical examination of the numerical data
to construct the butterfly. From the “skeleton” of the butterfly, finer
structures could be created by two transformation rules, which give the
spectrum by the L(left) and R(right) type or the C(center) type of the
“skeleton”. These rules, in principle, could generate the butterfly and the
spectrum of a particular field indefinitely, but they requires the knowledge
of the spectrum of child field parameters for a parent field and lack the
knowledge of electron fillings. The Harper’s difference equation is strikingly
similar to the Aubry-André (AA) model at the localization critical point and
sometimes are also called Aubry-André-Harper model. The incommensurate
magnetic field problem is also studied on layered systems with Moiré
patternDean  _et al._ (2013); Bistritzer and MacDonald (2011); Lu  _et al._
(2021).

In this work, we aim for providing self-contained rules to reconstruct the
fractal spectrum of the one-dimensional quasicrystal, the AA model Azbel
(1979); Aubry and André (1980); Roati  _et al._ (2008), to the thermodynamic
limit. Instead of viewing the fractal spectrum empirically, we associated the
rational periodicity with the number of electrons in a unit cell and the
irrational periodicity will generate the mini-bands pattern scaling with the
system size. We argue that the self-similar spectrum exists only for a
particular type of incommensurability (irrational periodicity). The fractal
filling of electrons depends only on the irrational periodicity and exists
beyond the critical point. Morerover, the power-law decay of gap sizes is
characteristic at the critical point of the AA model and makes the spectrum
fractal for proper incommensurability.

This paper is organized as follow. First, we introduce the Aubry-André model
and the multi-band theory for rational periodicity. The organizing rules of
the spectrum are summarized for arbitrary periodicities. Second, we show the
fractal spectrum for a typical irrational periodicity and derive the self-
similarity in the filling positions of gaps from the organizing rules. We
argue that the self-similar spectrum is limited to a certain type of
irrational periodicities. Third, we identify the scaling behavior of the sizes
of the emerging mini-gaps across different phases of the AA model and show
that the power-law scaling is characteristic at the critical point.
Conclusions and discussions will be given in the end.

The model and organizing rules of bands. A general one-dimensional tight-bind
model with a constant hopping parameter and periodic onsite potential could be
written as

$\begin{split}H=-t\sum_{j}(c_{j}^{\dagger}c_{j+1}+h.c.)+\sum_{j}f(j)c_{j}^{\dagger}c_{j},\end{split}$
(1)

where $f(j)$ is some periodic function with rational or irrational frequency
$\alpha$ so that $f(j+1/\alpha)=f(j)$. A typical function is the cosine
function,

$\begin{split}f(j)=\lambda\cos(2\pi\alpha j+\phi),\end{split}$ (2)

known as the Aubry-André model, where $\alpha$ is an irrational number for
quasicrystals, i.e. $\alpha=(\sqrt{5}-1)/2$, and $\phi$ is a random phase that
does not change the energy spectrum fundamentally. Without particular mention,
the default $\phi$ is set to zero. $\lambda$ is the strength of the quasi-
periodic potential relative to the hopping ($t=1$ by default). The Fourier
transformation of the AA model reveals the self-duality of the model Aubry and
André (1980) between the momentum and real space, where $\lambda$ is mapped to
$2t$ in the momentum space. This self-duality identifies the localization
transition point at $\lambda=2t$. For $\lambda<2t$, all the wavefunctions are
localized in the momentum space while for $\lambda>2t$, all the wavefunctions
are localized in the real space.

Before we start to dive into the irrational part of the AA model, it is
relatively easy to understand this model with rational frequencies,
$\alpha=p/q$ with $p,q$ are coprime integers. The simplest example is
$\alpha=1/2$ and the onsite potentials oscillates between $-\lambda,\lambda$
along the chain. We could solve this model analytically by considering the
unit cell of length $2d$ (assuming the nearest neighbor distance is $d$) and
each unit cell has two atoms. If we denote the annihilation operators for the
two atoms as $a_{j},b_{j}$,

$H=-t\sum_{j}(a_{j}^{\dagger}b_{j}+b_{j}^{\dagger}a_{j+1}+h.c.)+\lambda[\sum_{j}a_{j}^{\dagger}a_{j}-b_{j}^{\dagger}b_{j}],$
(3)

whose continuum model could be written as bilinear $[a_{k},b_{k}]^{\prime}$
with $k$ as the momentum indices

$\begin{split}H(k)&=\begin{pmatrix}\lambda&-2t\cos(kd)\\\
-2t\cos(kd)&-\lambda\end{pmatrix}\end{split}.$ (4)

It is easy to see that there are two bands from the continuum spectrum
$E_{\pm}=\pm\sqrt{\lambda^{2}+4t^{2}\cos^{2}(kd)},$ with $k=\frac{2\pi
n}{Ld},kd\in[\frac{2\pi}{L},2\pi].$ The gap appears at a half filling with
size $\Delta=2\lambda$. For an arbitrary rational frequency $\alpha=p/q$, we
have corresponding $q-$band continuum in the thermodynamic limit:

$\begin{split}H(k)=\begin{pmatrix}V_{1}&-te^{ikd}&0&\cdots&-te^{-ikd}\\\
-te^{-ikd}&V_{2}&-te^{ikd}&0&\cdots\\\ \vdots&\vdots&\vdots&\vdots&\vdots\\\
-te^{ikd}&0&\cdots&-te^{-ikd}&V_{q}\end{pmatrix}\end{split},$ (5)

where $V_{j}=\lambda\cos(2\pi\frac{p}{q}j+\phi),i=1,2,\cdots,q$. For some
periodicity with a certain phase, e.g. $\alpha=1/4,\phi=0$, there are bands
merging due to the same potentials at $j=1+4n,3+4n$. However, this symmetry
can be broken by changing the phase $\phi$.

From the above unit cell viewpoint, rules could be summarized for rational
frequencies: (i) For $\alpha=p/q$, the periodicity of the potential is $q$ and
it should generate $q$ bands for a general phase $\phi$. Namely, electrons are
filled equally into $q$ bands with gaps appearing at fillings
$1/p,2/p,\cdots,(p-1)/p$. (ii) For arbitrary $\phi$, the models between
$\alpha,1-\alpha$ differ by a phase shift of $2\phi$. Namely, the spectrums of
$\alpha$ and $1-\alpha$ are symmetric. Combining the above two rules, the gap
filling positions could be summarized for an arbitrary $\alpha$ whose
fractional representation is unknown, or irrational $\alpha$: (iii) The
spectrum is separated by bands of fillings
$\\{\alpha\\},\\{2\alpha\\},\cdots;\\{1-\alpha\\},\\{2(1-\alpha)\\}\cdots$
with $\\{\cdots\\}$ denotes the fractional part of the number.

Figure 1: The DOS and nearest state gaps of the AA model at different phases
with system size $L=10,000,\alpha=(\sqrt{5}-1)/2$ and $\phi=0$. Upper panel:
the density of states (DOS). Lower panel:the nearest state gaps with peaks
locating around
$i=[N\\{\alpha\\}],[N\\{2\alpha\\}],\cdots;[N\\{1-\alpha\\}],[N\\{2(1-\alpha)\\}],\cdots$,
where $[\cdots]$ denotes the largest integer smaller the number.

The rule (iii) is the most important rule for our later discussion about the
gap size scaling. Empirically, the sizes of gaps appear in a descending manner
by the above order. For instance, the largest gaps start at fillings
$\\{\alpha\\},\\{1-\alpha\\}$, and smaller gaps at
$\\{2\alpha\\},\\{2(1-\alpha)\\}$ and so on. By definition, for an irrational
$\alpha$, the two filling series will continues indefinitely and give infinite
number of gaps in the thermodynamic limit, regardless of the quasi-periodic
strength $\lambda$. As shown in Fig. 1, the filling of the first several sharp
gaps coincide at the same position for different $\lambda$ with the same
frequency. However, the densities of states (DOS) are quite different between
the critical point ($\lambda=2$) and the other two phases. Note that the sizes
of the nearest state gaps are quite different for different phases–the gap
sizes decay in a slower fashion at the critical point.

The self-similarity of the AA model at the critical point is indicated in Fig.
2. If we zoom in the center part of the full spectrum and rescale the center
part to the whole band size (red curves), there is a clear reconstruction of
the whole spectrum from the middle part in terms of both the DOS and the
nearest state gaps. Mismatches of the gap sizes are due to the finite lattice
size, especially the smaller size of the center band.

Figure 2: The full fractal spectrum at the critical point with system size
$L=10,000,\alpha=(\sqrt{5}-1)/2$ and averaged over $100$ random phases. Upper
panel: the density of states (DOS) of the full spectrum and the scaled center
spectrum (filling=$[0.3821,0.6181]$ with edge energy $\pm 0.1891$ and scaled
by factor $13.74=2.5975/0.1891$). Lower panel:the nearest state gaps of the
full and scaled spectrum.

The fractal filling positions of gaps. Now that the organizing rules of the
electron fillings is summarized, the next step is to uncover the
transformation of the self-similar spectrum in terms of the filling positions
of gaps. Since the filling positions depend on the frequency $\alpha$, the
self-similar relation of the fractal spectrum also depends on $\alpha$.
Considering the golden ration frequency $\alpha=(\sqrt{5}-1)/2$, the fractal
fillings could be generated by the following transformation steps:

  1. 1.

Take the initial spectrum as three bands separated by two gaps at fillings
$\alpha,1-\alpha$, two symmetric golden ratio fillings.

  2. 2.

Transform each current band into three more bands with relative filling
positions of the two gaps at $\alpha,1-\alpha$.

To be quantitative, the gap positions could be classified into layers: 2 gaps
in the first layer ($\alpha,1-\alpha$); $6$ gaps in the second layer
($\\{2\alpha\\},\\{3\alpha\\},\\{4\alpha\\},\\{2(1-\alpha)\\},\\{3(1-\alpha)\\},\\{4(1-\alpha)\\}$);
$18$ gaps in the third layers and so on. Thus, the $k$th layer will have
$3^{k-1}\cdot 2$ new gaps and spectrum with $k$ layers will have $3^{k}-1$
gaps in total.

Note that $\\{n\alpha\\}=1-\\{n(1-\alpha)\\}$ due to the identity
$\\{n\alpha\\}+\\{-n\alpha\\}=1$. If we sort the second layer,
$\\{3(1-\alpha)\\}\approx 0.1459,\\{2\alpha\\}\approx
0.2361,\\{4\alpha\\}\approx 0.4721,\\{4(1-\alpha)\\}\approx
0.5279,\\{2(1-\alpha)\\}\approx 0.7639,\\{3\alpha\\}\approx 0.8541$, each pair
locates within the previous bands
$[0,(1-\alpha)],[(1-\alpha),\alpha],[\alpha,1]$. Thanks to the golden ratio
relation, $\alpha^{2}=1-\alpha$, the second layer could be rewritten as

$\begin{split}\\{3(1-\alpha)\\}&=2-3\alpha=(1-\alpha)^{2},\\\
\\{2\alpha\\}&=2\alpha-1=(1-\alpha)\alpha,\\\
\\{4\alpha\\}&=4\alpha-2=(1-\alpha)+(2\alpha-1)(1-\alpha),\\\
\\{4(1-\alpha)\\}&=3-4\alpha=(1-\alpha)+(2\alpha-1)\alpha,\\\
\\{2(1-\alpha)\\}&=2-2\alpha=\alpha+(1-\alpha)^{2},\\\
\\{3\alpha\\}&=3\alpha-1=\alpha+(1-\alpha)\alpha.\end{split}$ (6)

From the above expression, we could see that the second layer of gaps happen
to locate at the relative position of $1-\alpha,\alpha$ of the previous band,
which corresponds to one transformation of the previous band. The above six
gaps are all at the second order positions with highest order $\alpha^{2}$.
The next two gap positions $\\{5\alpha\\}=5\alpha-3=0+(1-\alpha)^{2}\alpha$,
$\\{5(1-\alpha)\\}=4-5\alpha=\\{3\alpha\\}+(1-\alpha)^{3}$ are in the third
layer. The gaps $\\{n\alpha\\},\\{n(1-\alpha)\\}$ belong to $k$th layer if the
integer $n$ satisfies

$\frac{3^{k-1}-1}{2}<n<=\frac{3^{k}-1}{2}.$ (7)

This self-similarity could also be confirmed by the two transformation
function of the Hofstadter butterfly, which only works at the critical point.
In the fractal butterflyHofstadter (1976), the spectrum of frequency $\alpha$
could be separated into three parts: the left (L), right (R) and the center
(C) parts, where each part could be a scaled spectrum of another two
frequencies $\alpha^{\prime}$(L,R parts),$\beta^{\prime}$(C part) based on the
two functions

$\begin{split}\alpha^{\prime}&=\frac{1}{\alpha}-[\frac{1}{\alpha}]=\frac{\sqrt{5}-1}{2}=\alpha\\\
\beta^{\prime}&=(\frac{1}{\alpha}-2)^{-1}=-\frac{\sqrt{5}+3}{2}\to\frac{1-\sqrt{5}}{2}=1-\alpha,\end{split}$
(8)

where $[x]$ denotes the greatest integer less than or equal to $x$. Thus, with
the golden ratio frequency, the new L, R and C parts are generated by the
same/conjugate incommensurate $\alpha$ and the self-similar transformations
are confirmed.

Now that the fractal fillings for $\alpha=(\sqrt{5}-1)/2$ is revealed, can we
extend the self-similar transformations to arbitrary irrational frequencies
$\alpha$? From the rules of gap filling positions, the frequency should be
able to transform $\\{n\alpha\\}$ into higher order representation
$\alpha^{k}$, namely the irrational frequency for a fractal spectrum should be
at least a solution of a polynomial equation. From the perspective of the
butterfly transformations, the irrational frequency should return to the same
value after several iterations of transformations. Another example for the
fractal spectrum is $\alpha=\frac{\sqrt{2}-1}{2}$. We could find that self-
similar pattern formed around the center of the spectrum (the gap position
$2\alpha=\frac{1}{2}-(\frac{1}{2}-\alpha)^{2}$, or
$\alpha^{\prime},\beta^{\prime}=\alpha$ after two transformation of the two
function.) However, the frequency $\alpha=\ln(2)$ could never give such self-
similar relations.

The scaling of gap sizes. From the above analysis, the fractal filling
positions of gaps only concerns the special irrational numbers. However, the
fractal spectrum only manifests itself at the critical point $\lambda=2$ of
the AA model. The scaling of gap sizes with higher multiples $n$ (or higher
layers $k$) for various strength $\lambda$ becomes critical. From the self-
similar relations at the critical point, we could derive that the gap sizes
should follows a power-law decay. From the initial spectrum, the three bands
(L, C, R parts) are just scaled spectrum of the whole spectrum by two
empirical factors

$\begin{split}f_{L/R}&=\frac{E_{max}-E_{side}}{2E_{max}}=0.0729,\\\
f_{C}&=\frac{E_{mid}}{E_{max}}=0.1392,\end{split}$ (9)

where $E_{max}=2.5975$,$E_{side}=1.8742$, $E_{mid}=0.18935$ denotes the half
band width of the whole spectrum, starting energy of the side band and the
half band width of the center band, respectively. Each new layer of gaps could
be generated by multiplying $f_{L/R}$ of the previous gaps for L, R parts and
$f_{C}$ for C part. Since each current gap will be transformed into three new
gaps, we could verify the scaling by the self-consistent relation:

$\begin{split}E_{max}&=\sum_{k=0}^{\infty}(E_{side}-E_{mid})(2f_{L/R}+f_{C})^{k}\\\
&=\frac{E_{side}-E_{mid}}{1-(2f_{L/R}+f_{C})},\end{split}$ (10)

where $(E_{side}-E_{mid})$ is the size of the first two gaps.

If we sort the nearest state gaps by descending order with index $n$, the
scaling of gap size are shown in Fig. 3. In both the localized and delocalized
phase, the gap sizes quickly converge to the generic energy resolution
$\delta=1/L$. The critical point gives a clear power-law scaling of gap sizes
to a scale lower than the energy resolution, with a power law

$\delta_{n}\sim n^{-1.98}.$ (11)

Note that the gap sizes decay in an exponential fashion as a function of the
layer order $k$, but the sorted order index $n$ relates to the layer order in
a logarithmic way, $k=1+[\log_{3}(n)]$, for $\alpha=(\sqrt{5}-1)/2$. In the
large $n$ (or $k$) limit, the dominating term is
$(\sqrt[3]{f_{L/R}^{2}f_{C}})^{k}$, which gives

$\delta_{n}\sim(f_{L/R}^{2}f_{C})^{k/3}\sim n^{-2.19}.$ (12)

The deviation of the power-law comes from the finite layer
$k=[\log_{3}(1,000)]=6$ and the finite lattice size. The power-law scaling of
gaps makes it possible for the self-similar transformations among small and
large scales through simple multiplication factors. However, this power-law
behavior does not limited by this particular irrational frequency. It appears
at different irrational frequencies, with and without self-similar structures,
as shown in the lower panel of Fig. 3. In other words, the power-law scaling
is characteristic at the critical point while the fractal spectrum is not
universal.

Figure 3: The scaling behavior of sorted nearest state gaps of the half
spectrum. Only the first $1,000$ gap sizes are shown out of total $L/2=5,000$
gaps. The gaps are averaged over $100$ random phases $\phi$. The horizontal
dashed line marks $\delta=1/L$, the energy resolution scale. The dashed lines
are linear fits of the sorted gaps as a function of order index. Upper panel:
The gap scaling in different phases. Lower panel: The gap scaling of various
frequency at the critical point $\lambda=2$.

In conclusion, we reveal the fractal spectrum of the AA model from the aspect
of the fractal filling. When a periodic potential is applied to a lattice
system, the band will be separated into multiply bands according to the
periodicity of the potential, which could be mapped to a multi-atom system
with translational symmetry. When the periodicity is an irrational number, the
translational symmetry is broken and the system is separated into infinite
mini-bands. The electron fillings of these mini-bands are just multiples of
the frequency $\alpha$ and its counterpart $1-\alpha$ with a fold into the
range $[0,1]$.

With specific irrational frequency, e.g. $\alpha=(\sqrt{5}-1)/2$, the filling
positions of gaps are fractal: more gaps could be generated by iteratively
inserting large scale structure into mini-bands. However, these fractal
fillings are limited to certain types of irrational numbers, solutions of
polynomial equations.

The fractal fillings depend only on the irrational frequencies of the
potential, but the fractal spectrum requires the mini-gap sizes decay in a
power-law fashion which depends on the strength of the quasi-periodic
potential. The critical point of the AA model happens to give rise this power-
law decay of gap sizes.

This work separates the study of the fractal spectrum into two independents
parts. The fractal spectrum of the AA model requires both a proper
incommensurate frequency and the power-law decay of the emerging gaps. The
first gives the fractal filling of the spectrum; the second makes the fractal
filling manifested in the DOS. This study of the fractal spectrum shall shad
lights on studies of incommensurate systems.

Acknowledgements. The author is grateful to Jedediah H. Pixley for encouraging
this project as an independent work. The author also thanks Eric Gawiser for
taking this project as a course project in his astrophysics class and thanks
Ghanashyam Khanal for coding neural networks to conduct density estimation on
the fractal DOS. This work was supported by NSF CAREER Grant No. DMR1941569.

## References

  * Cano  _et al._ (2021) J. Cano, S. Fang, J. H. Pixley, and J. H. Wilson, Phys. Rev. B 103, 155157 (2021).
  * Wilson  _et al._ (2020) J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley, Phys. Rev. Research 2, 023325 (2020).
  * Fu  _et al._ (2020) Y. Fu, E. J. König, J. H. Wilson, Y.-Z. Chou, and J. H. Pixley, npj Quantum Materials 5, 1 (2020).
  * Fu  _et al._ (2021) Y. Fu, J. H. Wilson, and J. H. Pixley, Phys. Rev. B 104, L041106 (2021).
  * Chou  _et al._ (2020) Y.-Z. Chou, Y. Fu, J. H. Wilson, E. J. König, and J. H. Pixley, Phys. Rev. B 101, 235121 (2020).
  * Pixley  _et al._ (2018) J. H. Pixley, J. H. Wilson, D. A. Huse, and S. Gopalakrishnan, Phys. Rev. Lett. 120, 207604 (2018).
  * An  _et al._ (2021) F. A. An, K. Padavić, E. J. Meier, S. Hegde, S. Ganeshan, J. H. Pixley, S. Vishveshwara, and B. Gadway, Phys. Rev. Lett. 126, 040603 (2021).
  * Thiem and Chalker (2015) S. Thiem and J. T. Chalker, Phys. Rev. B 92, 224409 (2015).
  * Deguchi  _et al._ (2012) K. Deguchi, S. Matsukawa, N. K. Sato, T. Hattori, K. Ishida, H. Takakura, and T. Ishimasa, Nature materials 11, 1013 (2012).
  * Wu  _et al._ (2019) A.-K. Wu, S. Gopalakrishnan, and J. H. Pixley, Phys. Rev. B 100, 165116 (2019).
  * Hofstadter (1976) D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * Hsu and Falicov (1976) W. Y. Hsu and L. M. Falicov, Phys. Rev. B 13, 1595 (1976).
  * Harper (1955) P. G. Harper, Proceedings of the Physical Society. Section A 68, 874 (1955).
  * Azbel’ (1964) M. Y. Azbel’, Eksp. Teor. Fiz. 46, 939 (1964).
  * Rauh (1974) A. Rauh, physica status solidi (b) 65, K131 (1974).
  * Xu (1986) H. Xu, Journal of Physics C: Solid State Physics 19, L795 (1986).
  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, _et al._ , Nature 497, 598 (2013).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Phys. Rev. B 84, 035440 (2011).
  * Lu  _et al._ (2021) X. Lu, B. Lian, G. Chaudhary, B. A. Piot, G. Romagnoli, K. Watanabe, T. Taniguchi, M. Poggio, A. H. MacDonald, B. A. Bernevig, _et al._ , Proceedings of the National Academy of Sciences 118 (2021).
  * Azbel (1979) M. Y. Azbel, Phys. Rev. Lett. 43, 1954 (1979).
  * Aubry and André (1980) S. Aubry and G. André, Ann. Israel Phys. Soc 3, 18 (1980).
  * Roati  _et al._ (2008) G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature 453, 895 (2008).