# GloVeInit at SemEval-2020 Task 1: Using GloVe Vector Initialization for
Unsupervised Lexical Semantic Change Detection

Vaibhav Jain
Delhi Technological University
Delhi, India
<EMAIL_ADDRESS>

###### Abstract

This paper presents a vector initialization approach for the SemEval2020 Task
1: Unsupervised Lexical Semantic Change Detection. Given two corpora belonging
to different time periods and a set of target words, this task requires us to
classify whether a word gained or lost a sense over time (subtask 1) and to
rank them on the basis of the changes in their word senses (subtask 2). The
proposed approach is based on using Vector Initialization method to align
GloVe embeddings. The idea is to consecutively train GloVe embeddings for both
corpora, while using the first model to initialize the second one. This paper
is based on the hypothesis that GloVe embeddings are more suited for the
Vector Initialization method than SGNS embeddings. It presents an intuitive
reasoning behind this hypothesis, and also talks about the impact of various
factors and hyperparameters on the performance of the proposed approach. Our
model ranks 13th and 10th among 33 teams in the two subtasks. The
implementation has been shared publicly.111github.com/vaibhav29498/GloVeInit

##  1 Introduction and Background

11footnotetext: This work is licensed under a Creative Commons Attribution 4.0
International License. License details:
http://creativecommons.org/licenses/by/4.0/.

Lexical Semantic Change (LSC) Detection is an active research topic in the
field of natural language processing, and has been applied for diachronic
(across time) and synchronic (across domains) tasks [Schlechtweg et al.,
2019]. Previously limited to manual ”close-reading” approaches, the
availability of large-scale corpora have allowed the use of computational
methods for this task [Tahmasebi et al., 2018]. This topic has found various
applications in various disciplines such as improving information retrieval
from historical documents [Morsy and Karypis, 2016], preventing cross-domain
ambiguity in requirements elicitation interviews [Jain et al., 2020], and
studying the impact of societal and cultural changes on word meanings and
usage [Tahmasebi and Risse, 2017].

LSC detection involves the use of two corpora $C_{1}$ and $C_{2}$ which, in
the diachronic case, belong to different time periods $t_{1}$ and $t_{2}$
respectively. The various approaches found in the literature usually involve
the construction of a word embedding space specific to each corpus. The
embeddings can be constructed through count-based methods such as Positive
Pointwise Mutual Information and Random Indexing, or predictive methods such
as Skip-Gram with Negative Sampling (SGNS). Most of these models are
stochastic in nature which means that separately trained embedding models live
in their own space. In order to project them onto a _unified space_ ,
alignment techniques such as vector initialization and orthogonal Procrustes
are used. The LSC of a word is then quantitatively determined by measuring the
contextual dissimilarity between the word’s representations [Tahmasebi et al.,
2018].

Recently, there have been efforts to evaluate these various methods by
comparing their results with manually-annotated data [Schlechtweg et al.,
2019, Ahmad et al., 2020]. The SemEval-2020 Task 1 is one such effort which is
based on LSC detection in corpora of English, German, Latin, and Swedish
languages [Schlechtweg et al., 2020]. Its aim is to provide an evaluation
framework for unsupervised LSC detection systems by comparing their results
against a ground truth, as annotated by native speakers or scholars. It
consists of two subtasks: Given two corpora $C_{1}$ and $C_{2}$ (for time
periods $t_{1}$ and $t_{2}$), for a set of target words,

  1. 1.

decide which words lost or gained senses between $t_{1}$ and $t_{2}$, and
which ones did not; as annotated by human judges.

  2. 2.

rank them according to their degree of LSC between $t_{1}$ and $t_{2}$ as
annotated by human judges. A higher rank means stronger change.

##  2 System Overview

###  2.1 Vector Initialization

The Vector Initialization (VI) alignment method, which was first used by ?),
involves the embedding space for $t_{1}$ to be trained independently on
$C_{1}$. It is then used to initialize the embedding space for $t_{2}$ which
is then subsequently trained from $C_{2}$. The underlying idea is that the
embedding for a word $w$ will get considerably updated if it is used within
different contexts in $C_{1}$ and $C_{2}$, otherwise it will receive only a
slight update. A study by ?), which applied the VI method on SGNS embeddings,
found it to perform significantly weaker on LSC detection tasks than other
methods such as orthogonal Procrustes (OP). They attributed this to the
sensitivity of the VI method to the frequency of a word in $C_{2}$. The high
frequency of a word in $C_{2}$ will result into the word’s embedding getting
frequent updates away from its initial state. This leads to a large divergence
between the word’s representations even if it does not undergo any LSC from
$t_{1}$ to $t_{2}$.

In a recent shared task on LSC detection in the German language, a modified VI
approach achieved the third-best result and outperformed certain OP approaches
[Ahmad et al., 2020]. Instead of only initializing on the word embeddings
obtained from $C_{1}$, the modified approach initializes on the complete SGNS
model which includes the hidden layer. However, it also theoretically suffers
from sensitivity to high frequency.

###  2.2 Comparison of GloVe and SGNS

GloVe (Global Vectors) and SGNS are unsupervised algorithms for obtaining word
embedding spaces [Pennington et al., 2014, Mikolov et al., 2013]. GloVe makes
use of co-occurrence matrix $X$; its $(i,j)$ entry, $X_{ij}$ is the number of
times the word $w_{j}$ appears in the context of the word $w_{i}$ (as defined
by the window-size $L$). It trains the word embeddings by minimizing the cost
function

$J=\sum^{V}_{i=1}\sum^{V}_{j=1;j\neq i}f(X_{ij})(u^{T}_{j}v_{i}-\log
X_{ij})^{2},$ (1)

where $V$ is the vocabulary size and $u,v\in\mathbb{R}^{D}$ are the word and
context word vectors respectively. The final word embeddings can be obtained
by summing or averaging the two.

On the other hand, SGNS is a predictive model which relies on a shallow two-
layer neural network which, given a word, predicts the set of its context
words. To avoid the expensive softmax function in the training objective,
negative sampling is used by drawing a few negative samples from a noise
distribution.

An important distinction between GloVe and SGNS from this paper’s point of
view is the number of updates that are made to a word embedding during
training. In a single epoch, the number of updates to the SGNS embedding of a
word is roughly equal to the number of words that appear in its context
throughout the corpus. Hence, the number of updates is proportional to the
frequency of the word which has an upper-bound of the total number of words in
the corpus. However, this relationship is not exactly linear because of
downsampling of frequent words and negative sampling. In the GloVe model, the
number of updates received in a single epoch is equal to the number of
_distinct_ context words. This number is limited by an upper-bound of
vocabulary size $V$, which is usually much less than the total number of words
in the corpus.

As discussed in Section 2.1, the VI method can falsely give a high LSC score
to words with high frequency. Our hypothesis is that the GloVe model is more
suited for the VI method due to its lesser sensitivity to high frequency. To
give an indication of the extent of this difference, the relationship between
a word’s frequency and the number of updates to its embedding is depicted in
Figure 1 for both GloVe and SGNS models trained on the $C_{2}$ corpora of all
the four languages. These results are based on the official GloVe
implementation222github.com/stanfordnlp/GloVe and the
gensim333radimrehurek.com/gensim implementation of SGNS. The models were
trained for a single epoch and only words with a frequency of greater than or
equal to 5 were considered. It is clear that the proportion of updates for
frequent and rare words is comparatively more balanced in the case of GloVe.
On the other hand, there is a high rate of growth in the number of updates
with respect to the frequency in the case of SGNS, which indicates that it is
biased towards estimating a high LSC for frequent words.

(a) CCOHA 1960-2010 (English)

(b) BZ and ND 1946-1990 (German)

(c) LatinISE 0-2000 (Latin)

(d) KubHist 1895-1903 (Swedish)

Figure 1: Plot of frequency of a word (x-axis) vs the number of updates to its
embedding (y-axis). The subcaptions include the name of the corpora $C_{2}$,
the corresponding time-period $t_{2}$, and its language.

##  3 Results and Experimental Setup

The dataset provided by the competition organizers consists of the corpora
$C_{1}$ and $C_{2}$, and the list of target words for four languages: English,
German, Latin, and Swedish [Schlechtweg et al., 2020]. All of the corpora were
already prepossessed: they were in tokenized and lemmatized form with
punctuation marks and one-word sentences removed.

Subtask 1 is concerned with identifying the loss or gain of one or more word
sense(s), where subtask 2 tests a model’s ability to detect fine-grained
changes in the two sense frequency distributions. For example, consider the
word _cell_ which historically referred to either a chamber or the smallest
unit of an organism. However, the relative usages of both these senses have
decreased with time, with mobile phone being the predominant meaning now. For
subtask 2, the ground-truth ranking of the target words is determined by
calculating the Jensen-Shannon divergence between their normalized sense
frequency distributions from $t_{1}$ and $t_{2}$ [Donoso and Sánchez, 2017]. A
submission is scored by its Spearman’s rank-order correlation coefficient
against the ground-truth ranking.

We used cosine distance as the metric to calculate the distance between a
word’s vectors in diachronic spaces, and every word with a distance of more
than 0.55 (for German) or 0.45 (for other languages) were classified to have
gained of lost word sense(s).

In all of our experiments, only words having a frequency of at least 5 were
considered. Our best-performing solution during the evaluation phase achieved
an accuracy of 60% in subtask 1 and a score of 0.352. It was obtained by
training word embeddings with a dimensionality of 50 on both $C_{1}$ and
$C_{2}$ for 60 epochs each with a window-size of 10. The subtask-2 score was
close to but less than the modified SGNS-based VI approach discussed in
Section 2.1, which was proposed by the team _IMS_ and received a score of
0.372. We ranked 13th and 10th out of 33 participants in the two subtasks
respectively.

Further experiments on subtask-2 were conducted after the competition for
analyzing the impact of hyperparameters like window-size $L$ and embedding
dimensionality $d$ and are reported in Table 1. Number of training epochs was
limited to 20 for models with window-size 10 because of their high
computational requirements. The results suggest that training the models for
higher number of epochs can produce better results if the embedding
dimensionality is high, but can backfire in the opposite case. Using a larger
window-size improves the average result in most cases.

L | Epochs | d | Scores
---|---|---|---
English | German | Latin | Swedish | Average
5 | 20 | 5 | 0.012 | 0.366 | 0.391 | 0.15 | 0.23
10 | -0.144 | 0.362 | 0.394 | 0.114 | 0.182
20 | 0.193 | 0.402 | 0.368 | 0.244 | 0.302
50 | 0.278 | 0.446 | 0.315 | 0.229 | 0.317
100 | 0.226 | 0.312 | 0.254 | 0.186 | 0.244
60 | 5 | -0.018 | 0.354 | 0.432 | 0.121 | 0.222
10 | -0.095 | 0.283 | 0.379 | 0.108 | 0.169
20 | 0.212 | 0.41 | 0.377 | 0.324 | 0.331
50 | 0.228 | 0.464 | 0.329 | 0.366 | 0.347
100 | 0.269 | 0.312 | 0.276 | 0.309 | 0.291
10 | 20 | 5 | 0.011 | 0.334 | 0.485 | 0.218 | 0.262
10 | -0.036 | 0.351 | 0.409 | 0.146 | 0.218
20 | 0.144 | 0.479 | 0.397 | 0.187 | 0.302
50 | 0.198 | 0.477 | 0.275 | 0.267 | 0.304
100 | 0.199 | 0.349 | 0.257 | 0.227 | 0.258
Table 1: Results of the post-evaluation experiments

###  3.1 Error Analysis

We defined the ranking error for each target word as the difference between
its predicted rank and its true rank divided by the total number of target
words. It lies in the range $(-1,1)$ and has an ideal value of zero. A
positive value indicates that the model overestimated the LSC for a word, and
a negative value indicates otherwise. A regression plot between the relative
frequency of the target words and their ranking error as per our best
performing model is depicted in Figure 2. A word’s relative frequency is
defined as its frequency divided by the total frequency of all words in the
corpus.

Figure 2: Regression plot between the frequency of the target words and their
ranking error

There is a statistically significant correlation between the ranking error and
the relative frequency, with the result of the Spearman’s rank-order
correlation test being $\rho=0.186$ and $p=0.020$. This indicates that the
GloVe model’s relative insensitivity to high frequency can lead to it
assigning a rank lower than the true rank to such words. We believe that the
frequency of the target words is not large enough for SGNS’ sensitivity to
high frequency to have a major impact on the results. This explains our
model’s inferior performance as compared to SGNS-based model in this task.

##  4 Conclusion and Future Work

In this paper, we reported our work in the SemEval2020 Task 1. We proposed a
GloVe-based VI approach which achieved the 11th rank out of 32 participating
teams. We gave a theoretical reasoning behind why GloVe is less sensitive
towards high frequency than SGNS and thus more suited for VI method, and
empirically showed the magnitude of the difference between the number of
updates to word embeddings in the two models. We believe that despite a lower-
than-expected performance in this competition, our work presents a good case
for the suitability of GloVe model when corpora of larger size are involved.
However, proving this quantitatively is a challenging task because of the
limitations associated with manual annotation-based evaluation.

Planned future work includes the study of how techniques such as dimension-
wise mean-centering and length normalization, which have proved beneficial in
OP-based approaches, can be applied for VI method.

## Acknowledgments

We would like to thank the task organizers for preparing the dataset and
organizing the competition, Jens Kaiser (member of the team _IMS_) from
University of Stuttgart for sharing the results from his work on the SGNS-
based VI alignment method, and the anonymous reviewers for their detailed and
helpful feedback.

## References

  * [Ahmad et al., 2020]  Adnan Ahmad, Kiflom Desta, Fabian Lang, and Dominik Schlechtweg.  2020\.  Shared task: Lexical semantic change detection in german.  CoRR, abs/2001.07786.
  * [Donoso and Sánchez, 2017]  Gonzalo Donoso and David Sánchez.  2017\.  Dialectometric analysis of language variation in twitter.  In Proceedings of the Fourth Workshop on NLP for Similar Languages, Varieties and Dialects (VarDial), pages 16–25, Valencia, Spain, April. Association for Computational Linguistics.
  * [Jain et al., 2020]  Vaibhav Jain, Ruchika Malhotra, Sanskar Jain, and Nishant Tanwar.  2020\.  Cross-domain ambiguity detection using linear transformation of word embedding spaces.  In Joint Proceedings of REFSQ-2020 Workshops, Doctoral Symposium, Live Studies Track, and Poster Track co-located with the 26th International Conference on Requirements Engineering: Foundation for Software Quality (REFSQ 2020), volume 2584 of CEUR Workshop Proceedings, Pisa, Italy. CEUR-WS.org.
  * [Kim et al., 2014]  Yoon Kim, Yi-I Chiu, Kentaro Hanaki, Darshan Hegde, and Slav Petrov.  2014\.  Temporal analysis of language through neural language models.  In Proceedings of the ACL 2014 Workshop on Language Technologies and Computational Social Science, pages 61–65, Baltimore, MD, USA, June. Association for Computational Linguistics.
  * [Mikolov et al., 2013]  Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean.  2013\.  Distributed representations of words and phrases and their compositionality.  In Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 2, NIPS’13, page 3111–3119, Red Hook, NY, USA. Curran Associates Inc.
  * [Morsy and Karypis, 2016]  Sara Morsy and George Karypis.  2016\.  Accounting for language changes over time in document similarity search.  ACM Trans. Inf. Syst., 35(1), September.
  * [Pennington et al., 2014]  Jeffrey Pennington, Richard Socher, and Christopher Manning.  2014\.  Glove: Global vectors for word representation.  In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 1532–1543, Doha, Qatar, October. Association for Computational Linguistics.
  * [Schlechtweg et al., 2019]  Dominik Schlechtweg, Anna Hätty, Marco Del Tredici, and Sabine Schulte im Walde.  2019\.  A wind of change: Detecting and evaluating lexical semantic change across times and domains.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pages 732–746, Florence, Italy, July. Association for Computational Linguistics.
  * [Schlechtweg et al., 2020]  Dominik Schlechtweg, Barbara McGillivray, Simon Hengchen, Haim Dubossarsky, and Nina Tahmasebi.  2020\.  SemEval-2020 Task 1: Unsupervised Lexical Semantic Change Detection.  In To appear in Proceedings of the 14th International Workshop on Semantic Evaluation, Barcelona, Spain. Association for Computational Linguistics.
  * [Tahmasebi and Risse, 2017]  Nina Tahmasebi and Thomas Risse.  2017\.  On the uses of word sense change for research in the digital humanities.  In Research and Advanced Technology for Digital Libraries, pages 246–257, Cham. Springer International Publishing.
  * [Tahmasebi et al., 2018]  Nina Tahmasebi, Lars Borin, and Adam Jatowt.  2018\.  Survey of computational approaches to lexical semantic change.  CoRR, abs/1811.06278.


#  Dynamical preparation of an atomic condensate in a Hofstadter band

Han Fu  James Franck Institute, University of Chicago, Chicago, IL 60637, USA
Department of Physics, University of Chicago, Chicago, IL 60637, USA  Andreas
Glatz  Materials Science Division, Argonne National Laboratory, Argonne, IL
60439, USA  Department of Physics, Northern Illinois University, DeKalb, IL
60115, USA  F. Setiawan  Pritzker School of Molecular Engineering, University
of Chicago, Chicago, IL 60637, USA  James Franck Institute, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA  Kai-Xuan Yao  James Franck Institute, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA  Enrico Fermi Institute, University of Chicago,
Chicago, IL 60637, USA  Zhendong Zhang  James Franck Institute, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA  Enrico Fermi Institute, University of Chicago,
Chicago, IL 60637, USA  Cheng Chin  James Franck Institute, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA  Enrico Fermi Institute, University of Chicago,
Chicago, IL 60637, USA  K. Levin  James Franck Institute, University of
Chicago, Chicago, IL 60637, USA  Department of Physics, University of Chicago,
Chicago, IL 60637, USA

###### Abstract

The creation of a Hamiltonian in the quantum regime which has non-trivial
topological features is a central goal of the cold-atom community, enabling
widespread exploration of novel phases of quantum matter. A general scheme to
synthesize such Hamiltonians is based on dynamical modulation of optical
lattices which thereby generate vector potentials. At the same time the
modulation can lead to heating and serious difficulties with equilibration.
Here we show that these challenges can be overcome by demonstrating how a
Hofstadter Bose-Einstein condensate (BEC) can be dynamically realized, using
experimental protocols. From Gross-Pitaevskii simulations our study reveals a
complex, multistage evolution; this includes a chaotic intermediate ``heating"
stage followed by a spontaneous reentrance to the BEC. The observed behavior
is reminiscent of evolution in cosmological models.

##  I Introduction

One of the great challenges in the field of ultracold atoms is to realize a
topological phase of a quantum, many-body system. While a number of novel
Hamiltonians have been realized experimentally Jotzu  _et al._ (2014);
Aidelsburger  _et al._ (2014); Miyake  _et al._ (2013), often based on
artificial gauge fields Lin  _et al._ (2009); Eckardt (2017); Goldman  _et
al._ (2014); Galitski  _et al._ (2019); Dalibard  _et al._ (2011); Struck  _et
al._ (2013); Goldman and Dalibard (2014); Dalibard (2015); Aidelsburger  _et
al._ (2018); Barbiero  _et al._ (2019), observing collective physics has
remained challenging. Creating a Hamiltonian with nontrivial topological
properties such as the iconic Harper-Hofstadter model Thouless  _et al._
(1982) and addressing it within the quantum regime will enable wide ranging
explorations of topological phases Kane and Mele (2005); Qi and Zhang (2011);
Hasan and Kane (2010). This has implications for atomic and condensed matter
physics and other sub-disciplines as well.

Almost all schemes for arriving at these Hamiltonians are dynamical in nature
Goldman and Dalibard (2014); Quelle  _et al._ (2017); Struck  _et al._ (2012);
Bukov  _et al._ (2015); Bukov and Polkovnikov (2014). In cold atomic gases
they involve the introduction of time-dependent optical lattices which
generate artificial gauge fields. Unfortunately, this dynamical engineering
has an important adverse consequence: heating Bilitewski and Cooper (2015);
Rubio-Abadal  _et al._ (2020), which presents impediments for reaching the
quantum regime. Successfully implementing this time-dependent or ``Floquet"
engineering in the quantum regime is, hence, a central goal of our larger
community. Most urgent is to identify the pathways involved in these
successful realizations. Importantly, there are no known fundamental barriers
for arriving in the quantum regime of the classic Hofstadter model. The MIT
groupKennedy  _et al._ (2015) has reported evidence for a Bose condensate in a
Hofstadter band. The Munich group Aidelsburger  _et al._ (2014), which has
simulated this Hamiltonian and observed topological features Aidelsburger  _et
al._ (2014), has, however, met difficulties in reaching the ground state.

It should be emphasized that Floquet generation of topological bands involves
multi-band participation with band crossing and inversion, abrupt band minima
transitions and complex patterns of flux penetration. In more conventional,
non-topological systems Eckardt (2017) in which new band minima are created by
Floquet engineering the changes in the bandstructure are more continuous. This
derives from the fact that a finite critical shaking amplitude is required to
obtain new band minima, which, in turn controls the onset of the equilibration
process. In the Hofstadter case, by contrast, an infinitesimally small shaking
amplitude will shift the ground state minima and abruptly initiate
equilibration. As a consequence we find that the resulting dynamics leads to a
chaotic, intermediate heating stage concomitant with the introduction of flux
and the related reorganization of the complex condensate phase pattern. The
novel dynamics makes BEC formation in Floquet-engineered topological systems
both richer and more complicated.

A central goal of this paper is to elucidate how such a superfluid can be
successfully realized given these complications. This involves characterizing
the requisite dynamical pathways. Instrumental to our paper is, then, the time
evolution. Here we emphasize our rather unexpected finding involving highly
chaotic behavior at intermediate times en route to re-condensation in a
topological band. We emphasize that understanding these phenomena should prove
an enormous benefit to the cold-atom, as well as to the solid-state Oka and
Kitamura (2019); Hunt  _et al._ (2013); Dean  _et al._ (2013); Ponomarenko
_et al._ (2013) and photonics communities Owens  _et al._ (2018) with a shared
interest in Floquet engineering.

In our theoretical investigation into the time sequence involved in Hofstadter
band condensation, we focus on three challenges which the system must address.
First, there is heating from the direct application of the Floquet drive. A
second challenge arises from the sudden change of the many-body ground state.
Meeting this challenge requires that the condensate wavefunction quickly
develop a specific and complex phase pattern. A third challenge comes from
accommodating interparticle interaction effects which are required for
equilibration, but not generally compatible with analytical predictions based
on Floquet engineering.

Figure 1: Characteristics of the Floquet-engineered Hofstadter model. (a)
Illustration of the ideal Hofstadter Hamiltonian which is the target model,
only approximately realized through a lattice-shaking protocol based on Eq.
(1) which is used here. $J_{x}^{\prime},\,J_{y}^{\prime}$ are the $x$ and $y$
axis tunneling parameters, and the phase $\phi^{y}_{i,j}=2\alpha\pi
i+\frac{\pi}{2}j$. The effective flux $\Phi$ inside each square cell is $1/4$
times the flux quantum. (b) Pre-shaking energy bandstructure at $k_{y}=0$
based on Eq. (1), where the energy $E$ is in recoil units $E_{R}$. Here
$\omega$ is the modulation frequency which couples the bands. Unless noted
otherwise, throughout the paper we use Aidelsburger  _et al._ (2014);
Bilitewski and Cooper (2015),
$V_{y}=6E_{R},\,V_{x}=10E_{R},\,V_{yl}=0.81E_{R}$, with shaking frequency
$\omega=0.72E_{R}/\hbar$, where the recoil energy
$E_{R}=\hbar^{2}(\pi/a)^{2}/2m$, and $a$ is the lattice constant of the
underlying square lattice, and $m$ is the atomic mass. (c) Floquet-engineered
Hofstadter bands at $\kappa=0.58\hbar\omega$. The energy separation between
states ${\bf{k}}=0$ and ${\bf{k}}=(\pm\pi/2a,0)$ is roughly $0.0025E_{R}$.
Here and throughout the paper, the ground state of this Floquet-engineered
Hofstadter bandstructure is indicated by `min'.

Our paper reports a rich set of dynamical processes en route to forming a BEC
in a Floquet-Hofstadter bandstructure. These observations are derived from
Gross-Pitaevskii (GP) simulations. Despite the aforementioned challenges, we
are able to provide a large body of evidence supporting the emergence of
condensation in a Hofstadter lattice. Moreover, our analysis shows that a
substantial fraction of the atoms are in the ground state. Interestingly we
observe a multistage dynamics, which has features in common with models of
cosmological evolution Allahverdi  _et al._ (2010). This starts with a
coherent series of oscillations of the population and is followed by a chaotic
``heating" stage, which is accompanied by an abrupt injection of magnetic flux
App ; Price  _et al._ (2016). We assume no dissipation in our simulations but,
surprisingly observe that the system is ultimately able to spontaneously relax
into the new ground state where condensation occurs. The intermediate heating
stage, in particular, is found to be essential, enabling reentrance to this
new condensate. When the system reaches the steady state, we are able to
extract occupations of different Hofstadter bands, thus characterizing a small
number of excitations which coexist with the BEC in the ground state.

##  II Floquet-Hofstadter theory

We follow the approach used by the Munich group Aidelsburger  _et al._ (2014);
Bilitewski and Cooper (2015) for Floquet engineering of the Hofstadter
Hamiltonian. This involves loading bosons into a two-dimensional optical
lattice with a potential which includes both a periodically oscillating
contribution $V_{os}$, and a static superlattice $V_{st}$. We consider a
square lattice having lattice constant $a$ superposed on an additional lattice
with constant $2a$ in the $y$ direction. The lattice potentials are given by

$\displaystyle V_{st}=$ $\displaystyle V_{x}\sin^{2}\left(\pi
x/a\right)+V_{y}\sin^{2}\left(\pi y/a\right)+V_{yl}\sin^{2}\left(\pi
y/2a\right),$ (1) $\displaystyle V_{os}=$
$\displaystyle\kappa\left[\sin\left(\pi/4+\pi
y/2a\right)\cos\left(\phi_{0}+\omega t-2\alpha\pi x/a\right)\right.$
$\displaystyle\left.+\cos\left(\pi/4+\pi y/2a\right)\sin\left(\phi_{0}-\omega
t-2\alpha\pi x/a\right)\right],$

where $V_{x},\,V_{y},\,V_{yl}$ represent the strengths of the respective
components in the static lattice, $\kappa$ is the amplitude of the oscillating
lattice which is identically zero before we turn on the shaking and
$\phi_{0}=\pi/4$. The presence of a time-dependent potential $V_{os}$ enables
the atoms to tunnel in the $y$ direction and acquire a position-dependent,
Aharonov-Bohm-like phase. Here $\alpha$ is the ratio between the flux per unit
cell of the square lattice $\Phi$ to the flux quantum $\Phi_{0}$, with
$\Phi/\Phi_{0}\equiv\alpha=1/4$. Our GP simulations are based on a Hamiltonian
which includes two body interactions and the single particle contributions
from both kinetic energy terms and lattice potentials which directly implement
Eq. (1). We consider a two-dimensional system.

Figure 2:  Abrupt changes in energy bands and wavefunctions due to shaking,
obtained from numerical calculations using Floquet theory. (a) Energy bands at
$k_{y}=0$ without shaking; $\bf{k}$ is the magnetic zone quasi-momentum. This
figure is identical to Fig. 1(b), now replotted in the framework of Floquet
theory. (b) Energy bands at $\kappa=0.1\hbar\omega$. Panels (c) and (d)
respectively show ground-state wavefunctions in the momentum ($\bf{p}$) space
at $\kappa=0$ and $0.1\hbar\omega$. These are labeled as `min' in the Floquet-
Hofstadter bands in (a) and (b). The abrupt change of characteristic momenta
as $\kappa$ varies reflects a first-order-like transition. Color codes
indicate contribution from the initial bands to the Floquet-Hofstadter bands;
intermediate colors in (b) represent band hybridization. Note that Fig. 1(c)
and Fig. 2(b) represent slightly different parameter sets, with the latter
chosen for pedagogical purposes to illustrate more clearly the dramatic change
in the band structure that ensues even at very small $\kappa$.

When the modulation energy $\hbar\omega$ is much larger than the effective
tunneling parameters, the system approaches the ideal Hofstadter Hamiltonian
Bukov and Polkovnikov (2014); Aidelsburger (2014); Bilitewski (2016). This
Hamiltonian has only nearest-neighbor tunneling in the $x$\- and $y$-direction
denoted by $J^{\prime}_{x}$ and $J^{\prime}_{y}e^{\imath\phi^{y}_{i,j}}$
respectively, where $\phi^{y}_{i,j}=2\alpha\pi i+\frac{\pi}{2}j$ and
$J_{x}^{\prime},\,J_{y}^{\prime}$ are real. The coordinates here are
$(x,y)=(i,j)a$. This ideal case, schematically illustrated in Fig. 1(a),
should be contrasted with the Floquet-Hofstadter realization based on Eq. (1).
Our theory implements the full dynamical Hamiltonian, which naturally includes
higher-order terms in $J^{\prime}_{x}/\hbar\omega$ and
$J^{\prime}_{y}/\hbar\omega$.

We characterize this latter Floquet-Hofstadter Hamiltonian through the
resulting band structure. In the absence of shaking ($\kappa=0$), the band
structure obtained from Eq. (1), is shown in Fig. 1(b). By contrast when
$\kappa$ assumes the experimental value Aidelsburger  _et al._ (2014);
Bilitewski and Cooper (2015) ($\kappa=0.58\hbar\omega)$ a very different band
dispersion emerges which is presented in Fig. 1(c).

It should be stressed that once a shaking amplitude $\kappa\neq 0$ is applied,
regardless of how small $\kappa$ is, there is a dramatic change of the ground
state. We contrast the bandstructures for the two situations: in the absence
of shaking [Fig. 2(a)] and at a small shaking amplitude
$\kappa=0.1\hbar\omega$ [Fig. 2(b)]. Here the energy minima in the ground band
shift their position abruptly from the band center ${\bf{k}}=0$ to the band
edge ${\bf{k}}=(\pm\frac{\pi}{2a},0)$. Note that, these two quasi-momenta
${\bf{k}}=(\pm\frac{\pi}{2a},0)$ are connected by a reciprocal vector and thus
correspond to a unique state; this is henceforth called the ``Floquet-
Hofstadter ground state" 111The Floquet-Hofstadter ground state is defined as
the lowest energy state in the ground band. This is the band in which the
${\bf{k}}=0$ state is adiabatically connected to the ground state in the
static system. .

The ground-state wavefunctions also exhibit a discontinuous change as can be
seen by comparing their behavior without shaking and with shaking at
$\kappa=0.1\hbar\omega$ through their distribution in momentum space, see
Figs. 2(c) and 2(d), respectively. In these figures the wavefunctions are
expanded in terms of ${\bf{p}}={\bf{k}}+{\bf{G}}$, where $\bf{G}$ are the
reciprocal wavevectors of the oscillating lattice potential. Before shaking,
the atoms are confined to ${\bf{p}}=(0,0)$. (The two extra spots appearing in
Fig. 2(c) are associated with higher reciprocal vectors.) In the presence of
lattice shaking a new set of four characteristic momenta emerge, represented
by ${{\bf{p}}=(\pm\frac{\pi}{2a},0),\,(0,\pm\frac{\pi}{2a}})$. Importantly, a
macroscopic population of these four momentum states can serve as a signature
that particles are occupying the ground state of the Floquet-Hofstadter band.

We simulate the dynamics of the atoms in the Floquet lattice through a GP
numerical procedure which uses a Graphics Processing Unit-based quasispectral,
split-step method to solve the GP equation based on fast Fourier transforms
Pathria and Morris (1990). Here we include a small, nonzero interparticle
interaction potential $U_{0}=7.5\times 10^{-4}E_{R}$. In our numerical
simulations, we start with a condensate in the static lattice $V_{st}$ and
linearly ramp up the shaking amplitude $\kappa$ in $V_{os}$ of Eq. (1). After
the ramp, $\kappa$ is held constant; from the GP simulations we are able to
examine the full evolution of the time-dependent wavefunction in both real and
momentum space.

Figure 3: Quench dynamics of Bose condensates undergoing transition from
conventional bands (with $\kappa=0$) to Hofstadter bands (at
$\kappa=0.58\hbar\omega$) with $\kappa$ linearly ramped from zero to
$0.58\hbar\omega$ within $30\,\mathcal{T}$ where $\mathcal{T}$ is the Floquet
period. (a) Time dependence of particle populations in four characteristic
momentum groups labeled using the color code in the inset: dashed line (red),
dash-dotted line (purple), thick solid line (green), and thin solid line
(blue). The transfer of boson populations between different groups indicates a
three-stage evolution. After initial oscillations in the first stage, a
`heating' state emerges which then spontaneously transitions to the final
condensation stage. (b) Particle population at $p_{y}=\pi/2a$ for
$t=720\,\mathcal{T}$ in the `heating' stage. The inset is the corresponding
image in the full momentum space within the same $2\pi/a\times 2\pi/a$
Brillouin zone as in the inset of (a), where the blue dashed box indicates the
relevant vertically integrated region. (c) Counterpart of (b) at
$t=1230\,\mathcal{T}$ in the condensation stage. The transition from broad
distribution in (b) to sharp peaks along the $p_{x}$ direction in (c) provides
some evidence for condensate formation. This analysis shows that a sizeable
(about $40\%$) fraction of the atoms is condensed.

##  III Characterizing the evolutionary pathways

Our simulations reveal a rich dynamics when the atoms transfer to the Floquet-
Hofstadter band. We observe three distinct evolutionary stages, as we follow
the momentum-space populations, see Fig. 3(a). Below we outline the key
features of each stage. In the first stage (from $0$ to $300\,\mathcal{T}$),
we see a period of coherent oscillations which involves transient occupations
of higher bands. A complicated dynamics then ensues within the second stage
(from $300$ to $900\,\mathcal{T}$). Here the population becomes widely
distributed over different states and different bands, see Fig. 3(b). We refer
to this second stage as the ``intermediate heating" stage, where interesting,
highly chaotic behavior occurs. This time period reflects the non-adiabatic
evolution and is reminiscent of the ``preheating" and ``turbulent" stages
associated with inflationary models of cosmology Allahverdi  _et al._ (2010);
Micha and Tkachev (2003, 2004). It is during this second stage, as the
wavefunction begins to develop a new and complex phase pattern, that we
observe a sudden onset of flux penetration App .

By contrast, in the third stage (beginning around $900\,\mathcal{T}$), the
population starts to settle into the Floquet-Hofstadter ground state. This
appearance of population accumulation into the ground state is suggestive of
Bose condensation. One sees that a rather sharp momentum distribution emerges
during this time, see, for example, Fig. 3(c). This matches that of the ground
state shown in Fig. 2(d). We emphasize that the evolution occurs spontaneously
in our simulations App , which are calculated without dissipation (consistent
with experimental conditions). Notably, this transition into the final ground
state is only possible in the presence of two-body interactions which drive
collisions and subsequent relaxation into a new set of momenta.

Figure 4: Comparison of real-space phase correlation function and vortex
structure between predicted Floquet-Hofstadter ground state (left) and final
dynamically evolved state from simulations (right) at $t=6000\,\mathcal{T}$.
(a) Absolute value of phase correlation functions given by
$g_{1}(\Delta{\bf{r}})=\langle e^{\imath\phi(i,j)}e^{-\imath\phi(i+\Delta
i,j+\Delta j)}\rangle$, where $\langle\dots\rangle$ denote averaging over
different ensembles and different $(i,j)$ positions with fixed relative
displacement $\Delta{\bf{r}}=(\Delta i,\Delta j)a$. Here $\phi(i,j)$ is the
local phase of the wavefunction at ${\bf{r}}=(x,y)=(i,j)a$. This shows a
finite spatial correlation length. (b) Distribution of vortices (blue dots)
and antivortices (red crosses). In the simulations, the checkerboard
arrangements are present in both the distribution of the phase correlations
and that of vortices. These are the predicted signatures of the Hofstadter
BEC.  Figure 5: (a) Stroboscopic plot of particle population $N_{\bf{p}}(t)$
at four characteristic momentum spots of the Floquet-Hofstadter ground state:
${\bf{p}}=(-\pi/2a,0)$ (dark blue: thick lower solid line), $(0,\pi/2a)$ (red:
thick upper solid line), $(\pi/2a,0)$ (light blue: thin lower solid line), and
$(0,-\pi/2a)$ (green: thin upper solid line). These spots correspond to the
red dots in the inset of Fig. 3(a). We define
$N_{\bf{p}}(t)=|\psi_{\bf{p}}(t)|^{2}$ where $\psi_{\bf{p}}$ is the
wavefunction expansion at momentum $\bf{p}$. Dashed lines show the average
values of the corresponding population curve with the same color. On average,
the two spots along the $y$ (or $x$) axis are approximately equally occupied.
The population of the two vertical spots is $1.70$ times that of the two
horizontal spots, close to the predicted value 1.67. (b) Power spectrum of the
wavefunction summed over the four characteristic spots
$N_{\textrm{tot}}(\varepsilon)=\sum_{\bf{p}}|\psi_{\bf{p}}(\varepsilon)|^{2}$,
where $\psi_{\bf{p}}(\varepsilon)$ is the Fourier transform of
$\psi_{\bf{p}}(t)$ from $t=5900$ to $6900\,\mathcal{T}$. At this time the
system appears to reach a dynamical equilibrium. The highest peak is
normalized to unity. The lower peaks (labeled 1, 2, 3) show weak occupation of
excited states whose quasimomenta are indicated in the inset.

##  IV Evidence for Condensation

Motivated by these suggestions of condensation in momentum space, and
recognizing that these sorts of analyses are complicated, we turn to more
direct evidence for a condensate through studies of spatial phase coherence.
We wait for time $6000\,\mathcal{T}$ when the atomic population is fully
settled into a steady state in momentum space. We then evaluate the phase
correlation in real space, and observe a long-range phase coherence which
extends over 10 sites, see Fig. 4(a). This provides additional and more direct
evidence for a BEC in the Floquet-Hofstadter ground state (referred to as a
``Hofstadter BEC"). At the same time we observe a well organized distribution
of vortices and antivortices. This distribution displays a checkerboard
pattern which matches theoretical predictions determined from the ground
state. Moreover, this vortex checkerboard pattern 222These vortices are
intrinsic and should be contrasted with vortex patterns associated with a
rotating trap as in Phys. Rev. Lett. 92, 020403. is associated with a similar
correlation length, see Fig. 4(b).

Following this strong evidence for a Hofstadter BEC, we next investigate in
which Floquet bands the condensate resides. To this end, we study the time-
dependent wavefunctions at the four characteristic momenta which appear as the
red spots in the inset of Fig. 3(a). The populations associated with these
four spots exhibit oscillatory behavior when viewed in a stroboscopic fashion
333The time slice within each cycle is chosen to be at zero phase such that
$t$ in Eq. (1) is $0+n\mathcal{T}$ where $n$ are integer multiples., see Fig.
5(a). Moreover, the time-averaged populations agree well with predictions App
derived from the ground state in the lowest Floquet-Hofstadter band.

These time-dependent oscillations reveal coherent superpositions involving
occupations of excited states. We are able to extract the energy spectrum of
these excitations from the Fourier transform of the wavefunctions in the time
domain, see Fig. 5(b). All Fourier spectra show the same set of peaks
consisting dominantly of the ground state with a few excited states. Comparing
with the Floquet-Hofstadter band structure, we can identify three excited
states from the energies of the weaker peaks App , see the inset of Fig. 5(b).
We speculate that, these excitations likely relate to the finite coherence
length and defects seen in Figs. 4(a) and 4(b); presumably they arise from the
nonadiabatic dynamics in the evolution to the Floquet-Hofstadter ground state.

##  V Dynamics of Hofstadter BEC formation

While the preceding sections have summarized our central results, it is useful
to understand in more detail the dynamics we observe in our GP simulations.
Particularly noteworthy in this regard is the very novel intermediate heating
stage.

###  V.1 Numerical details of the Gross Pitaevskii simulations

The numerical code we use employs Pathria and Morris (1990) GPU-based parallel
computing and is designed such that it conserves the particle number. The
general GP equation is:

$\displaystyle\imath\hbar\partial_{t}\psi({\bf{r}},t)=$ $\displaystyle
e^{\imath\gamma}\left[-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu+V+g_{int}|\psi({\bf{r}},t)|^{2}\right]\psi({\bf{r}},t),$
(2)

where the damping constant $\gamma$ is set to zero so that our simulations are
dissipationless. Here $\mu$ (set to unity) is the chemical potential,
$V=V_{st}+V_{os}$ is the total potential term (see Eq. (1) for $V_{st}$ and
$V_{os}$), and $g_{int}$ is the interaction strength which determines the
interaction energy $U_{0}=g_{int}n_{0}$. Here $n_{0}$ is the mean particle
density. Since $U_{0}$ is directly tuned in our simulations, we focus on this
parameter instead of $g_{int}$.

###  V.2  Magnetic-Brillouin-zone entropy

To more quantitatively characterize these dynamics, in addition to the
population curves shown in Fig. 3, we introduce an effective time-dependent
``entropy" $S_{\scriptscriptstyle{\mathrm{MBZ}}}$ calculated using states in
the first magnetic Brillouin zone. This serves to quantify the disorder in the
momentum distribution, and is defined by the occupation probability associated
with different ${\bf{p}}$ (momentum) spots. We caution that this ``entropy"
relates to how widely the particle distribution spreads in momentum space.
This does not represent a thermodynamical definition of entropy. It serves to
describe the sharpness of the momentum peaks over time. As such this
``entropy" can, at intermediate times, decrease with time. Eventually,
however, as in Fig. 6 the system enters the long time heating period where the
``entropy" monotonically increases.

Figure 6: Evolution of the Floquet system.  (a) Time-dependent momentum-space
distribution of the driven system obtained from simulations. (b) Entropy as a
function of time. Several different time frames are labeled by using the same
letters as in panel (a). Note that the longer-time heating, which is
associated with an upturn in the entropy, is evident from panel (b). (c) The
scaling behavior of the entropy with different interaction strengths:
$U_{0}=4.5\times 10^{-3}E_{R}$ (blue), $2.08\times 10^{-3}E_{R}$ (purple), and
$7.5\times 10^{-4}E_{R}$ (red). At $t=1500\,\mathcal{T}$, the sequence of
lines from top to bottom is blue, red, and purple. The scaled time $t_{s}$ is
calculated in such a way that $t_{s}\propto 1/\sqrt{U_{0}}$ and $t_{s}=t$ for
$U_{0}=7.5\times 10^{-4}E_{R}$. Figure 7: Evolution of the $\bf{p}$-space
distribution within and near the intermediate heating stage. At the beginning
of the intermediate ``heating" stage we see that the momentum- ($\bf{p}$-)
space distribution forms streaks in the horizontal direction; these clear up
to form sharp spots at later times when condensation into the Floquet
engineered Hofstadter BEC begins.  Figure 8: Time evolution of the vortex
distribution within and near the intermediate heating stage. The time periods
here are the same as in Fig. 7. This figure illustrates the evolution of the
wavefunction phase. The blue circles (red crosses) represent the vortices
(antivortices) located at positions consistent with predictions based on the
Floquet-Hofstadter ground state, while the filled cyan circles (purple stars)
indicate the vortex (antivortex) dislocations.

We define

$S_{\scriptscriptstyle{\mathrm{MBZ}}}=-\sum_{\bf{p}}\rho({\bf{p}})\ln\rho({\bf{p}})$
(3)

where $\rho({\bf{p}})=N_{\bf{p}}/N_{t}$ is the ratio of the particle number
$N_{\bf{p}}$ at momentum $\bf{p}$ to the total number $N_{t}$. Consistent with
the $\bf{p}$-space evolution in Fig. 6(a), we see that the entropy change is
also clearly divided into three stages, see Fig. 6(b). Soon after the initial
stage where the entropy is relatively stable, the entropy enters a `heating'
phase where it exhibits a rapid, exponential-like growth. This is associated
with a clear maximum in the entropy $S_{\scriptscriptstyle{\mathrm{MBZ}}}$. We
believe this rapid growth is rather generic, as we have seen it in simulating
other simpler Floquet systems, where it has been associated with an
``inflaton" picture Feng  _et al._ (2018). This picture inverts the usual
Bogoliubov description of the excitation spectrum, to describe a collection of
selectively amplified momentum modes which are at lower energy than the $p=0$
initial (unstable) state.

Interaction effects drive this behavior. The analysis of growth exponents in
simpler systems Feng  _et al._ (2018) suggests a universality where the
characteristic time scales in the second stage vary as $\propto
1/\sqrt{U_{0}}$. With this in mind, Fig. 6(c) shows the Floquet entropy
presented in terms of rescaled time variables for three different values of
the interaction energy $U_{0}$. This scaling with $\sqrt{U_{0}}$ provides an
adequate but imperfect fit to the power-law dependence in the scaling. Notably
this is consistent with an ``inflaton" model described elsewhere Feng  _et
al._ (2018).

We find the system appears to reach dynamical equilibrium in the third stage
near $t=6000\,\mathcal{T}$. The entropy value in this time domain is rather
stable; nevertheless, after $t=6000\,\mathcal{T}$,
$S_{\scriptscriptstyle{\mathrm{MBZ}}}$ begins to slowly increase. This can be
interpreted as heating in the long-time limit, which is also expected to occur
experimentally.

###  V.3 Analysis of the intermediate heating stage

We would like to draw particular attention to the intermediate ``heating"
stage we observe. This represents a crucial (albeit, transient) step in the
evolutionary dynamics in which there appears to be chaotic behavior, as seen
from Fig. 3. In this section we focus on this behavior by tracking specific
features in the evolution of the system through a sequence of figures. Whether
this chaotic state represents true ``turbulence" or not, it should be noted
that the GP dynamics is associated with weak quantum turbulent behavior
Tsatsos  _et al._ (2016) in non-equilibrated systems when a persistent source
of energy is applied, along with some degree of intrinsic or inevitable
dissipation and many-body interactions.

Fig. 7 illustrates how the momentum space distribution of the condensate
wavefunction evolves within and near the intermediate heating stage. The
dispersing or spreading out of the momentum peaks suggests highly chaotic
behavior which is observed over an extended time period. Here the
characteristic momentum peaks exhibit streaks along the horizontal direction,
beginning around $450\,\mathcal{T}$ and persisting for approximately another
$300\,\mathcal{T}$. After this, new momentum peaks associated with the new
Hofstadter condensate appear. It should be noted that the asymmetry between
the $x$ and $y$ directions which leads to the streaks, reflects the gauge used
to implement the artificial vector potential.

Even in our more detailed numerical studies App  in which high momentum states
are filtered out of the GP numerics, where we see a very ``purified" dynamical
evolution, we always find an intermediate chaotic heating stage. Strikingly
here the numerical filtering (associated with higher energy band occupation)
is able to remove most of the disorder from our momentum and real space plots,
except during this chaotic evolutionary stage.

Our numerical simulations enable us to more systematically investigate the
dynamics to determine how the system effects the transition from a
conventional condensate to one with the highly complex phase pattern of the
Hofstadter BEC. We saw in Fig. 4 that this introduction of phase is reflected
in a checkerboard vortex-antivortex pattern. We now exploit this pattern to
probe how phase coherence is dynamically established. This is illustrated in
Fig. 8 which indicates in cyan and purple how the vortex and anti-vortex
dislocations evolve against the background checkerboard pattern. The time
sequence is the same as for Fig. 7. What is striking is the suddeness of flux
penetration. These vortex-antivortex pairs with a large fraction of
dislocations initially appear precipitously at $450\,\mathcal{T}$. They then
rapidly reorganize as the dislocations are removed and the extended
checkerboard pattern is systematically developed. Interestingly, this time
frame where flux abruptly penetrates is roughly the same as the onset of
spreading out of sharp momentum peaks found in the initial condensate.

##  VI Conclusions

In conclusion, the work in this paper addresses the important and complicated
question of how one can successfully guide a wide class of Floquet engineered
systems Goldman and Dalibard (2014); Aidelsburger  _et al._ (2018); Wu  _et
al._ (2016) into the quantum regime. From our simulations, we show that a
dynamical conversion of a regular BEC into the Hofstadter ground state can be
realized with high efficiency. The dynamics involves an intriguing chaotic
``heating" stage during which ``magnetic" flux rapidly penetrates.

More generally, this paper addresses a need in the quantum gas community to
prepare novel quantum matter with Floquet engineering. Concerns about heating
are widespread not only for dynamical preparation of topological matter
444Alternative schemes, besides Floquet engineering, have recently been
discussed in the literature Sharma  _et al._ (2020); Wang  _et al._ (2021),
which may introduce some form of heating as well. but more generally to
surmount barriers McKay and DeMarco (2011) for reaching the quantum regime.
Quite intriguingly, the system is seen to overcome these challenges and the
way it does so is in many ways reminiscent of evolution in cosmological
models; this involves a similar time progression including an intermediate
turbulence Allahverdi  _et al._ (2010) en route to equilibration.

## Acknowledgements

We acknowledge L. Feng, L. W. Clark, Weihan Hsiao, B. M. Anderson, Nigel
Cooper, T. Bilitewski, I. S. Aronson, and A. Polkovnikov for helpful
discussions. We acknowledge support by the U. S. Department of Energy, Office
of Basic Energy Sciences, under contract number DE-SC0019216, the Army
Research Office under Grant No. W911NF-15-1-0113, and the University of
Chicago Materials Research Science and Engineering Center, funded by the
National Science Foundation under Grant No. DMR-1420709. A. Glatz acknowledges
support from the U. S. Department of Energy, Office of Science, Basic Energy
Sciences, Materials Sciences and Engineering Division. F. Setiawan
acknowledges support from the Army Research Office under Grant No.
W911NF-19-1-0328.

##  Appendix A Simulating the ideal Hofstadter model by Floquet-engineered
Hamiltonian

The general Harper-Hofstadter model we simulate with Floquet engineering is

$\displaystyle H_{iHH}=$
$\displaystyle-\sum_{i,j}\left(J_{x}^{\prime}e^{\imath\phi_{i,j}^{x}}\hat{a}^{\dagger}_{i+1,j}\hat{a}_{i,j}+J_{y}^{\prime}e^{\imath\phi^{y}_{i,j}}\hat{a}^{\dagger}_{i,j+1}\hat{a}_{i,j}+\mathrm{h.c.}\right),$
(A.1)

which can be approached using the Floquet Hamiltonian [with the lattice
potential given by Eq. (1) of the main text] in the ideal limit of
$J^{\prime}_{x}/\hbar\omega\rightarrow 0$ and
$J^{\prime}_{y}/\hbar\omega\rightarrow 0$. Here, the tunneling phases in $x$
and $y$ directions $\phi_{i,j}^{x}$ and $\phi_{i,j}^{y}$ are associated with
the applied vector field in each direction. The same filling factor can be
associated with different gauges Kennedy  _et al._ (2015). To calculate the
corresponding band structure in a chosen gauge (used in Ref. Bilitewski and
Cooper, 2015), we use the magnetic translation operator to identify the
eigenstates Aidelsburger (2014). The results are shown in Fig. 9, where panel
(a) shows the 3D band structure, and panel (b) shows the 2D color contour plot
of the lowest band. There are four different ground states in the lowest band.
This is to be contrasted with the unique ground state of the Floquet-
engineered Hamiltonian realized by implementing Eq. (1) with a moderately
large modulation frequency. The four-fold degeneracy in the ground states of
the ideal Hofstadter model is lifted in the Floquet-engineered case by higher-
order terms in $J^{\prime}_{x}/\hbar\omega$ and $J^{\prime}_{y}/\hbar\omega$
which hybridize states connected by wavevectors introduced by the oscillating
lattice.

Figure 9: Band structure and wavefunction for the ideal Hofstadter
Hamiltonian. (a) 3D plot of the band structure for the tunneling parameters
$J_{x}^{\prime}=J_{y}^{\prime}=1E_{R}$. (b) Color contour plot of the lowest
band. There are four different degenerate ground states labeled as min1 at
${\bf{k}}=(-\pi/2a,0)$, min2 at ${\bf{k}}=(0,-\pi/2a)$, min3 at
${\bf{k}}=(0,0)$, and min4 at ${\bf{k}}=(-\pi/2a,-\pi/2a)$.  Figure 10:
Ground-state wavefunction of the Floquet band. (a) Momentum-space expansion of
the ground state in terms of the reciprocal vectors
${\bf{G}}=n_{1}{\bf{G_{1}}}+n_{2}{\bf{G_{2}}}$. (b) Real-space representation
of the wavefunction showing both phase and amplitude. Note that both the phase
and amplitude distributions display a checkerboard pattern.

We now present more details on the band structure predicted by Floquet theory.
As shown by Fig. 1 in the main text, there is only one unique ground state
corresponding to two equivalent degenerate states in the lowest band. The
ground state can be expanded in terms of the total momentum
${\bf{p}}={\bf{k}}+{\bf{G}}$ at ${\bf{k}}=(-\frac{\pi}{2a},0)$. The
corresponding amplitude of the wavefunction at different $\bf{G}$ is presented
in Fig. 10(a) where ${\bf{G}}=n_{1}{\bf{G_{1}}}+n_{2}{\bf{G_{2}}}$ with
${\bf{G_{1,2}}}=(\frac{\pi}{2a},\pm\frac{\pi}{2a})$ being wavevectors of the
oscillating potential $V_{os}$ and $n_{1,2}$ being integers. The distribution
of the ground-state wavefunction has four dominant peaks in $\bf{G}$, which
correspond to the four characteristic $\bf{p}$ spots shown by Fig. 2 in the
main text. The relative phases of the ground-state wavefunction at the four
spots are
$\theta_{l}-\theta_{r}\approx\pi,\theta_{t}-\theta_{b}\approx\pi,\theta_{b}-\theta_{r}\approx
2.04$ radians, where the subscripts correspond to the left ($l$), right ($r$),
top ($t$), and bottom ($b$) momentum spots, respectively. We present an
overlay picture of the density and phase distributions in Fig. 10(b) for the
ground-state wavefunction in the Floquet-engineered case.

Figure 11: Effects of increasing the lattice depths (a) Lowest-energy band for
the ideal Hofstadter Hamiltonian calculated using the parameters in Fig. 9.
(b) Lowest-energy band for the Floquet Hamiltonian with very deep lattices.
Here $V_{y}=20E_{R},\,V_{x}=22E_{R},\,V_{yl}=1.0625E_{R}$,
$\hbar\omega=0.99664E_{R}$, and $\kappa=0.58\hbar\omega$. (c) Lowest Floquet
band with lattice depths used in the main text. The energy dispersion is shown
for the first magnetic Brillouin zone. And the energy is measured with respect
to the lowest energy $E_{min}$. The energy units (denoted by the arbitrary
unit [arb]) and $E_{min}$ are different for different cases.  Figure 12: Real-
space wavefunction in larger systems.  (a)Absolute value of density
correlation function $g_{2}^{n}(\Delta i,\Delta j)=\langle n(i,j)n(i+\Delta
i,j+\Delta j)\rangle/\langle n\rangle^{2}$ where $n(i,j)=|\psi(i,j)|^{2}$ is
the local density. Here, $\langle\cdots\rangle$ denotes averaging over
different ensembles and different $(i,j)$ positions with fixed relative
displacement $(\Delta i,\Delta j)$. (b)Absolute value of phase correlation
function $g_{1}$. The definition is given in the main text. (c) Distribution
of vortices (blue dots) and antivortices (red crosses). While the unique
ground state yields constant density correlation even at large distance, the
phase correlation reveals a finite size for the physically coherent regions.
The vortex structure is associated with vortex defects rather than
dislocations; this reflects the absence of vortices and antivortices. Figure
13: Illustration of flux penetration seen at early times within intermediate
heating stage.

By tuning the shaking frequency and the lattice depth in the Floquet
Hamiltonian, we can either go to the limit of an ideal Hofstadter band or stay
with the current Floquet band. In the following, we consider an intermediate
case between the two scenarios by proper choice of the parameters. We find an
intermediate state that exhibits two minima in the lowest band, see Fig. 11.
We can clearly see a trend of fewer minima in the lowest band with shallower
lattices and smaller shaking frequencies. One should note, however, though the
minima in Fig. 11(c) seem to be surviving minima in Fig. 11(a) at
${\bf{k}}=(\pm\frac{\pi}{2a},0)$, in fact, the ground states are quite
different for these two cases.

##  Appendix B Comparison between Floquet prediction and GP simulation results

###  B.1 Comparison of relative phases at the four characteristic $\bf{p}$
spots for the ground state

In the main text, we have already seen from the simulations that the particle
distribution is peaked at four characteristic $\bf{p}$ spots within the
magnetic Brillouin zone (see Fig. 3); this agrees well with the predictions
shown in Fig. 2(d), which are calculated by diagonalizing the Floquet
Hamiltonian. To be more quantitative, we look at the ratio between the
wavefunction amplitudes at the four spots. The time-averaged populations from
the simulations shown by Fig. 5(a) are consistent with our predictions
extracted from Fig. 10(a). The occupation at ${\bf{p}}=(0,\pm\frac{\pi}{2a})$
is $1.70$ times that at ${\bf{p}}=(\pm\frac{\pi}{2a},0)$, close to the
predicted value 1.67.

We can also check the relative phases between these characteristic spots. From
our simulations, we can directly obtain the wavefunction expansion in $\bf{p}$
space. We find that the relative phases are as follows:
$\langle\theta_{l}-\theta_{r}\rangle\approx
3.10,\langle\theta_{t}-\theta_{b}\rangle\approx
3.15,\langle\theta_{b}-\theta_{r}\rangle\approx 1.14$. These are consistent
with our prediction presented in Sec. I, except that there is approximately a
difference of order unity in $\langle\theta_{b}-\theta_{r}\rangle$. This
discrepancy, though quite robust and also present in our filtered simulations
(to be discussed in Appendix C. 2, where the BEC is much cleaner), seems to be
related to interaction effects, as when we decrease the interaction energy to
$2.5\times 10^{-4}E_{R}$, the phase difference increases significantly to
1.44.

###  B.2 Comparison of frequency spectrum for both the ground state and
excited states

In Fig. 5 of the main text, we have seen that besides the dominant ground
state in the observed BEC, there are also excited states. Table 1 presents a
summary of the energy comparisons which allow us to identify some of these
excited states. When we introduce a high-frequency filtration in our GP
dynamics (discussed in Appendix C. 2), we see a negligibly small occupation of
these higher bands.

Table 1: Identification of the ground state (min) and 3 excited states. The
table compares energies from simulations ($\varepsilon_{sim}$) and predictions
($\varepsilon_{pre}$) for the ground state and excited states in the Floquet-
Hofstadter band. As in the main text, min denotes the `ground state', while 1,
2, 3 refer to the same states as appear in Fig. 5.

$\begin{array}[]{|c|c|c|c|}\hline\cr\textrm{Spot index}&\textrm{Band
index}&\varepsilon_{sim}\,\textrm{modulo}\,\hbar\omega&\varepsilon_{pre}\,\textrm{modulo}\,\hbar\omega\\\
\hline\cr\textrm{min}&1&-0.011&-0.012\\\ \hline\cr 1&2&0.054&0.057\\\
\hline\cr 2&3&0.070&0.068\\\ \hline\cr 3&1&0.014&0.013\\\
\hline\cr\end{array}$

###  B.3 Comparison between the real-space correlation functions and vortex
structure

The BEC we obtain from simulations is a combination of mostly the ground state
and a small number of excited states, which is also reflected in the finite
coherence length in real space (see Fig. 4). To check this more thoroughly, we
look at an even bigger system ($4\times 4$ times larger). In Fig. 12, we show
the real-space density and phase correlation functions of the BEC obtained
from our GP simulations for such a system. Indeed, we see defects are present
leading to a finite range for the spatial phase coherence. They appear more
directly (mostly as vacancies) in the vortex-antivortex distribution.

It is important to note that the positions of the vortices and antivortices in
the distribution are not interchangeable as there is a single, non-degenerate
ground state. This can be viewed as a chiral asymmetry associated with an
artificial magnetic field. The defects we observe are manifested as an absence
of a vortex or antivortex as distinct from dislocations, and we believe they
derive, at least in part, from the excited states co-existing with the
condensate 555We note that, this alternating pattern is also found in the
initial magnetic Brillouin zone, but, importantly, it disappears during a
large portion of the evolution time, until well past the end of the inflation
stage.

##  Appendix C Dynamics of the Hofstadter BEC formation from GP simulations

###  C.1 Flux penetration dynamics

In Fig. 8 of the main text we discussed the rather precipitous appearance of
magnetic flux. Here in this appendix we look at this behavior in a more
refined way focusing on the transient period where the vortices first appear.
This is illustrated in Fig. 13 which shows how the system rapidly evolves from
one state with uniform phase to another with a complex phase pattern. The
color scheme is the same as in Fig. 8 and the cyan (purple) coloration labels
vortex (antivortex) dislocations. This figure illustrates the very transient
appearance of these dislocations. Their initial density is very high which
presumably represents the onset of flux penetration; it then rapidly decreases
as the checkerboard pattern of organized vortices associated with the
Hofstadter BEC begins to emerge. We speculate that this transient high density
of dislocations might indicate some degree of turbulent behavior which arises
due to the onset of an artificial vector potential.

###  C.2 Removing high-frequency contributions: purifying the BEC

Figure 14: Effects of high energy filtering in the simulations (a) Density
image in $\bf{p}$-space at different times. (b) Population curves, as defined
in Fig. 3. At $t=500\,\mathcal{T}$, the sequence of lines from top to bottom
is red, purple, blue, and green. (c) Entropy as a function of time. We find
that with a moderate filtration, the system reaches a cleaner BEC state with
the excitations almost completely gone. Importantly, the intermediate heating
stage is still present, and persists for a shorter period of time. The final
condensation fraction is greatly enhanced. Finally there is no sign of long
time heating as seen from the entropy which does not increase in the long-time
limit.

The dynamics we simulate in our GP equation do not involve energy dissipation,
once the initial state is established 666By starting from a random
configuration, the system evolves to the proper initial equilibrium state for
a negative phase angle $\gamma<0^{\circ}$ in Eq. (2), see Ref. Scherpelz  _et
al._ , 2014.. We stress that our simulations without dissipation are
consistent with experimental conditions in atomic systems. We might expect
that in actual experiments, the system picks up some higher momentum
($\bf{p}$) excitations over longer time scales. Indeed, these are presumably
responsible for the final heating stage. High $\bf{p}$ excitations do occur in
our simulations.

It is also informative, then, to compare the behavior of the system when these
high-energy excitations are removed, as is often done when studying the
stochastic GP equation Weiler  _et al._ (2008). For this reason, we apply a
high-momentum filter at each time step of the numerical integration, i.e.,
multiplying the Fourier transform of the order parameter by a Gaussian
function. Physically, this process may represent intrinsic losses, such as
those due to three-body and other collisions.

When the filtration is weak (where the momentum threshold above which the
modes will be removed is high), we find the behavior is generally unaffected
except that the higher-band excited states are no longer present, and the peak
structure of the target BEC becomes sharper as shown in Fig. 14.

This can be seen more clearly through a comparison with the unfiltered case.
Panels (a) and (c) in Fig. 14 can be contrasted with panels (a) and (b) in
Fig. 6, while Fig. 14(b) can be contrasted with Fig. 3(a). This comparison
reveals those features arising from higher-energy states, presumably deriving
from the role of the higher bands. For the most part the early-time evolution
is similar. The oscillations which are present without filtration are greatly
diminished, thus suggesting that these may come from higher band occupation.
We stress that the chaotic, intermediate heating stage is still present. We
also observe that the longer-time heating (seen in the entropy plot) vanishes
with filtration and consequently the BEC is more stable. Because the presence
of this longer-time heating appears more realistic, we conclude that the
unfiltered case is the more physical.

###  C.3 Role of Kibble-Zurek mechanism

Figure 15: Comparison of coherence between systems with different ramping
rates at $t=6000\,\mathcal{T}$. (a) Vortex (blue dot) and antivortex (red
cross) distribution of real space wavefunction. (b) Absolute value of phase
correlation function $g_{1}$ for real space wavefunction. The ramp time is
$300\,\mathcal{T}$ (left panel) and $30\,\mathcal{T}$ (right panel),
respectively.

Here, we want to briefly discuss the role of the Kibble-Zurek (KZ) mechanism.
When a dynamical system crosses a critical point like a phase transition point
by ramping a key parameter such as the shaking amplitude, the correlation
length in the ordered phase is determined by the ramping rate. The slower the
rate is, the bigger the domain size or the correlation length. Here, one might
wonder if similar effects are responsible for defects observed in this paper.
We have a similar transition point, but the change of the band structure and
thus the ground-state wavefunction is abrupt (and more first-order like) at
the critical shaking amplitude $\kappa=0$.

Our results suggest an absence of important effects associated with the KZ
mechanism. This can be seen by comparing the real-space distribution of the
vortices between cases with different ramping rates, see Fig. 15. While the
right panel in Fig. 15 corresponds to a very fast ramp (essentially a quench),
the left panel corresponds to a slow process (with a ramping period 10 times
as long). The domain size and phase correlation length are found to be
comparable, indicating the lack of important KZ effects here.

Figure 16: Evolution of the $p$-space distribution for stronger inter-particle
interactions: $U_{0}=0.009E_{R}$. Here we increase the repulsion by a factor
of 10 compared to that used in the main text. We see that the system seems to
never reach the Hofstadter BEC and remains in a chaotic state throughout the
simulation.

###  C.4 Strong interaction effect: absence of condensation

It is important to investigate the effects associated with the interaction
strength $U_{0}$, since the validity of Floquet predictions is based on
assuming that such interactions are negligible. GP simulations allow the
simultaneous incorporation of Floquet engineering along with interaction
effects. Our results show that with a moderately large $U_{0}$ the
evolutionary behavior tends to be very noisy without clearing up, see Fig. 16.
This behavior suggests the failure to form a BEC. Indeed this is consistent
with observations in Ref. Miyake  _et al._ , 2013, where, when the collision
rate is too high, this is seen to seriously disturb the single-particle band
structure.

## References

  * Jotzu  _et al._ (2014) G. Jotzu _et al._ , Nature 515, 237 (2014).
  * Aidelsburger  _et al._ (2014) M. Aidelsburger _et al._ , Nature Physics 11, 162 (2014).
  * Miyake  _et al._ (2013) H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 185302 (2013).
  * Lin  _et al._ (2009) Y.-J. Lin _et al._ , Nature 462, 628 (2009).
  * Eckardt (2017) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Reports on Progress in Physics 77, 126401 (2014).
  * Galitski  _et al._ (2019) V. Galitski _et al._ , Physics Today 72(1), 38 (2019).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Struck  _et al._ (2013) J. Struck _et al._ , Nature Physics 9, 738 (2013).
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard, Phys. Rev. X 4, 031027 (2014).
  * Dalibard (2015) J. Dalibard, arXiv:1504.05520 (2015).
  * Aidelsburger  _et al._ (2018) M. Aidelsburger, S. Nascimbene, and N. Goldman, Comptes Rendus Physique 19, 394 (2018), quantum simulation / Simulation quantique.
  * Barbiero  _et al._ (2019) L. Barbiero, C. Schweizer, M. Aidelsburger, E. Demler, N. Goldman, and F. Grusdt, Science Advances 5: eaav7444 (2019).
  * Thouless  _et al._ (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Quelle  _et al._ (2017) A. Quelle, C. Weitenberg, K. Sengstock, and C. M. Smith, New Journal of Physics 19, 113010 (2017).
  * Struck  _et al._ (2012) J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Phys. Rev. Lett. 108, 225304 (2012).
  * Bukov  _et al._ (2015) M. Bukov, L. D'Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015).
  * Bukov and Polkovnikov (2014) M. Bukov and A. Polkovnikov, Phys. Rev. A 90, 043613 (2014).
  * Bilitewski and Cooper (2015) T. Bilitewski and N. R. Cooper, Phys. Rev. A 91, 063611 (2015).
  * Rubio-Abadal  _et al._ (2020) A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Phys. Rev. X 10, 021044 (2020).
  * Kennedy  _et al._ (2015) C. J. Kennedy, W. C. Burton, W. C. Chung, and W. Ketterle, Nature Physics 11, 859 (2015).
  * Oka and Kitamura (2019) T. Oka and S. Kitamura, Annual Review of Condensed Matter Physics 10, 387 (2019).
  * Hunt  _et al._ (2013) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Science 340, 1427 (2013).
  * Dean  _et al._ (2013) C. R. Dean _et al._ , Nature 497, 598 (2013).
  * Ponomarenko  _et al._ (2013) L. A. Ponomarenko _et al._ , Nature 497, 594 (2013).
  * Owens  _et al._ (2018) C. Owens, A. LaChapelle, B. Saxberg, B. M. Anderson, R. Ma, J. Simon, and D. I. Schuster, Phys. Rev. A 97, 013818 (2018).
  * Allahverdi  _et al._ (2010) R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine, and A. Mazumdar, Annual Review of Nuclear and Particle Science 60, 27 (2010).
  * (31) See appendix.
  * Price  _et al._ (2016) R. M. Price, D. Trypogeorgos, D. L. Campbell, A. Putra, A. Valdés-Curiel, and I. B. Spielman, New Journal of Physics 18, 113009 (2016).
  * Aidelsburger (2014) M. Aidelsburger,  _Artificial gauge fields with ultracold atoms in optical lattices_ , Ph.D. thesis, Ludwig-Maximilians-Universit$\ddot{\textrm{a}}$t M$\ddot{\textrm{u}}$nchen (2014).
  * Bilitewski (2016) T. Bilitewski,  _Interacting atoms in time-dependent potentials and artificial gauge fields_ , Ph.D. thesis, University of Cambridge (2016).
  * Note (1) The Floquet-Hofstadter ground state is defined as the lowest energy state in the ground band. This is the band in which the ${\bf{k}}=0$ state is adiabatically connected to the ground state in the static system.
  * Pathria and Morris (1990) D. Pathria and J. Morris, Journal of Computational Physics 87, 108 (1990).
  * Micha and Tkachev (2003) R. Micha and I. I. Tkachev, Phys. Rev. Lett. 90, 121301 (2003).
  * Micha and Tkachev (2004) R. Micha and I. I. Tkachev, Phys. Rev. D 70, 043538 (2004).
  * Note (2) These vortices are intrinsic and should be contrasted with vortex patterns associated with a rotating trap as in Phys. Rev. Lett. 92, 020403\.
  * Note (3) The time slice within each cycle is chosen to be at zero phase such that $t$ in Eq. (1) is $0+n\mathcal{T}$ where $n$ are integer multiples.
  * Feng  _et al._ (2018) L. Feng, L. W. Clark, A. Gaj, and C. Chin, Nature Physics 14, 269 (2018).
  * Tsatsos  _et al._ (2016) M. C. Tsatsos, P. E. Tavares, A. Cidrim, A. R. Fritsch, M. A. Caracanhas, F. E. A. dos Santos, C. F. Barenghi, and V. S. Bagnato, Physics Reports 622, 1 (2016).
  * Wu  _et al._ (2016) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Science 354, 83 (2016).
  * Note (4) Alternative schemes, besides Floquet engineering, have recently been discussed in the literature Sharma  _et al._ (2020); Wang  _et al._ (2021), which may introduce some form of heating as well.
  * McKay and DeMarco (2011) D. McKay and B. DeMarco, Reports on Progress in Physics 74, 054401 (2011).
  * Note (5) We note that, this alternating pattern is also found in the initial magnetic Brillouin zone, but, importantly, it disappears during a large portion of the evolution time, until well past the end of the inflation stage.
  * Note (6) By starting from a random configuration, the system evolves to the proper initial equilibrium state for a negative phase angle $\gamma<0^{\circ}$ in Eq. (2), see Ref. Kathy.
  * Weiler  _et al._ (2008) C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, and B. P. Anderson, Nature 455, 948 (2008).
  * Sharma  _et al._ (2020) V. Sharma, S. Choudhury, and E. J. Mueller, Phys. Rev. A 101, 033609 (2020).
  * Wang  _et al._ (2021) X.-Q. Wang, G.-Q. Luo, J.-Y. Liu, W. V. Liu, A. Hemmerich, and Z.-F. Xu, Nature 596, 227 (2021).
  * Scherpelz  _et al._ (2014) P. Scherpelz, K. Padavić, A. Rançon, A. Glatz, I. S. Aranson, and K. Levin, Phys. Rev. Lett. 113, 125301 (2014).