# Lightweight Diffusion Models with Distillation-Based Block Neural
Architecture Search

Siao Tang
Tsinghua University
Xin Wang
Tsinghua University
Corresponding author Hong Chen
Tsinghua University
Chaoyu Guan
Tsinghua University
Yansong Tang
Tsinghua University
Wenwu Zhu *
Tsinghua University

###### Abstract

Diffusion models have recently shown remarkable generation ability, achieving
state-of-the-art performance in many tasks. However, the high computational
cost is still a troubling problem for diffusion models. To tackle this
problem, we propose to automatically remove the structural redundancy in
diffusion models with our proposed Diffusion Distillation-based Block-wise
Neural Architecture Search (DiffNAS). Specifically, given a larger pretrained
teacher, we leverage DiffNAS to search for the smallest architecture which can
achieve on-par or even better performance than the teacher. Considering
current diffusion models are based on UNet which naturally has a block-wise
structure, we perform neural architecture search independently in each block,
which largely reduces the search space. Different from previous block-wise NAS
methods, DiffNAS contains a block-wise local search strategy and a retraining
strategy with a joint dynamic loss. Concretely, during the search process, we
block-wisely select the best subnet to avoid the unfairness brought by the
global search strategy used in previous works. When retraining the searched
architecture, we adopt a dynamic joint loss to maintain the consistency
between supernet training and subnet retraining, which also provides
informative objectives for each block and shortens the paths of gradient
propagation. We demonstrate this joint loss can effectively improve model
performance. We also prove the necessity of the dynamic adjustment of this
loss. The experiments show that our method can achieve significant
computational reduction, especially on latent diffusion models with about 50%
MACs and Parameter reduction.

##  1 Introduction

The diffusion model is one of the generative models which learns the
underlying data distribution by adding noise and iteratively denoising [8, 31,
18]. Diffusion models have shown remarkable generative ability, outperforming
advanced GANs and VAEs on a series of benchmarks [10, 4]. Diffusion models
have been applied to a variety of tasks and achieved great success, such as
text-to-image generation [26], text-to-video generation [37], graph generation
[34], and speech synthesis [2].

However, a critical limitation of diffusion models is their high computational
cost, which is mainly from two aspects. On the one hand, diffusion models need
a long sequence of noise prediction steps (e.g., 1000 steps for DDPM [8]) to
generate a batch of high-quality outputs. On the other hand, the computation
amount of each noise prediction is large because of the high computational
overhead of the network architecture. To tackle the first problem, substantial
works [18, 20, 28, 21, 35, 22] engage in reducing sampling steps by finding
more effective sampling trajectories or knowledge distillation.

Figure 1:  Concept comparison between DiffNAS and previous computation
reduction works. Previous works mainly focus on efficient sampling or
quantization. From another important perspective, we propose to remove the
structural redundancy of diffusion models by automatically searching for the
most efficient architecture.

We focus on the second problem, i.e., the high computational overhead brought
by large network architecture. Quantization [29, 16] and pruning methods [5]
have been explored to compress diffusion models. From an orthogonal
perspective, we engage in removing the structural redundancy for diffusion
models by automatically searching for the most efficient architecture. Figure
1 shows the concept comparison between our method and previous works.
Specifically, we propose a Diffusion Distillation-based Block-wise Neural
Architecture Search (DiffNAS) framework, which searches for the smallest
architecture whose generative quality is better or on par with a pretrained
teacher, through the proposed block-wise local search strategy during NAS
evaluation and the proposed dynamic joint loss during subnet retraining. Our
method is orthogonal to quantization or pruning methods and can be combined
with them to achieve further compression. Furthermore, we use the performance
of the baseline model (that needs to be compressed) as the constraint for our
method. This constraint guarantees the obtained architecture can theoretically
have better performance than the baseline, which previous quantization or
pruning methods can’t ensure.

Concretely, considering the pyramidal block-wise architecture of U-Net [27]
that most diffusion models [8, 31, 26, 32, 24, 33] are based on, we adopt the
block-wise NAS paradigm [13, 14] and regard each resolution level as a block.
We use the baseline model expected to be compressed as the teacher. We use the
intermediate distillation loss, i.e., $L2$ loss, to supervise the block-wise
training of the supernet. The loss measures the distance between the teacher
output and the student output at each block.

During the NAS evaluation, we use the distillation loss to measure the
distance of subnets to the teacher, since a smaller loss means a better
ability to imitate the teacher. Previous works [13] use a global search
strategy to search for the expected architecture under certain constraints.
However, we demonstrate the global search strategy leads to unfairness across
different blocks, because it fails to balance the loss value magnitude of
these blocks. Therefore, we propose to block-wisely select the smallest
subnets which have a lower evaluation loss compared to the subnet having the
same architecture as the teacher. This loss constraint theoretically
guarantees the searched subnet can have better performance than the teacher.
Finally, we obtain the complete subnet by concatenating the searched
architectures of each block.

During the retraining, we observe that retraining the searched subnet only
with the original noise prediction loss leads to poor performance than the
teacher. One main reason is that only using the original loss can not maintain
the consistency between supernet training and subnet retraining. Specifically,
the distillation loss actually reflects the ability of these subnets to
imitate the teacher when training supernet and evaluating. This indicates we
also need the teacher to guide the retraining process of subnets. Therefore,
instead of retraining the subnet without distillation as previous works [13],
we propose to train with a dynamic joint loss combining the distillation loss
and the original noise prediction loss. Specifically, we gradually increase
the weight of the noise prediction loss and decrease the weight of the
distillation loss. The early retraining can be regarded as a warm-up stage
when the student quickly learns the knowledge from the teacher. In the later
retraining, the student needs to be refined, so the noise prediction loss
becomes more important since it directly determines the accuracy of noise
prediction. The distillation loss maintains the consistency between supernet
training and retraining. Moreover, it further improves model performance since
the intermediate feature maps provide information at different resolution
levels and shorten the paths of gradient propagation. We conduct experiments
on various datasets with different resolutions and computation amounts.

In conclusion, our contributions are summarized as follows:

  1. 1.

To compress diffusion models, we propose to remove the structural redundancy
by automatically searching for the most efficient architecture.

  2. 2.

We propose DiffNAS, a framework to compress diffusion models without degrading
performance. DiffNAS contains a novel block-wise search strategy to avoid
unfairness brought by global search, and a novel dynamic joint loss for NAS
retraining to maintain the consistency between supernet training and subnet
retraining.

  3. 3.

The results show that our method achieves about 30% MACs reduction on pixel-
wise models, as well as about 50% MACs and Parameters reduction on latent
diffusion models, while maintaining on-par or even achieving better
performance.

##  2 Related Work

Neural Architecture Search (NAS). Neural Architecture Search (NAS) aims to
explore more efficient architectures than manually designed ones, replacing
the efforts of human experts to design models. NAS has three essentials, i.e.,
search space, search strategy, and evaluation strategy. Early works [41, 42,
1] sample an architecture candidate by RNN or evolutionary algorithm and
evaluate its performance after training from scratch. The type of NAS costs a
large amount of computation and is thus impractical. Afterward, most works
employ a weight-sharing supernet to represent search space. Gradient-based
methods [36, 17] assign choosing weights to architecture options and optimize
these weights with the net parameters jointly. One-shot methods disentangle
supernet training with architecture search, usually training the supernet with
path sampling or path dropout. when evaluating, subnets directly inherit
parameters from the trained supernet. To reduce search space and improve the
accuracy of model ranking, block-wise NAS methods are proposed [13, 14]. They
use a distance loss between the feature maps of the teacher and the student to
guide block-wise supernet training and evaluation. In our work, we simulate
the same supernet training process as previous works [13]. However, to cope
with the unfairness and the inconsistency problems, our method proposes a
novel block-wise search strategy and a retraining strategy with a dynamic
joint loss.

Model Compression. Model compression is a significant technology to lighten
large neural networks. Common methods include Network Pruning, Knowledge
Distillation (KD), and Quantization. Pruning aims to remove the relatively
less important parameters. Knowledge distillation aims to leverage the
knowledge of a teacher network to guide the training of a more compact student
network. Quantization refers to the process of using fewer bits to represent
float point tensors and perform model inference. In practice, pruning might be
limited by the support of inference platforms and hardware. Quantization can
not cope with the redundancy from architecture. Different from the perspective
of using few bits or cutting off unimportant weights, our method engages in
searching for appropriately small architectures. Therefore, our method is
orthogonal to quantization and pruning, which means it can also be combined
with pruning or quantization to further compress diffusion models.

##  3 Preliminaries

Diffusion Models. Diffusion models are probabilistic generative models that
learn the underlying data distribution by adding noise to the input and
gradually denoising the noisy variable. The initial idea of diffusion models
can trace back to [30]. In recent years, diffusion models have shown powerful
performance in generative tasks and attracted a significant amount of
attention. Besides noise prediction, diffusion models have several variants,
e.g., data prediction [25] and velocity prediction [9] models. In this paper,
we discuss the most frequently-used form, i.e., noise prediction models.

Diffusion models define a forward process that adds noise to the original
input:

$\displaystyle
q\left(\boldsymbol{x}_{t}|\boldsymbol{x}_{0}\right)=\mathcal{N}\left(\alpha_{t}\boldsymbol{x}_{0},\sigma^{2}_{t}\boldsymbol{I}\right),$

where $\alpha_{t}$ and $\sigma_{t}$ are differentiable functions of $t$, which
are determined by certain noise schedule. The formula can be reparametered as
$\boldsymbol{x}_{t}=\alpha_{t}\boldsymbol{x}_{0}+\sigma_{t}\boldsymbol{\epsilon}$,
where $\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. The
signal-to-noise-ratio (SNR) $\alpha_{t}/\sigma_{t}$ is strictly decreasing
with respect to $t$, such that the $q(x_{t}|x_{0})$ is nearly a unit Gaussian
when $t$ is large enough. The neural network $\boldsymbol{\epsilon_{\theta}}$
is trained to predict the true noise $\epsilon$ given $\boldsymbol{x}_{t}$ and
$t$ , i.e., minimizing:

$\displaystyle\mathbf{E}_{\boldsymbol{x}_{0},\boldsymbol{\epsilon},t}[w_{t}\left\|\boldsymbol{\epsilon}_{\theta}\left(\alpha_{t}\boldsymbol{x}_{0}+\sigma_{t}\boldsymbol{\epsilon},t\right)-\boldsymbol{\epsilon}\right\|_{2}^{2}].$

At inference time, diffusion models generate new samples by iteratively
denoising $\boldsymbol{x_{T}}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ using
some kind of sampling strategies, e.g., DDIM [31], PNDM [18] and DPM-Solver
[20].

Block-wise NAS. Block-wise NAS is firstly proposed by Liu et al. [13]. The
motivation is to reduce search space and promote the accuracy of model ranking
by narrowing the performance gap between proxy subnets with shared weights and
the retrained stand-alone ones. The gap means the ranks of subnets inheriting
weights from the supernet are not consistent with the ranks of the retrained
ones. Some works [3, 15] show the gap is narrowed as the search space
decreases. Let $N$ denote the whole supernet, block-wise NAS divides $N$ into
$n$ blocks:

$\displaystyle\mathcal{N}=\mathcal{N}_{n-1}\ldots\mathcal{N}_{i+1}\circ\mathcal{N}_{i}\cdots\circ\mathcal{N}_{0}.$

Then, block-wise NAS train each block of the supernet separately:

$\displaystyle\mathcal{W}_{i}^{*}=$
$\displaystyle\min_{\mathcal{W}_{i}}\mathcal{L}_{\text{train
}}\left(\mathcal{W}_{i},\mathcal{A}_{i};\mathbf{X_{i}},\mathbf{Y_{i}}\right),\quad
i=0,1\cdots n-1,$

where $\mathcal{A}_{i}$ denotes the search space of the $i$-th block and
$\mathcal{W}_{i}$ denotes the parameters. $\mathbf{X_{i}}$ denotes the input
of the $i$-th block, which is the output of the $(i-1)$-th block of the
teacher. $\mathbf{Y_{i}}$ denotes the groundtruth of the $i$-th block, which
is the output of the $i$-th block of the teacher. Block-wise NAS can extremely
reduce the search space, from $\prod_{i=0}^{n-1}C^{d_{i}}$ to $C^{d_{i}}$,
where $C$ denotes the number of the optional operations in each layer and
$d_{i}$ denotes the number of layers in the $i$-th block.

##  4 Methodology

Figure 2 illustrates the whole framework of our method. Given a target
pretrained model that is expected to be lightened, we first block-wisely train
a supernet containing a weight-sharing search space, with distillation as the
supervision. The search space includes the architecture of the teacher. Then
we block-wisely search for the smallest subnet with better performance than
the teacher. Finally, we retrain the subnet with a dynamic joint loss
combining the knowledge distillation loss and the original loss. In this
section, we will describe the three stages in detail. The time and memory
costs of our method are not more than training a baseline, which we show
details in the Appendix (Sec. 8).

###  4.1 Block-wise NAS Training

Most diffusion models use UNet-based architecture which has several levels of
resolutions for feature encoding and decoding. The encoder extracts features
at different abstraction levels, from low-level textures to high-level
semantics, while the decoder gradually recovers the information at different
levels. Considering this naturally block-wise architecture, we divide the
whole UNet into blocks, regarding each resolution level as a block. In other
words, we perform NAS independently for each block. we use the $L2$ loss
between the block output of teacher and student to supervise NAS training for
each block. Let $T=\\{E^{t}_{i},D^{t}_{i}\\}$ and $\
S=\\{E^{s}_{i},D^{s}_{i}\\},i=0,1,...,n-1$ respectively denote the teacher and
the student with $n$ resolution levels, where $E_{i}$ and $D_{i}$ denote the
encoder and the decoder at $i$-th resolution level. For simplicity, we ignore
the notation of the middle block that is between the encoding stage and the
decoding stage. The training objective for each block is written below:

$\displaystyle\mathcal{W}_{i}^{*}=\min_{\mathcal{W}_{i}}$
$\displaystyle\mathcal{L}_{\text{train
}}\left(\mathcal{W}_{i},\mathcal{A}_{i};\mathbf{X_{i}},\mathbf{Y_{i}}\right),\quad
i=0,1\cdots 2n-1,$
$\displaystyle\mathcal{L}_{\text{train}}=||\text{Block}_{\alpha}(\mathbf{X_{i}})-\mathbf{Y_{i}}||^{2}_{2}\
,$
$\displaystyle\mathbf{X_{i}}=\begin{cases}\mathbf{Y_{i-1}},\text{Block}_{i}\in\\{E_{i}^{s}\\}\\\
\text{Cat}(\mathbf{Y_{i-1}},\mathbf{Skip_{i}}),\text{Block}_{i}\in\\{D_{i}^{s}\\},\end{cases}$

where $\mathcal{W}_{i}$ denotes the net parameters, $\mathcal{A}_{i}$ denotes
the architecture search space of $i$-th block, and $\alpha$ denotes the
selected subnet. $\mathbf{Y_{i}}$ is the ground-truth feature map which is the
output of the $i$-th block of the teacher. $\mathbf{X_{i}}$ denotes the input
of the $i$-th block, which comes from the teacher. When it is an encoder
block, $\mathbf{X_{i}}$ is the teacher output of the $(i-1)$-th block, or else
it is the concatenation of the skip features and the teacher output of the
$(i-1)$-th block. $\mathcal{L}_{\text{train}}$ is a L2 loss with mean
reduction.

The whole supernet training procedure is shown in the left part of Figure 2.
Each block represents a resolution level. We use the random single-path
strategy [6] for subnet sampling. It means that for each training step, we
randomly select one operation at each layer and thus obtain a random
architecture to be trained. Note that for diffusion models, we observe that
the $L2$ loss remains almost unchanged after training some epochs, even though
further training continues to improve the generative performance. Therefore,
we only need much fewer epochs for NAS supernet training compared to training
a standalone diffusion model from scratch.

Figure 2: The illustration of our method DiffNAS, which consists of three
stages, i.e., supernet training, block-wise search, and subnet retraining. In
the supernet training, different blocks represent different levels of
resolution. For each block, the intermediate feature maps of the teacher are
used to supervise the training of the supernet. The input features of each
student block also come from the previous outputs of the teacher. The mark ’+’
denotes the concatenation operation. Then we block-wisely search for the
smallest subnet with better performance than the teacher, and concatenate the
searched subnet in each block as the final net. Finally, we retrain the
searched subnet with a dynamic joint loss. The black line denotes the path of
the distillation loss, while the blue dashed line denotes the path of the
original loss.

###  4.2 Block-wisely Search for the Smallest Student

After supernet training, we conduct the evaluation on the training dataset. We
search for the model that has the least computation cost with on-par or even
better performance than the teacher. Previous works [13] use a global strategy
to search for the appropriate subnets. Since different blocks have different
scales of feature values, they use a relative loss such as the relative $L1$
loss to normalize the loss values of different blocks. The relative $L1$ loss
is to divide the original $L1$ loss by the standard deviation among all
elements of the ground-truth feature map. Then they use the summation of these
losses of each block as the final loss to evaluate the subnet performance.
However, we observe that the relative loss also leads to unfairness, because
some blocks have much larger losses than others. It means that blocks with
much larger losses dominate the final loss and thus dominate the search
process. In this case, the summed loss is meaningless. Therefore, instead of
the global search requiring the balance of losses, we propose to block-wisely
select the subnet. Specifically, we use the subnet which has the same
architecture as the teacher as the baseline. Let $\alpha_{base}$ denote it.
For each block, we select the smallest subnet that has a better evaluation
loss compared to $\alpha_{base}$. Finally, we obtain the complete subnet by
concatenating the searched architectures of each block. For each block, the
search strategy is shown as follows:

$\displaystyle\alpha_{i}^{*}=\underset{\alpha_{i}\in\mathcal{A}_{i}}{\arg\min}\
\text{Cost}(\alpha_{i}),\ i=0,1\cdots 2n-1,$ $\displaystyle\ s.t.\
\mathcal{L}_{\text{val}}\left(w_{i},\alpha_{i};\mathbf{X_{i}},\mathbf{Y_{i}}\right)\leq
r\mathcal{L}_{\text{val}}\left(w_{i},\alpha_{base};\mathbf{X_{i}},\mathbf{Y_{i}}\right),$
$\displaystyle\qquad\text{where}\quad\mathcal{L}_{\text{val}}=||\text{Block}_{\alpha}(\mathbf{X_{i}})-\mathbf{Y_{i}}||^{2}_{2}\
,$

where Cost(.) represents the computation cost function with respect to the
architecture and $\alpha_{i}$ denotes the architecture candidate of the $i$-th
block. $\mathcal{L}_{\text{val}}$ denotes the validation loss. $r$ is a
coefficient that can relax the constraint if we set it to a value larger than
1. Smaller subnets tend to have larger loss values, so we can increase $r$
slightly to obtain a sequence of smaller subnets. In this way, we can control
the compression rate. Note that in our experiments, after retraining with our
proposed strategy, the subnet with a $r$ slightly greater than 1 (e.g., 1.02)
can also have better performance than the teacher. The final optimal subnet is
the concatenation of $\alpha_{i}^{*}$ of all blocks. Benefiting from the
reduced search space brought by block-wise NAS, we can exhaustively traverse
and evaluate all subnets in each block. To reduce the time costs, the
evaluation for each subnet only computes several batches rather than the whole
dataset, which largely speeds up the evaluation process. The experiments
indicate this partial evaluation can achieve good enough results. The
Algorithm 1 illustrates the whole search process, where ChooseBase denotes the
function that chooses the subnet with the same architecture as the teacher.

Input: Supernet, Cost function, block input $\mathbf{X}$, ground-truth
$\mathbf{Y}$, release coefficient $r$

Output: The best lightweight model

Define SearchBlock(i):

BaseNet = ChooseBase(Supernet[i])

$\text{loss}_{\text{base}}$ =
$||\text{BaseNet}(\mathbf{X_{i}})-\mathbf{Y_{i}}||_{2}^{2}$

$\text{loss}_{\text{min}}$ =$\text{cost}_{\text{min}}$ = INF

$\text{net}_{\text{optimal}}$ = None

for  _Subnet in Supernet[i]_ do

$\mathbf{\hat{Y}_{i}}$ = Subnet($\mathbf{X_{i}}$)

loss = $||\mathbf{\hat{Y}_{i}}-\mathbf{Y_{i}}||^{2}_{2}$

cost = Cost(Subnet)

if  _$\text{loss}\leq r\text{loss}_{\text{base}}$_ then

if  _$cost <{cost}_{min}$_ then

$\text{cost}_{\text{min}}$ = cost

$\text{loss}_{\text{min}}$ = loss

$\text{net}_{\text{optimal}}$= Subnet

else if  _$cost=={cost}_{min}$_ then

if  _$loss <{loss}_{min}$_ then

$\text{loss}_{\text{min}}$ = loss

$\text{net}_{\text{optimal}}$= Subnet

end if

end if

end for

Subnet = []

for  _$i\ in\ range(block\\_num)$_ do

Subnet.append(SearchBlock(i))

end for

Subnet = Concatenate(Subnet)

Return Subnet

Algorithm 1 Search for the smallest Subnet.

###  4.3 Subnet Retraining with Joint Dynamic Loss

When retraining the searched subnet from scratch, we observe that only using
the noise prediction leads to poor performance. In supernet training and
evaluation, We use the distance between the student and the teacher to measure
the quality of subnets. In other words, we actually measure the ability to
imitate the teacher. Therefore, retraining the selected best student with the
guidance of the teacher is important, which maintains the consistency between
subnet evaluation and retraining. We propose a dynamic joint loss combining
the distillation loss and the noise prediction loss, formulated as follows:

$\displaystyle\mathcal{L}_{\text{retrain}}=\gamma\cdot(1-\beta)\mathcal{L}_{\text{dis}}+\beta\mathcal{L}_{\text{ori}}\
,$
$\displaystyle\mathcal{L}_{\text{ori}}=||\epsilon_{\theta}\left(\alpha_{t}x+\sigma_{t}\epsilon,t\right)-\epsilon||_{2}^{2}\
,$
$\displaystyle\mathcal{L}_{\text{dis}}=\sum_{i}\mathcal{L}_{\text{train}}=\sum_{i}||\text{Block}_{i}(\mathbf{X_{i}})-\mathbf{Y_{i}}||^{2}_{2}\
,$

where $\beta$ is the weight that balances the distillation loss and the noise
prediction loss. $\beta$ is set to 0 at the beginning and finally equals to 1.
We use a linear or a step function to adjust $\beta$. The experiments show
both functions are effective, and we recommend the step function as the first
choice to save time. $\mathcal{L}_{\text{ori}}$ is computed by passing
$\alpha_{t}x+\sigma_{t}\epsilon$ and $t$ through all the layers, while
$\mathcal{L}_{\text{dis}}$ is computed block-wisely whose inputs come from the
teacher just like the NAS training. This means in the initial training, we
mainly use the guidance of the teacher. The guidance can provide more
informative objectives, since the intermediate feature maps provide
information on different abstraction levels and shorten the paths of gradient
propagation. This can be regarded as a warm-up stage where the student quickly
learns the knowledge from the teacher. In the later training, the student
needs to be refined, so $\mathcal{L}_{ori}$ is more important since it
directly determines the accuracy of noise prediction. We also use a parameter
$\gamma$ to balance the learning rates of $\mathcal{L}_{\text{dis}}$ and
$\mathcal{L}_{\text{ori}}$.

Model | MACs (G) | ddim-50 | ddim-100 | dpm-20 | dpm-25
---|---|---|---|---|---
Teacher | 6.05 | 4.75 | 4.25 | 6.41 | 4.18
S1(linear) | 4.18(-31%) | 4.66 | 4.32 | 5.08 | 4.12
S2(linear) | 4.51(-26%) | 4.56 | 4.17 | 4.69 | 3.92
S2(step) | 4.51(-26%) | 4.59 | 4.20 | 5.29 | 3.97
Random(linear) | 4.527(-26%) | 5.41 | 4.85 | 4.62 | 3.95
Random(step) | 4.468(-26%) | 4.69 | 4.33 | 5.65 | 4.21
No-Dis-S2 | 4.510(-26%) | 5.24 | 4.81 | 5.10 | 4.02
Fixed-Loss-S1 | 4.175(-31%) | 9.56 | - | - | 11.97

Table 1: CIFAR-10 image generation with pixel-wise model.

Model | MACs (G) | dpm-20 | dpm-25 | pndm-50
---|---|---|---|---
Teacher | 23.88 | 4.46 | 4.33 | 5.66
S1(linear) | 17.98(-25%) | 4.16 | 4.34 | 5.47
S1(step) | 17.98(- 25%) | 3.73 | 4.03 | 5.53

Table 2: CelebA image generation with pixel-wise model.

Model | MACs (G) | dpm-20 | dpm-25 | pndm-50
---|---|---|---|---
Teacher | 62.05 | 13.61 | 10.82 | 12.86
S1(step) | 46.38(-25%) | 12.34 | 10.51 | 13.93

Table 3: LSUN-church 128x128 image generation with pixel-wise model.

Dataset | Model | MACs (G) | Params (M) | FID | sFID
---|---|---|---|---|---
LSUN-church | Teacher | 18.64 | 294 | 5.26 | 24.85
| S1(step) | 10.12(-46%) | 123(-58%) | 5.75 | 25.34
CelebA-HQ | Teacher | 96.02 | 273 | 7.85 | 16.75
| S1(step) | 53.06(-45%) | 149(-46%) | 8.49 | 17.97

Table 4: LSUN-church 256x256 and CelebA-HQ 256x256 generation with latent
diffusion models. Use DDIM sampler of 50 steps.

##  5 Experiments

###  5.1 Implementation Details

Search space. For simplicity, we only search the convolution kernel size of
Resnet blocks which accounts for the majority of the total computation. Our
method can also be suitable for larger search spaces such as including layers
and channel numbers. The standard UNet models use 3x3 convolution, so we take
kernel sizes of {1, 3, 5} as search space. We construct a supernet for each
block with each layer consisting of these three options.

Supernet training. For pixel-wise models, on CIFAR-10 and CelebA, we use the
checkpoints provided by the DDPM and DDIM repositories respectively as the
teacher. On LSUN-church, we train our own baseline model as the teacher. For
latent models, we use the checkpoints provided by the official repository
[26]. The supernet training for different blocks can be done in parallel to
speed up, because the input of each block comes from the teacher rather than
the previous block. We find that the $L2$ loss tends to be unchanged after
some epochs of training, even though further training continues to improve
model performance. Therefore, we only need to train the supernet for much
fewer epochs than training a baseline from scratch, and then we can do the
evaluation. For example, training a teacher model for LSUN-church needs nearly
450,000 steps, while training the supernet only needs 1,750 steps for each
block. When evaluating each subnet, we find computing for several batches is
enough which can save much time.

Subnet retraining. At the beginning of subnet retraining, we set $\beta$ to 0,
which means we only use the distillation loss. In this stage, the distillation
loss transfers knowledge from the teacher to the subnet. Then we increase
$\beta$ to assign a larger weight to the original noise prediction loss. Let
$\beta$-$steps$ denote the number of training steps when $\beta$ increases to
1. We test both the linear and the step functions. The linear function means
$\beta$ increases from 0 to 1 linearly. The step function means $\beta$
abruptly changes to 1 at $\beta$-$steps$. We test these two functions on
CIFAR-10 and CelebA. The experiments show that compared to the linear
strategy, the step strategy performs similarly on CIFAR-10 and noticeably
better on CelebA. Both of them surpass the teacher and the ablation models. we
recommend the step function as the preferred choice since it is more time-
efficient. We observe that the training is not sensitive to the value of
$\gamma$.

Figure 3: The curves of FID on CIFAR-10 with respect to MACs. The dashed lines
denote the FID of the teacher. Different colors represent different samplers.

###  5.2 Results

We use the most frequently-use metric FID [7] to measure the quality of image
generation. We also compute sFID for latent diffusion models to better measure
spatial similarity for high-resolution images. We test the classic sampler
DDIM [31] as well as advanced fast samplers, DPM-Solver [20] and PNDM [18]. We
use MACs (Multiply-Accumulate Operations) reduction to measure the compression
ratio of computational cost. More details about evaluation settings can be
found in the Appendix (Sec. 7.2).

The results on pixel-wise diffusion models. Table 1 shows the FID of our
lightweight models and the teacher on CIFAR-10. S1 is the searched model
corresponding to $r=1.02$, while S2 corresponds to $r=1.00$. “Random” is a
randomly selected architecture, while “No-Dis-S2” is a S2 model that is
retrained without the joint loss. “Fixed-Loss-S1” is a S1 model trained with
fixed rather than dynamic $\beta$. we will discuss these three models in the
ablation study. The results show that S2 with the linear function
significantly outperforms the teacher under all samplers, and S1 with the
linear function also outperforms the teacher except for ddim-100. We also find
that the S2 with the step function achieves a very similar FID to the S2 with
the linear function, except for dpm-20. Also, the S2 (step) noticeably
outperforms the teacher. As for compression ratio, S1 achieves a 31% MACs
reduction, while S2 achieves a 26% reduction. Table 2 shows the FID
measurement on CelebA. S1 corresponds to $r=1.00$. Both the linear and the
step function achieve better FID than the teacher. The step one performs
better than the linear one. S1 achieves a 25% MACs reduction. Table 3 shows
the FID measurement on LSUN-church 128x128. S1 corresponds to $r=1.03$. The S1
with the step function performs better than the teacher in general.

The results on latent diffusion models. Our method performs better on latent
diffusion models, which achieve about 50% MACs reduction and also Parameter
reduction. Table 4 shows the results on LSUN-church 256x256 and CelebA-HQ with
the latent diffusion models. Our compressed model has comparable performance
to the teacher with significant MACs reduction (46% and 45%) and Parameter
reduction (58% and 46%).

As mentioned earlier, we can obtain a series of search models with different
compression ratios by slightly adjusting $r$. Smaller subnets tend to have
larger evaluation losses, and vice versa. We test a series of models with
different MACs on CIFAR-10. Figure 3 shows the relations between FID and
subnet MACs on CIFAR-10. The smaller MACs corresponds to a searched model
under a larger $r$. The results indicate that the performance of subnets
increases with respect to MACs, which accords with the intuition that more
complex models usually have better performance.

Figure 4 shows the generated images of the original and our compressed latent
models on LSUN-church 256x256. Figure 5 shows the results on CelebA-HQ
256x256. These qualitative results show there is no visual difference between
our compressed model and the baseline. Samples of pixel-wise models can be
found in the Appendix.

Figure 4: The generated images samples on LSUN-church 256 of our lightweight
model and the teacher. The sampling strategy is DDIM with 50 steps. Figure 5:
The generated images samples on CelebA-HQ 256 of our lightweight model and the
teacher. The sampling strategy is DDIM with 50 steps.

###  5.3 Ablation Study

The effectiveness of NAS. Previous works [38, 40] suggest that a number of NAS
methods are no better than or struggle to outperform a random architecture. So
it is important to prove that our method can surpass the random architecture.
On CIFAR-10, We randomly select a model whose MACs nearly equals to S2 and
train it with the same setting as S2. Table 1 shows the comparison between the
random models and S2. The table shows that S2 (linear) outperforms Random
(linear) significantly except for a similar FID at dpm-20. S2 (step)
noticeably outperforms Random (step) under all samplers. These results prove
the effectiveness of our DiffNAS, i.e., DiffNAS really obtains a noticeably
better architecture than the random one.

The effectiveness of the dynamic joint loss. As mentioned before, we propose a
joint loss to maintain the consistency between supernet training and subnet
retraining, provide guidance, and shorten paths of gradient propagation. On
CIFAR-10, we train a model referred to as “No-Dis-S2“ that represents a S2
model retrained without the distillation loss. The results in Table 1 show
that both the S2 with the linear joint loss and the S2 with the step function
noticeably outperform “No-Dis-S2“. This proves the necessity of
$\mathcal{L}_{\text{dis}}$. Further, in order to prove the necessity to
dynamically adjust the weight $\beta$, we conduct another experiment where we
assign the same fixed weight to $\mathcal{L}_{\text{ori}}$ and
$\mathcal{L}_{\text{dis}}$, which is denoted by Fix-Loss-S1 in Table 1. The
results show the model has very poor performance, which proves the necessity
to reduce the weight of $\mathcal{L}_{\text{dis}}$ in the later steps of
training.

The effectiveness of block-wise search. As mentioned before, our method
proposes to block-wisely search for the most appropriate subnets to avoid the
unfairness among blocks of the global search strategy adopted by previous
works. This unfairness is due to the large differences of loss values among
different blocks, even after being balanced. For example, we observe that some
blocks have loss values about 0.50, while some other blocks have loss values
about 0.02. Therefore, the global search strategy which summarizes the losses
of all blocks will be dominated by the blocks that have large losses and
ignore the blocks with small losses. To prove this proposition, we conduct an
experiment that uses an evolutionary algorithm to globally search for the
smallest subnet, on the CIFAR-10 dataset. We use the relative $L2$ loss shown
below to balance the loss values among blocks:

$\displaystyle\mathcal{L}_{\text{val}}=\frac{||\text{Block}_{i}(\mathbf{X_{i}})-\mathbf{Y_{i}}||^{2}_{2}}{\text{var}(\mathbf{Y_{i}})}\
,$

where $\text{var}(\mathbf{Y_{i}})$ denotes the variance of the ground-truth.
After performing the evolutionary algorithm for 30 epochs, we obtain a subnet
that has a better evaluation loss than the teacher. However, this subnet has
an obviously poor performance compared to the teacher performance in Table 1,
with the FID of 5.08, 8.23, and 5.74 on ddim-50, dpm-20, and dpm-25
respectively. This proves that our block-wise search is much more effective
than the global one.

##  6 Conclusion

To reduce the computational cost of diffusion models, we propose DiffNAS which
automatically searches for the smallest architecture while maintaining the
performance as much as possible. DiffNAS contains a novel block-wise search
strategy and a novel retraining strategy with a dynamic joint loss. The block-
wise search strategy can avoid the unfairness of loss values among different
blocks. The dynamic joint loss can maintain the consistency between supernet
training and subnet retraining, and transfer teacher knowledge to improve
performance. Empirical results show that our method achieves about 30% MACs
reduction on pixel-wise models, as well as about 50% MACs and Parameters
reduction on latent diffusion models, while maintaining on-par or even
achieving better performance.

## References

  * Baker et al. [2016] Bowen Baker, Otkrist Gupta, Nikhil Naik, and Ramesh Raskar.  Designing neural network architectures using reinforcement learning.  In _International Conference on Learning Representations_ , 2016.
  * Chen et al. [2021] Nanxin Chen, Yu Zhang, Heiga Zen, Ron J Weiss, Mohammad Norouzi, Najim Dehak, and William Chan.  Wavegrad 2: Iterative refinement for text-to-speech synthesis.  _arXiv preprint arXiv:2106.09660_ , 2021.
  * Chu et al. [2021] Xiangxiang Chu, Bo Zhang, and Ruijun Xu.  Fairnas: Rethinking evaluation fairness of weight sharing neural architecture search.  In _Proceedings of the IEEE/CVF International Conference on computer vision_ , pages 12239–12248, 2021.
  * Dhariwal and Nichol [2021] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in Neural Information Processing Systems_ , 34:8780–8794, 2021.
  * Fang et al. [2023] Gongfan Fang, Xinyin Ma, and Xinchao Wang.  Structural pruning for diffusion models.  _arXiv preprint arXiv:2305.10924_ , 2023.
  * Guo et al. [2020] Zichao Guo, Xiangyu Zhang, Haoyuan Mu, Wen Heng, Zechun Liu, Yichen Wei, and Jian Sun.  Single path one-shot neural architecture search with uniform sampling.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XVI 16_ , pages 544–560. Springer, 2020.
  * Heusel et al. [2017] Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  _Advances in neural information processing systems_ , 30, 2017.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Ho et al. [2022a] Jonathan Ho, William Chan, Chitwan Saharia, Jay Whang, Ruiqi Gao, Alexey Gritsenko, Diederik P Kingma, Ben Poole, Mohammad Norouzi, David J Fleet, et al.  Imagen video: High definition video generation with diffusion models.  _arXiv preprint arXiv:2210.02303_ , 2022a.
  * Ho et al. [2022b] Jonathan Ho, Chitwan Saharia, William Chan, David J Fleet, Mohammad Norouzi, and Tim Salimans.  Cascaded diffusion models for high fidelity image generation.  _J. Mach. Learn. Res._ , 23(47):1–33, 2022b.
  * Karras et al. [2017] Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen.  Progressive growing of gans for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_ , 2017.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Li et al. [2020a] Changlin Li, Jiefeng Peng, Liuchun Yuan, Guangrun Wang, Xiaodan Liang, Liang Lin, and Xiaojun Chang.  Block-wisely supervised neural architecture search with knowledge distillation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1989–1998, 2020a.
  * Li et al. [2021] Changlin Li, Tao Tang, Guangrun Wang, Jiefeng Peng, Bing Wang, Xiaodan Liang, and Xiaojun Chang.  Bossnas: Exploring hybrid cnn-transformers with block-wisely self-supervised neural architecture search.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 12281–12291, 2021.
  * Li et al. [2020b] Xiang Li, Chen Lin, Chuming Li, Ming Sun, Wei Wu, Junjie Yan, and Wanli Ouyang.  Improving one-shot nas by suppressing the posterior fading.  In _Proceedings of the IEEE/CVF Conference on computer vision and pattern recognition_ , pages 13836–13845, 2020b.
  * Li et al. [2023] Xiuyu Li, Long Lian, Yijiang Liu, Huanrui Yang, Zhen Dong, Daniel Kang, Shanghang Zhang, and Kurt Keutzer.  Q-diffusion: Quantizing diffusion models.  _arXiv preprint arXiv:2302.04304_ , 2023.
  * Liu et al. [2018] Hanxiao Liu, Karen Simonyan, and Yiming Yang.  Darts: Differentiable architecture search.  _arXiv preprint arXiv:1806.09055_ , 2018.
  * Liu et al. [2021] Luping Liu, Yi Ren, Zhijie Lin, and Zhou Zhao.  Pseudo numerical methods for diffusion models on manifolds.  In _International Conference on Learning Representations_ , 2021.
  * Liu et al. [2015] Ziwei Liu, Ping Luo, Xiaogang Wang, and Xiaoou Tang.  Deep learning face attributes in the wild.  In _Proceedings of International Conference on Computer Vision (ICCV)_ , 2015.
  * Lu et al. [2022] Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan Li, and Jun Zhu.  Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps.  _Advances in Neural Information Processing Systems_ , 35:5775–5787, 2022.
  * Luhman and Luhman [2021] Eric Luhman and Troy Luhman.  Knowledge distillation in iterative generative models for improved sampling speed.  _arXiv preprint arXiv:2101.02388_ , 2021.
  * Lyu et al. [2022] Zhaoyang Lyu, Xudong Xu, Ceyuan Yang, Dahua Lin, and Bo Dai.  Accelerating diffusion models via early stop of the diffusion process.  _arXiv preprint arXiv:2205.12524_ , 2022.
  * Nash et al. [2021] Charliegenera Nash, Jacob Menick, Sander Dieleman, and Peter Battaglia.  Generating images with sparse representations.  In _International Conference on Machine Learning_ , pages 7958–7968. PMLR, 2021.
  * Nichol and Dhariwal [2021] Alexander Quinn Nichol and Prafulla Dhariwal.  Improved denoising diffusion probabilistic models.  In _International Conference on Machine Learning_ , pages 8162–8171. PMLR, 2021.
  * Ramesh et al. [2022] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.  Hierarchical text-conditional image generation with clip latents.  _arXiv preprint arXiv:2204.06125_ , 2022.
  * Rombach et al. [2022] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 10684–10695, 2022.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18_ , pages 234–241. Springer, 2015.
  * Salimans and Ho [2021] Tim Salimans and Jonathan Ho.  Progressive distillation for fast sampling of diffusion models.  In _International Conference on Learning Representations_ , 2021.
  * Shang et al. [2023] Yuzhang Shang, Zhihang Yuan, Bin Xie, Bingzhe Wu, and Yan Yan.  Post-training quantization on diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1972–1981, 2023.
  * Sohl-Dickstein et al. [2015] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International Conference on Machine Learning_ , pages 2256–2265. PMLR, 2015.
  * Song et al. [2020a] Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  _arXiv preprint arXiv:2010.02502_ , 2020a.
  * Song and Ermon [2020] Yang Song and Stefano Ermon.  Improved techniques for training score-based generative models.  _Advances in neural information processing systems_ , 33:12438–12448, 2020.
  * Song et al. [2020b] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2020b.
  * Vignac et al. [2022] Clement Vignac, Igor Krawczuk, Antoine Siraudin, Bohan Wang, Volkan Cevher, and Pascal Frossard.  Digress: Discrete denoising diffusion for graph generation.  _arXiv preprint arXiv:2209.14734_ , 2022.
  * Watson et al. [2021] Daniel Watson, Jonathan Ho, Mohammad Norouzi, and William Chan.  Learning to efficiently sample from diffusion probabilistic models.  _arXiv preprint arXiv:2106.03802_ , 2021.
  * Wu et al. [2019] Bichen Wu, Xiaoliang Dai, Peizhao Zhang, Yanghan Wang, Fei Sun, Yiming Wu, Yuandong Tian, Peter Vajda, Yangqing Jia, and Kurt Keutzer.  Fbnet: Hardware-aware efficient convnet design via differentiable neural architecture search.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 10734–10742, 2019.
  * Wu et al. [2022] Jay Zhangjie Wu, Yixiao Ge, Xintao Wang, Weixian Lei, Yuchao Gu, Wynne Hsu, Ying Shan, Xiaohu Qie, and Mike Zheng Shou.  Tune-a-video: One-shot tuning of image diffusion models for text-to-video generation.  _arXiv preprint arXiv:2212.11565_ , 2022.
  * Yang et al. [2019] Antoine Yang, Pedro M Esperança, and Fabio M Carlucci.  Nas evaluation is frustratingly hard.  In _International Conference on Learning Representations_ , 2019.
  * Yu et al. [2015] Fisher Yu, Ari Seff, Yinda Zhang, Shuran Song, Thomas Funkhouser, and Jianxiong Xiao.  Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop.  _arXiv preprint arXiv:1506.03365_ , 2015.
  * Yu et al. [2019] Kaicheng Yu, Christian Sciuto, Martin Jaggi, Claudiu Musat, and Mathieu Salzmann.  Evaluating the search phase of neural architecture search.  In _International Conference on Learning Representations_ , 2019.
  * Zhong et al. [2018] Zhao Zhong, Junjie Yan, Wei Wu, Jing Shao, and Cheng-Lin Liu.  Practical block-wise neural network architecture generation.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2423–2432, 2018.
  * Zoph and Le [2016] Barret Zoph and Quoc Le.  Neural architecture search with reinforcement learning.  In _International Conference on Learning Representations_ , 2016.

Supplementary Material

##  7 Extended Experimental Settings

###  7.1 Datasets

To demonstrate the generalization of our method, we conduct experiments on
several datasets with different image resolutions.

CIFAR-10 [12]. CIFAR-10 is a popular image classification dataset that
contains 60,000 32x32 RGB images from 10 different classes. CIFAR-10 is a
popular baseline for computer vision tasks such as image generation.

CelebA [19] and CelebA-HQ [11]. CelebA is a large-scale face attributes
dataset with 202,599 celebrity images. CelebA-HQ is a higher-resolution
dataset developed from CelebA. They are also commonly used in image generation
tasks.

LSUN-church [39]. LSUN-church is part of the Large-Scale Scene Understanding
(LSUN) datasets, consisting of 126,527 outdoor church images. LSUN-church
covers a wide range of church architectures and styles, making it diverse and
challenging. LSUN-church is a frequently-used dataset in image generation to
measure the performance of generative models.

###  7.2 Evaluation Settings

We use FID and sFID to measure the quality of generative models. FID is an
important metric used in image generation tasks to measure the similarity
between generated images and real images. A lower FID indicates higher
similarity and always means better quality of generative models. sFID [23] is
another frequently-used metric, a variant of FID, which can better measure
spatial similarity. For all pixel-wise models, we generate 50,000 samples and
compute the FID score between these generated images and the whole training
dataset. For latent diffusion models, we generate 10,000 samples and calculate
both FID and sFID. For the measurement of MACs and parameter amount, we use
the THOP package.

##  8 Computational Costs of Our Method

Our method only requires a little extra time cost and no extra memory cost
compared to training a baseline. Specifically:

  * •

The supernet training stage is fast. As reported in the paper, we only need a
few steps to train each block, e.g., 1,750 steps using pixel-wise diffusion
and 2,000 steps using latent diffusion, on LSUN-Church. Besides, because we
only need to train a block once a time, the memory cost is much less than
training a baseline. We cost about 6 hours on two A100 GPUs when testing on
Church 256x256 with the latent diffusion model.

  * •

The search stage is fast. As discussed in the paper, we evaluate each subnet
in each block with several batches, so the stage is fast. Besides, this stage
doesn’t need training, so it is also memory-efficient. We cost about 1.2 hours
in a single A100 GPU when testing on Church with the latent diffusion model

  * •

The retraining stage also does not require extra memory cost. Our method
retrains the compressed subnet, so it needs less GPU memory cost than training
a baseline. We cost about 20G memory (batch size=32 for each GPU) when testing
on Church with the latent diffusion model, while training a baseline costs
about 25G memory. The stage only needs a little extra time cost to do teacher
inference at the beginning when the beta is less than 1. Moreover, We also
find the number of converged epochs is less than the baseline.

##  9 Qualitative Results of Pixel-wise Diffusion Models

We show the generated images of pixel-wise diffusion models on CelebA 64x64
and LSUN-Church 128x128, i.e., Figure 6 and 7. These qualitative results show
there is no visual difference between our compressed model and the baseline.

Figure 6: The generated images on CelebA of our lightweight model and the
teacher. The sampling strategy is a DPM-Solver of 20 steps. Figure 7: The
generated images on LSUN-church 128. The sampling strategy is a DPM-Solver of
25 steps.


#  Emergent Time Crystal with Tunable Period
in the Uniaxial Random Field XY Magnet

S. Basak  Department of Physics and Astronomy, Purdue University, West
Lafayette, IN 47907, USA  Purdue Quantum Science and Engineering Institute,
West Lafayette, IN 47907, USA  K. A. Dahmen  Department of Physics, University
of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  E. W. Carlson
<EMAIL_ADDRESS>Department of Physics and Astronomy, Purdue University,
West Lafayette, IN 47907, USA  Purdue Quantum Science and Engineering
Institute, West Lafayette, IN 47907, USA

###### Abstract

The addition of uniaxial random fields to the XY model induces an order-by
disorder transition, in which the XY magnet develops a spontaneous
magnetization in the direction perpendicular to the uniaxial random field.
Here, we use simulations to explore the robustness of this phase transition
with respect to a rotating driving field. We find that the order-by-disorder
transition is robust, persisting to finite applied field. In the vicinity of
the critical driving field strength, a time crystal emerges, in which the
period of the limit cycles becomes an integer $n>1$ multiple of the driving
period. Because $n$ increases with system size, the period of the time crystal
can be engineered. This period multiplication cascade is reminiscent of that
occuring in amorphous solids subject to oscillatory shear near the onset of
plastic deformation, and of the period bifurcation cascade near the onset of
chaos in nonlinear systems, suggesting it is part of a larger class of
phenomena in transitions of dynamical systems. Applications include magnets,
electron nematics, and quantum gases.

##  I Introduction

The XY model, in which interacting spins are confined to rotate within a
plane, has been a staple of statistical mechanics and condensed matter
studies, having been applied to a broad range of physical systems including
planar magnets, superfluids, superconductors, two-dimensional melting, nematic
liquid crystals, and electron nematics, among others. Lacour-Gayet, P. and
Toulouse, G. (1974); Aoyama and Ikeda (2006); Wehr  _et al._ (2006); Ostlund
and Halperin (1981); Nelson and Halperin (1979); Halperin and Nelson (1978);
Young (1979); Lilly  _et al._ (1999); Abanin  _et al._ (2007); Girvin and
MacDonald  In two dimensions, the XY model exhibits a BKT transition to a
power-law ordered phase, yet with no long-range order.Berezinskii (1971);
Kosterlitz and Thouless (1973) As such, the addition of random fields to a
two-dimensional XY model is expected to result in even less order: Imry and Ma
argued that a ($d\leq 4$)-dimensional system with continuous order parameter
(with O(n) symmetry with $n\geq 2$) in the presence of random fields cannot
have long range order for any finite disorder strength. Imry and Ma (1975)

However, the addition of uniaxial random fields reduces the global symmetry of
the Hamiltonian, and the Imry-Ma argument no longer applies.Crawford (2013) In
this case, the low temperature phase has long-range order via an order-by-
disorder transition, in which XY spins align perpendicular to the random
fields.Crawford (2013); Minchau and Pelcovits (1985) This is a special case of
a more general class of order-by-disorder transition, where an n-dimensional
spin system orders in a (n-k)-dimensional subspace due to orthogonal
k-dimensional random fields. Crawford (2013); Fel’dman (1999); Wehr  _et al._
(2006); Abanin  _et al._ (2007)

In this paper, we consider the possibility of a nonequilibrium transition. We
use simulations to study the order-by-disorder transition in the presence of a
rotating driving field. By analyzing the avalanche size distribution as a
function of magnitude of applied driving field. We find evidence that the
system undergoes a continuous nonequilibrium phase transition at a critical
amplitude of the driving field. Once a limit cycle is established, we observe
that the period of the hysteresis loops become n-fold near a critical applied
field strength, where n is as large as 7 in our largest systems. We present
evidence that the period of the subharmonic entrainment is rigid against
perturbations in initial conditions, and against perturbations of the drive
field, indicating that a classical discrete time crystal emerges at
criticality.Yao  _et al._ (2018); Shapere and Wilczek (2012) We present finite
size scaling evidence that the period of these multi-period limit cycles will
diverge in the thermodynamic limit. An experimental test of this would be the
presence of non-repeatability in the response due to a rotating driving field
near the transition.

As discussed further in Sec. IV.1, there are several experimental systems
corresponding to the XY model into which uniaxial random field disorder can be
incorporated, whereby these ideas can be tested experimentally. These include
layers of Josephson junctions,Lacour-Gayet, P. and Toulouse, G. (1974)
superfluid in a uniaxially stressed aerogel,Aoyama and Ikeda (2006) ultracold
atoms in the presence of speckle radiation,Wehr  _et al._ (2006) uniaxially
stressed 2D Wigner crystals,Ostlund and Halperin (1981); Nelson and Halperin
(1979); Halperin and Nelson (1978); Young (1979) the half-integer quantum Hall
effect,Lilly  _et al._ (1999) and possibly the graphene quantum Hall
ferromagnet. Abanin  _et al._ (2007); Girvin and MacDonald

##  II Model

We consider the uniaxial random field XY model on a square lattice, in the
presence of a driving applied field $H$:

$\displaystyle\mathcal{H}=-J\sum_{\langle
i,j\rangle}\cos{(\theta_{i}-\theta_{j})}$ $\displaystyle-$
$\displaystyle\sum_{i}h_{i}\cos{(\theta_{i})}$ (1) $\displaystyle-$
$\displaystyle H\sum_{i}\cos{(\theta_{i}-\phi)},$

where $\vec{S}_{i}\equiv(\cos(\theta_{i}),\sin(\theta_{i}))$ is the XY spin on
each site $i$, and $J$ is the nearest neighbor interaction strength. The
second term arises from the interaction of a local random field along the
x-axis and the XY spins. We choose the random field $h_{i}$ at each site $i$
from a gaussian probability distribution of width $R_{x}$,
$P(h_{i})=exp[-h_{i}^{2}/(2R_{x}^{2})]/(\sqrt{2\pi R_{x}^{2}})$. The order
parameter is the magnetization per site
$\vec{m}=\frac{1}{N}\sum_{i=1}^{N}\vec{S}_{i}$, where $N=L\times L$ is the
number of sites.

We study this system under the influence of a rotating applied driving field
whose angle $\phi=\Omega t$ advances in time slowly compared to all other
timescales in the problem. The dynamics is quasi-static: after each small
increment of the driving field angle, the energy of the system is minimized.
(See Sec. V for details of the simulation method.) This type of dynamicsRegev
_et al._ (2015) presupposes that the system is connected to a heat bath which
prevents heating by the drive.

$\pi/2$$\pi$$3\pi/2$$2\pi$$5\pi/2$-1-0.500.51$\pi/2$$\pi$$3\pi/2$$2\pi$$5\pi/2$$\pi/2$$\pi$$3\pi/2$$2\pi$$5\pi/2$-1-0.500.51-1-0.500.51$\pi/2$$\pi$$3\pi/2$$2\pi$$5\pi/2$-1-0.500.51(a)

$\phi(\vec{H})$

$\phi(\vec{H})$$H=0.150J$$H=0.070J$$H=0.041J$$H=0.020J$(b)

$m_{y}(\phi)$

$m_{x}(\phi)$(c)(d)

(a)  Figure 1:  Steady state response to rotating applied field. A system of
size $N=160\times 160$ with $R_{x}=0.5J$ is started from an initial applied
field in the $y$-direction. The initial spin configuration is aligned with the
applied field, then relaxed acccording to Eqns. 3 as described in the text,
after which the applied field is rotated counterclockwise as denoted in panel
(b). Panels (a), (c), and (d) show the response once steady state is reached
under the driving field. Panel (a) shows the response of the magnetization in
the $x$ direction, while panel (d) shows the response of the magnetization in
the $y$ direction. Panel (c) is a parametric plot of $m_{y}$ vs $m_{x}$. In
all panels, the arrows denote the state of the system when the driving field
is at an angle $\phi=\pi$, i.e. aligned along the $x$ direction. For driving
field strength $H=0.041J$, the response of the system has double the period of
the driving field. The open arrow on this trace denotes the state of the
system at driving field angle $\phi=\pi$ during every other cycle of the
driving field.

##  III Results

###  III.1 Behavior of the Limit Cycles

Fig. 1 shows the rich behavior of the limit cycles in rotating driving field,
as a function of the magnitude of the driving field $H$ at intermediate
disorder strength $R_{x}=0.5J$. Panel (b) shows the sense of the driving
field, which is held at constant magnitude, but rotated counterclockwise, i.e.
$\phi$ increases in time as $\phi=\omega t$ in the $\omega\rightarrow 0$
limit, starting from $\phi=\pi/2$. Fig. 1 (a) shows a plot of $m_{x}$ vs. the
angle $\phi$ of the applied field. Panel (d) shows a plot of $m_{y}$ vs. the
angle $\phi$ of the applied field. Panel (c) shows the combined parametric
plot of magnetization $m_{x}$ in the $x$ direction, plotted against the
magnetization $m_{y}$ in the $y$ direction. The sense of the parametric plot
in panel (c) is counterclockwise. In each case, the system is started from a
locally stable configuration in applied field $\vec{H}||y$ at zero
temperature, which has been relaxed from an initially saturated state aligned
with the initial applied field. The transient response before the limit cycle
is not shown in this figure. We discuss the transient response in Sec. III.3.

For moderate disorder strength $R_{x}=0.5J$, we find that at small amplitudes
of the driving field, the spontaneous magnetization in the $y$ direction
remains robust. This is evident in the small hysteresis loops we find for
$H=0.02J$ as shown by the purple trace in the parametric plot Fig. 1(c). This
indicates that the system continues to display spontaneous symmetry breaking
in the $y$ direction, retaining its Ising ferromagnetic character in the
presence of weak rotating driving field.

As the magnitude of the applied field is increased, there is a change in
behavior from ferromagnetic to paramagnetic response. This is evident in the
large, almost circular hysteresis loop we find for larger $H=0.15J$, as shown
by the red trace in the parametric plot Fig. 1(c). This change is consistent
with either a crossover in behavior or a non-equilibrium phase transition at a
critical magnitude of the driving field. Note that the rotating hysteresis
loops at intermediate driving field strengths $H=0.041J$ and $H=0.07J$ have
rich structure: Numerous avalanches are evident in these traces. As we will
see in Sec. III.2, the avalanche structure provides further insight into the
question of whether the change from ferromagnetic to paramagnetic response is
a crossover or a phase transition. Perhaps the most intriguing feature of the
intermediate driving field regime is that in the blue trace ($H=0.041J$), the
limit cycle has double the period of the driving field. We find that limit
cycles often become multiperiodic at intermediate field strength, for large
enough system size. We explore this region of the phase diagram in more depth
in Sec. III.4.

$0$$5000$$10000$$15000$$20000$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$

Largest Avalanche ($\overline{|\delta\vec{M}|_{max}}$)

Rotating field strength
($H/J$)$N=64^{2}$$N=80^{2}$$N=100^{2}$$N=128^{2}$$N=160^{2}$

(a)  $T\rightarrow 0,\ R_{x}=0.5J$

$0$$1\times 10^{8}$$2\times 10^{8}$$3\times 10^{8}$$4\times 10^{8}$$5\times
10^{8}$$6\times 10^{8}$$7\times 10^{8}$$8\times
10^{8}$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$

Second moment of Avalanches ($\overline{\langle(\delta\vec{M})^{2}\rangle}$)

Rotating field strength
($H/J$)$N=64^{2}$$N=80^{2}$$N=100^{2}$$N=128^{2}$$N=160^{2}$

(b)  $T\rightarrow 0,\ R_{x}=0.5J$

Figure 2:  Avalanche statistics. The response of the magnetization to rotating
driving field often proceeds via avalanches, in which there is a discontinous
jump in the magnetization $\delta M$ in response to a small change
$\delta\phi$ of the driving field angle. In panel (a), we plot the size of the
largest avalanche $|\delta\vec{M}|_{\rm max}$ per limit cycle at each rotating
field strength, disorder-averaged, for a range of system sizes. Panel (b)
shows the disorder-average of the second moment $\delta\vec{M}$ of the
avalanche size distribution. The brackets $\left<~{}~{}\right>$ denote an
average over the limit cycle, and the overbar denotes a disorder average. By
both of these measures, the size of the avalanches grows with system size
implying divergent fluctuations at a critical field strength in the
thermodynamic limit. The vertical bars in both panels mark the peak value from
a running 3-point average. Within the resolution of the plot in panel (a),
these values are concident for sizes $N=80\times 80$ and $N=100\times 100$,
and for sizes $N=128\times 128$ and $N=160\times 160$. In panel (b), the peak
values are coincident for sizes $N=64\times 64$ and $N=160\times 160$.

$-1$$-0.5$$0$$0.5$$1$$-1$$-0.5$$0$$0.5$$1$

$m_{y}(\phi_{\vec{H}})$

$m_{x}(\phi_{\vec{H}})$${\rm Transient}$${\rm Limit\ Cycle}$

(a)  $R_{x}=0.5J;~{}H=0.058J;~{}N=64\times 64.$

$-1$$-0.5$$0$$0.5$$1$$-1$$-0.5$$0$$0.5$$1$

$m_{y}(\phi_{\vec{H}})$

$m_{x}(\phi_{\vec{H}})$${\rm Transient}$${\rm Limit\ Cycle}$

(b)  $R_{x}=0.5J;~{}H=0.048J;~{}N=100\times 100.$

$-1$$-0.5$$0$$0.5$$1$$-1$$-0.5$$0$$0.5$$1$

$m_{y}(\phi_{\vec{H}})$

$m_{x}(\phi_{\vec{H}})$${\rm Transient}$${\rm Limit\ Cycle}$

(c)  $R_{x}=0.5J;~{}H=0.046J;~{}N=160\times 160.$

$0$$\pi/4$$\pi/2$$3\pi/4$$\pi$$5\pi/4$$0.03$$0.04$$0.05$$0.06$$0.07$

Duration of Transient ($\overline{\Delta\phi}$)

Rotating field strength
($H/J$)$N=64^{2}$$N=80^{2}$$N=100^{2}$$N=128^{2}$$N=160^{2}$

(d)  $\ R_{x}=0.5J,\ R_{y}=0$

Figure 3: Transient response and multiperiod limit cycles near the transition.
Panels (a-c) show the initial transient response (black curves), followed by
multiperiodic limit cycles (rainbow curves). (a) Transient response and
multiperiodic limit cycle for one disorder configuration at $N=64^{2}$. Here,
the transient response lasts roughly half a cycle before a period-2 limit
cycle is established. (b) Transient response and multiperiodic limit cycle for
one disorder configuration at a larger system size $N=100^{2}$. Here, the
transient response lasts roughly one cycle before a period-2 limit cycle
appears. (c) Transient response and multiperiodic limit cycle for one disorder
configuration at an even larger system size $N=160^{2}$. Here, the transient
response lasts almost 1.5 cycles before a period-3 limit cycle is established.
(d) The disorder-averaged duration of the transient response, as a function of
$H$. The mean of the transient distribution function for each system size is
marked by a vertical line of the corresponding color.

###  III.2 Avalanches Near the Transition

In this section we focus on the characteristics of the avalanches that occur
near the transition from Ising ferromagnetic to paramagnetic response. We find
a rich avalanche structure at intermediate field strengths, as can be seen in
the blue and green traces in Fig. 1 ($H=0.041J$ and $H=0.07J$, respectively).
Notice that while the avalanches are apparent in both $m_{x}$ and in $m_{y}$,
they are most prominent in $m_{y}$, which serves as the order parameter in
this system. When magnetization is cast as an extensive quantity,
$\vec{M}=\sum_{i=1}^{N}\vec{S}_{i}=N\vec{m}$, then in the thermodynamic limit,
avalanches $\delta\vec{M}$ of diverging size accompany a second order phase
transition.

Fig. 2a plots the size of the largest avalanche $|\delta\vec{M}|_{\rm max}$ at
each rotating field strength, for a range of system sizes $N=L\times L$.
Results are averaged over several disorder configurations of the random field
at disorder strength $R_{x}=0.5J$, ranging from $75$ disorder configurations
for system size $N=64^{2}$, to $30$ disorder configurations for system size
$N=160^{2}$. (See Appendix V.2.) Notice that fluctuations as measured by the
largest avalanche diverge with increasing system size at a critical driving
field strength, $H_{c}(R_{x}=0.5J)$. We estimate the value of $H_{c}$ at
$R_{x}=0.5J$ as follows: For each system size, the peak value based on a
3-point average is indicated by the vertical bar. The corresponding peak value
of the applied field strength, averaged over all system sizes, is
$H_{c}=(0.0452\pm 0.0015)J$.

In Fig. 2b, we plot the second moment of all avalanches in each limit cycle,
$\left<(\delta\vec{M})^{2}\right>$ at each rotating driving field strength,
for a range of system sizes. Results are disorder averaged, using the same
number of disorder configurations as in Fig. 2a. Notice that this alternate
measure of fluctuations based on the second moment of the avalanche size
distribution is also consistent with the system undergoing a second order,
nonequilibrium phase transition at a critical driving field strength, $H_{c}$.
In this case, we find that $H_{c}(R_{x}=0.5J)=(0.0432\pm 0.0016)J$, in
agreement with the value of the critical field strength we find from Fig. 2a.

Figure 4:  Example of spin configurations during a period-2 limit cycle. Spin
configurations (a) and (c) are for the same angle $\phi$ of the driving field,
but the spin configuration is different the second time through the driving
cycle. Likewise, spin configurations (b) and (d) are for the same angle $\phi$
of the driving field, but the spin configuration is different the second time
through the driving cycle. For this particular disorder configuration and
system size, the spin configurations repeat every 2 periods of the driving
cycle. Here, the driving field strength is $H=0.04J$, and the system size is
160$\times$160\. See Supplementary Information for further examples, and links
to videos.

00.050.10.150.20.250.300.050.10.1500.050.100.0500.050.040.050.06

\makebox(0.0,0.0)[]{{} }

$7\times 2\pi$$5\times 2\pi$$4\times 2\pi$$3\times 2\pi$$2\times 2\pi$(a)
$N=160^{2}$

# Multi-period loops / # All loops

$\langle H_{lc}/J\rangle_{N}$(b) $N=128^{2}$

\makebox(0.0,0.0)[]{{} }

(c) $N=100^{2}$

\makebox(0.0,0.0)[]{{} }

(d) $N=80^{2}$

\makebox(0.0,0.0)[]{{} }

$H/J$(e) $N=64^{2}$

(a)

$10^{-4}$$10^{-3}$$10^{-2}$$160^{-2}$$128^{-2}$$100^{-2}$$80^{-2}$$64^{-2}$

$\langle H_{lc}/J\rangle_{N}-0.0422$

$N^{-1}$$1.2\times 10^{14}\times N^{-4.13}+0.0442$$1.2\times 10^{8}\times
N^{-2.63}+0.0439$$7.2\times 10^{2}\times N^{-1.27}+0.0422$$\langle
H_{lc}/J\rangle_{N}$(f)

(b)

Figure 5: Multiperiod limit cycles. Panels (a-e) show what fraction of limit
cycles exhibit multiperiodicity as a function of driving field strength $H$ at
$R_{x}=0.5J$. The smallest system size we simulated, $N=64\times 64$, is shown
on the bottom left in panel (e). System size increases from bottom to top in
the left panels, up to system size $N=160\times 160$. In the bar graphs,
period-2 limit cycles are shaded pink; period-3 limit cycles are purple;
period-4 limit cycles are green; the period-5 limit cycle is orange; and the
period-7 limit cycle is red. We did not observe any period-6 limit cycles. In
each bar graph, the vertical blue line is the mean of the distribution
function, $\left<H_{lc}\right>$ in units of $J$. (f) From the results of
panels (a-e), we plot $\langle H_{lc}\rangle_{N}$ vs. the inverse of system
size $N$ on a log-log scale. A power law fit of $\langle H_{lc}\rangle_{N}$
for the three largest system sizes is given by the red curve; the fit for the
four largest system sizes is given by the green curve; and the fit for all
calculated system sizes is given by the dark blue curve. The $y$-intercept is
consistent among all of these fits, yielding an average value of $\langle
H_{lc}\rangle_{N\to\infty}=(0.0434\pm 0.0020)J$ .

###  III.3 Transient Response

Fig. 3(a-c) shows how the magnetization responds to a rotating driving field
in the vicinity of the phase transition. There is a transient response before
the system settles into a limit cycle. A limit cycle is the steadily repeating
response in the magnetization due to a rotating driving field. While we find
that most limit cycles have the same period as the driving field, we find that
near the transition regime, limit cycles often have a longer period. We first
discuss the behavior of the transient response, before turning our attention
to the behavior of the multiperiodic limit cycles in Sec. III.4

The transient response in panels (a-c) of Fig. 3 is marked in black. In Fig.
3(d), we plot the duration of the transient response, as a function of $H$,
for various system sizes. The results shown have been averaged over several
disorder configurations. (See Appendix V.2 for details.) At high and low
strength of the driving field, the transient response becomes so negligible as
to be smaller than the symbol size on this graph. However, at intermediate
driving field strength, the transient response grows with increasing system
size. The fact that the transient response grows with increasing system size
is further corroboration that the system is undergoing a second order phase
transition. In Fig. 3, the mean of each transient distribution function is
denoted by a vertical line, color coded to the system size. The average of the
mean value of $H$ from these vertical lines is $\langle
H_{tr}\rangle_{N}=(0.0430\pm 0.0014)J$, consistent with our previous estimates
of $H_{c}(R_{x}=0.5J)$.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0$$160^{-2}$$128^{-2}$$100^{-2}$$80^{-2}$$64^{-2}$

Maximum periodicity ($n_{max}$)

Inverse System Size ($N^{-1}$)$2.3\times 10^{-4}\times N^{1.22}+1$

(a)  Maximum period of limit cycles.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$160^{-1}$$128^{-1}$$100^{-1}$$80^{-1}$$64^{-1}$

Maximum likelihood ($P_{max}[n>1]$)

Inverse System Length ($L^{-1}$)$3.0\times\exp{(-367.10\times L^{-1})}$

(b)  Maximum likelihood of multiperiod limit cycles.

Figure 6:  Trends of the multiperiodic behavior of the limit cycles with
increasing system size. In panel (a) we plot the maximum period of the limit
cycles observed in Fig 5(a-e), as a function of 1/N (purple circles). The
black line is a fit to the simulation results. The trend is toward divergence
of the period of limit cycles in the thermodynamic limit. In panel (b), we
plot the maximum likelihood of multiperiodic limit cycles, obtained from the
peak heights of the left-hand panels in Fig 5 (purple circles). The black line
is a fit to the simulation results. The trend is toward saturation of the
likelihood of multiperiod behavior in the thermodynamic limit.

###  III.4 Period Increase Near the Transition

We now turn our attention to the behavior of the limit cycles at intermediate
driving field strength. One of the most fascinating features of the limit
cycles in this regime is that some of them have a longer period than that of
the driving field. Fig. 3 shows some representative cases of this behavior.
Fig. 4 visualizes how the spin configurations respond to the driving field
during one of the period-2 limit cycles. Domain walls have dramatically
different configurations during the second cycle as opposed to the first cycle
of the driving field, suggesting a prominent role for domain wall pinning and
domain wall creep. More examples of such behavior, including links to videos
of spin configurations during multiperiod limit cycles, can be found in the
Supplementary Information.

In order to explore this behavior quantitatively, we studied several disorder
configurations near the transition, as a function of system size. Fig. 5 shows
a histogram of the likelihood of multiperiod limit cycles. For a given
magnitude of the driving field $H$ and a given system size $N$, we plot the
number of disorder configurations whose limit cycle has a period greater than
that of the driving field, divided by the number of all disorder
configurations studied at that $H$ and $N$. Starting from the bottom panel on
the lefthand side of Fig. 5, panel (e), the system size increases as one moves
to the next panel up the page, up to panel (a) which shows the largest system
we studied, $N=160\times 160$. Different color bars indicate the period of the
multiperiod behavior: pink indicates period doubling; blue shows period
tripling; period-4 limit cycles are denoted in green; yellow is for period-5,
and orange is for period-7. We did not observe any period-6 limit cycles,
although presumably these would appear at certain disorder configurations as
well.

The vertical blue bars mark the mean of the distributions in Fig. 5(a-e),
$\langle H_{lc}\rangle_{N}$. In Fig. 5(f), we plot $\langle H_{lc}\rangle_{N}$
vs $1/N$ on a log-log plot, in order to determine the limiting value $\langle
H_{lc}\rangle_{N\to\infty}$. Fits of the finite size scaling in Fig. 5(f) for
all system sizes, the four largest system sizes, and the three largest system
sizes yield a consistent value for $\langle H_{lc}\rangle_{N\to\infty}$ within
error bars. The average of these three methods yields $\langle
H_{lc}\rangle_{N\to\infty}=(0.0434\pm 0.0020)J$.

$0$$1\times 10^{8}$$2\times 10^{8}$$3\times 10^{8}$$4\times 10^{8}$$5\times
10^{8}$$6\times 10^{8}$$7\times 10^{8}$$8\times
10^{8}$$0.02$$0.04$$0.06$$0.08$$0.1$$0.12$$0.14$

Second moment of Avalanches ($\overline{\langle(\delta\vec{M})^{2}\rangle}$)

Rotating field strength
($H/J$)$H_{c}$$N=64^{2}$$N=80^{2}$$N=100^{2}$$N=128^{2}$$N=160^{2}$${\rm
Ising}$

${\rm Non-Repeatability}$

${\rm Paramagnetic}$${\rm FM}$${\rm XY}$

(a)  Figure 7:  This figure shows the phase diagram as a function of the
strength of the rotating field in a uniaxial Random field. The region where
the number of multiperiodic loops and the maximum periodicity increases with
system size is labeled as the expected region of non-repeatablity for
infinitely large systems. This region coincides with the region where the
largest avalanche occurs is this system where $\langle H_{c}\rangle=(0.0437\pm
0.0009)J$, which is marked by the vertical blue line.

###  III.5 Approach to Non-Repeatability

We find that at small system size, multiperiod behavior is rare. However, as
the system size is increased, and the disorder configurations can become
correspondingly more rich, the likelihood of multiperiod behavior increases.
In Fig. 6(a), we plot the maximum observed period of a limit cycle, vs. $1/N$.
The maximum period increases with increasing system size, in a manner
consistent with diverging period in the thermodynamic limit.

Notice also that the distribution in Fig. 5(a-e) grows in height with
increasing system size. For $N=160\times 160$, we find that $20-30\%$ of
disorder configurations in the range $H=(0.04-0.046)J$ display multiperiodic
behavior. To quantify these effects, we plot the maximum height of the
distributions in Fig. 5(a-e) in Fig. 6(b). This measure also shows sharp
increase with increasing system size. The fact that both the likelihood of
multiperiod behavior and the period of limit cycles steadily increase with
increasing system size points toward a thermodynamic limit in which the period
of limit cycles goes to infinity. If the period of a system diverges in the
thermodynamic limit, then the system has effectively entered a regime of non-
repeatability. We discuss further implications of this finding in the next
section.

##  IV Discussion

Using four different methods to quantify the fluctuations in the system (see
Table 1), we find evidence for a second order nonequilibrium phase transition
from spontaneous Ising ferromagnetism at low driving field strength, to XY
paramagnetism at high driving field strength. The critical field strength at
which this transition occurs is consistent across all methods we employed,
yielding an average value of $H_{c}=0.0437\pm 0.0009$, as denoted in the phase
diagram in Fig. 7.

Table 1: Critical field strength.  Method | Value of $H_{c}/J$
---|---
Largest avalanche of limit cycle | $0.0452\pm 0.0015$
Second moment of avalanches in limit cycle | $0.0432\pm 0.0016$
Duration of transient response | $0.043\pm 0.0014$
Finite size scaling of multiperiodic behavior | $0.0434\pm 0.0020$
Overall average of above methods | $0.0437\pm 0.0009$

We furthermore find that far from being irrelevant, disorder plays a prominent
role at the transition. Because the disordered energy landscape makes the
system highly susceptible to spatial fluctuations near the transition, there
is both a longer transient response and a longer period of limit cycles near
$H_{c}$. Remarkably, both the likelihood of multiperiod behavior and the
period of the limit cycles increases with no sign of saturation as system size
is increased. The trend we find is toward a thermodynamic limit in which limit
cycles never repeat. A large enough physical system at this critical point
should therefore display a regime of non-repeatability. As shown in Fig. 7,
the regime of non-repeatability in the thermodynamic limit coincides with the
nonequilibrium phase transition. The dependence of this simple model upon
history implies that experiments on XY systems in uniaxial random field are
particularly sensitive to disorder. Conflicting experimental results could
arise if hysteresis protocols are not closely monitored.

Similar behavior is predicted to occur in models of amorphous solids under
periodic shear stress.Regev  _et al._ (2013, 2015); Leishangthem  _et al._
(2017). In these systems, simulations revealed that under periodic shear, the
response of the system becomes multiperiodic, in a way that is consistent with
chaotic behavior at a critical shear amplitude. More work would be needed to
determine whether the multiperiodic cascade observed here is indicative of
chaotic behavior in the thermodynamic limit. Similar multiperiod cascades
signal the onset of chaos in nonlinear systems, suggesting that the
multiperiod cascades observed here and in periodically driven models of
amorphous solids are characteristic of a larger class of transitions in
dynamical systems.

On the other hand, we predict that finite size physical systems in the
vicinity of the nonequilibrium transition should display the characteristics
of a classical discrete time crystal,Yao  _et al._ (2018); Shapere and Wilczek
(2012) in which the discrete time translation symmetry imposed by the periodic
drive is broken in a way that leads to rigid subharmonic entrainment. We find
that the period of the response remains stable against perturbations in the
initial conditions and stable against low temperature fluctuations (see the
SI), indicating that the spontaneous breaking of the discrete time symmetry is
rigid. Yao et al. find that the critical endpoint between a classical discrete
time crystal and the disordered phase of a dissipative, coupled chain of
classical nonlinear pendula terminates in a critical point which is not in the
Ising universality class. Because the nonequilibrium transition we find here
is in the Ising universality class, this indicates that there is more than one
classical discrete time crystal universality class. The results here further
underscore the fact that long-range interactions are not a necessary
ingredient to stabilize a time crystal.Yu  _et al._ (2019)

The work in this paper was done at uniaxial random field strength
$R_{x}=0.5J$, with zero random field strength in the $y$-direction. Further
work is needed to obtain the full phase diagram as a function of random field
strengths $R_{x}$ and $R_{y}$.

###  IV.1 Applications to Physical Systems

The uniaxial random field XY model has been applied to many systems, including
layers of Josephson junctions,Lacour-Gayet, P. and Toulouse, G. (1974)
superfluid in a uniaxially stressed aerogel,Aoyama and Ikeda (2006) ultracold
atoms in the presence of speckle radiation,Wehr  _et al._ (2006) uniaxially
stressed 2D Wigner crystals,Ostlund and Halperin (1981); Nelson and Halperin
(1979); Halperin and Nelson (1978); Young (1979) and the half-integer quantum
Hall effecteLilly  _et al._ (1999) Uniaxial random field-induced has also been
discussed in connection with the graphene quantum Hall ferromagnet. Abanin
_et al._ (2007); Girvin and MacDonald  We discuss below a few of these systems
in which there is also a clear way to drive the system with a rotating field.

####  IV.1.1 Electron nematics

An electron nematic occurs when the electronic degrees of freedom
spontaneously break the rotational symmetry of the host crystal. Electron
nematics have been observed or proposed in several material systems, including
transition metal oxides like cuprate superconductors, manganites, nickelates,
and cobaltites; valley symmetry breaking systems like single and bilayer
graphene, elemental bismuth, and AlGaAs 2DEG’s, as well as strontium
ruthenates and iron pnictides.Kivelson  _et al._ (2003); Fradkin and Kivelson
(2010) In mapping any nematic to an XY model, there is a factor of two between
the physical angle of the nematic in the plane, and the natural angles in an
XY model. This is because a nematic is symmetric under $180^{o}$ rotation,
whereas the XY spins change sign under the same operation. For an electron
nematic arising out of a crystal with discrete C4 rotational symmery, the
nematic order parameter switches sign when the nematic rotates by $90^{o}$.
The uniaxial random fields we discuss in this paper can arise in these systems
if random orienting fields are strong only along the major crystalline axes.
Note that in this case, the order-by-disorder transition would induce the
electron nematic to orient along a direction which is diagonal to the major
crystalline axes.

Several external perturbations can be used as a driving field on an electron
nematic, including magnetic field, electric field, high currents, and uniaxial
stress.Carlson and Dahmen (2011) Note that similar symmetry considerations
apply to the driving field in these systems. For example, a rotating applied
magnetic field $\vec{B}=[B_{x},B_{y}]=B[{\rm cos}(\omega t)),{\rm sin}(\omega
t)]$ can be used to exert the rotating driving field of Eqn. 1 for the case of
a nematic, with the caveat that rotating the applied field by $90^{o}$ changes
the sign of the driving field:

$\vec{H}=[H_{x},H_{y}]=H[{\rm cos}(2\omega t),{\rm sin}(2\omega t)]~{}.$ (2)

####  IV.1.2 Quantum Gases

Random-field induced order has been proposed to happen in coupled Bose-
Einstein condensade systems.Wehr  _et al._ (2006) Theoretical and numerical
results on two-component Bose gases predict that by using a Raman field to
couple two internal states, uniaxial random field disorder can be produced.
The uniaxial nature is achieved by a Raman coupling with constant phase, while
the randomness is achieved through random strength of the Raman
field.Niederberger  _et al._ (2008); Sanchez-Palencia and Lewenstein (2010)
Similarly, a rotating driving field can be applied by a Raman coupling with
uniform strength, but rotating phase.

####  IV.1.3 Magnetic systems

While the mapping of a magnetic system with XY symmetry to Eqn. 1 is clear,
the realization of a uniaxial random field in these systems is less clear. It
may be possible to design a system in which epitaxial strain from a substrate
exerts random uniaxial fields on a 2D XY ferromagnet through magnetoelastic
coupling.

###  IV.2 Conclusions

We have shown that the order-by-disorder transition of the two-dimensional XY
model in the presence of a uniaxial random field persists up to a critical
strength of the rotating driving field. Near the critical driving field
strength, the response of the system has a period which is an integer multiple
$n>1$ of the driving field period. This spontaneous breaking of the discrete
time symmetry of the driving field indicates that a classical discrete time
crystal emerges at the critical point. The trend with increasing system size
is toward increasing period $n$, indicating both that the period of the time
crystal can be engineered in small systems, and also suggesting the onset of
what is effectively non-repeatability as $n\rightarrow large$ in the
thermodynamic limit. Similar multiperiod cascades signal the onset of chaos in
nonlinear systems, and signal the onset of irreversibility in periodically
driven models of plastic deformation, suggesting that multiperiod cascades are
characteristic of a larger class of transitions in dynamical systems.

##  V Methods

The magnetization $m_{y}$ in the $y$-direction at intermediate disorder
strength $R_{x}=0.5J$ remains ordered even in the presence of weak applied
transverse field $H_{x}$ (see Supplementary Information). Therefore, to begin
the hysteresis studies, we first initialize the system in a $y$-magnetized
state, by starting from the fully saturated $y$ magnetization, with the
driving field aligned along $y$, $\vec{H}||y$, then allow the system to
relaxSilveira and Kardar (1999) at that applied field. We take the angle
$\phi$ of the applied field to be $\phi={\rm Arctan}(H_{y}/H_{x})$, so the
initial direction of the applied field is $\phi=\pi/2$. After rotating the
applied field by an amount $\delta\phi(\vec{H})$, the spin configuration is
updated successively so as to minimize the energy, in the $\omega\rightarrow
0$ limit. After a transient response, the response of the system then settles
into a limit cycle.

Each time the applied field direction is updated, the energy is minimized on
each site by aligning the spin on each site with its effective field,
$\vec{h}^{eff}_{i}$. Hence the following update strategy is repeated until the
spin configuration converges to the nearest energy minimum:

$\displaystyle\vec{h}^{eff}_{i}(t)$ $\displaystyle=$ $\displaystyle
J\sum_{j\in\langle i,j\rangle}\vec{S}_{j}(t)+\vec{h}_{i}+\vec{H},$
$\displaystyle\vec{S}_{i}(t+1)$ $\displaystyle=$
$\displaystyle\frac{\vec{h}^{eff}_{i}(t)}{|\vec{h}^{eff}_{i}(t)|}$ (3)

This update mechanism is similar to Eqn. (2) of Ref. Shukla and Kharwanlang
(2011), however the effective on-site field in our case includes only the
instantaneous influence of nearest neighbors, whereas Ref. Shukla and
Kharwanlang (2011) is working in a mean-field limit. The update algorithm we
employ is described in more detail below, in Sec. V.1.

We continue to allow spins to relax under the influence of Eqns. 3 until a
limit cycle is reached, defined by $\\{\vec{S}_{i}\\}(\phi+2\pi
n)=\\{\vec{S}_{i}\\}(\phi)$. We use the following parameters in our
simulations: $\delta m_{cutoff}=10^{-4}$, $\delta\phi_{max}=2\pi\times
10^{-4}$, $\delta\phi_{min}=2^{-14}\times\delta\phi_{max}$. Hence the
avalanches ($\delta m$) are only well-defined within the precision of the
driving field angle, $\delta\phi_{min}=2\pi\times 6.1\times 10^{-9}$.

###  V.1 Spin Relaxation Method

The rotation of the driving field and subsequent relaxation of the spin
configuration is performed as follows. Starting from an initial spin state
$\\{\vec{S}_{i}\\}(\phi)$ for a given applied field direction $\phi={\rm
Arctan}(H_{y}/H_{x})$ and with $\delta\phi$ initially set to
$\delta\phi=\delta\phi_{max}$:

  1. 1.

Update $\phi\rightarrow\phi+\delta\phi$.

  2. 2.

Use Eqns. 3 to relax the spin configuration.

  3. 3.

If $\delta m>\delta m_{cutoff}$, then:

    1. (a)

If $\delta\phi=\delta\phi_{min}$, accept the new spin configuration and the
new $\phi$, and proceed to Step 1

    2. (b)

Else reject the changes. Set $\delta\phi\rightarrow\delta\phi/2$ and proceed
to Step 1

  4. 4.

Else accept the new spin configuration and the new $\phi$, and:

    1. (a)

If $\delta\phi=\delta\phi_{max}$ or $\delta m\geq\frac{\delta m_{cutoff}}{2}$,
proceed to Step 1.

    2. (b)

Else, set $\delta\phi\to 2\times\delta\phi$ and proceed to Step 1.

###  V.2 Disorder Averages

Table 2 reports the number of disorder configurations used in Figs. 2, 3, and
5.

Table 2: Number of disorder configurations used in Fig. 2 (a) and (b), Fig. 3(d), and Fig. 5 (a-e). Size ($N=L\times L$) | Configurations
---|---
$64\times 64$ | $75$
$80\times 80$ | $60$
$100\times 100$ | $50$
$128\times 128$ | $40$
$160\times 160$ | $30$

##  VI acknowledgements

We acknowledge helpful conversations with D. Campbell, S. Choudhury, G. A.
Csathy, E. Fradkin, T. C. Li, and C. Nayak. E.W.C. and S.B. acknowledge
support from NSF Grant No. DMR-1508236, Department of Education Grant No.
P116F140459, and XSEDE Grant Nos. TG-DMR-180098 and DMR-190014. S.B.
acknowledges support from a Bilsland Dissertation Fellowship. K.A.D.
acknowledges support from NSF Grant No. CBET-1336634. This research was
supported in part through computational resources provided by Information
Technology at Purdue, West Lafayette, Indiana.

## References

  * Lacour-Gayet, P. and Toulouse, G. (1974) Lacour-Gayet, P. and Toulouse, G., J. Phys. France 35, 425 (1974).
  * Aoyama and Ikeda (2006) K. Aoyama and R. Ikeda, Phys. Rev. B 73, 060504 (2006).
  * Wehr  _et al._ (2006) J. Wehr, A. Niederberger, L. Sanchez-Palencia, and M. Lewenstein, Phys. Rev. B 74, 224448 (2006).
  * Ostlund and Halperin (1981) S. Ostlund and B. I. Halperin, Phys. Rev. B 23, 335 (1981).
  * Nelson and Halperin (1979) D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979).
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
  * Young (1979) A. P. Young, Phys. Rev. B 19, 1855 (1979).
  * Lilly  _et al._ (1999) M. P. Lilly, K. B. Cooper, J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett. 82, 394 (1999).
  * Abanin  _et al._ (2007) D. A. Abanin, P. A. Lee, and L. S. Levitov, Phys. Rev. Lett. 98, 156801 (2007).
  * (10) S. M. Girvin and A. H. MacDonald, Perspectives in Quantum Hall Effects: Novel Quantum Liquids in Low-Dimensional Semiconductor Structures, ed. Das Sarma and Pinczuk .
  * Berezinskii (1971) V. Berezinskii, Soviet Physics JETP 32, 493 (1971).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Imry and Ma (1975) Y. Imry and S.-k. Ma, Phys. Rev. Lett. 35, 1399 (1975).
  * Crawford (2013) N. Crawford, EPL (Europhysics Letters) 102, 36003 (2013).
  * Minchau and Pelcovits (1985) B. J. Minchau and R. A. Pelcovits, Phys. Rev. B 32, 3081 (1985).
  * Fel’dman (1999) D. É. Fel’dman, Journal of Experimental and Theoretical Physics 88, 1170 (1999).
  * Yao  _et al._ (2018) N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, ArXiv , 1801.02628 (2018).
  * Shapere and Wilczek (2012) A. Shapere and F. Wilczek, Phys. Rev. Lett. 109, 160402 (2012).
  * Regev  _et al._ (2015) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Nature Communications 6, 47 (2015).
  * Regev  _et al._ (2013) I. Regev, T. Lookman, and C. Reichhardt, Physical Review E 88, 062401 (2013).
  * Leishangthem  _et al._ (2017) P. Leishangthem, A. D. S. Parmar, and S. Sastry, Nature Communications 8, 14653 (2017).
  * Yu  _et al._ (2019) W. C. Yu, J. Tangpanitanon, A. W. Glaetzle, D. Jaksch, and D. G. Angelakis, Physical Review A 99, 033618 (2019).
  * Kivelson  _et al._ (2003) S. A. Kivelson, I. P. Bindloss, E. Fradkin, V. Oganesyan, J. M. Tranquada, A. Kapitulnik, and C. Howald, Reviews of Modern Physics 75, 1201 (2003).
  * Fradkin and Kivelson (2010) E. Fradkin and S. A. Kivelson, Science 327, 155 (2010).
  * Carlson and Dahmen (2011) E. W. Carlson and K. A. Dahmen, Nature Communications 2, 379 (2011).
  * Niederberger  _et al._ (2008) A. Niederberger, T. Schulte, J. Wehr, M. Lewenstein, L. Sanchez-Palencia, and K. Sacha, Physical Review Letters 100, 5 (2008).
  * Sanchez-Palencia and Lewenstein (2010) L. Sanchez-Palencia and M. Lewenstein, Nature Physics 6, 87 (2010).
  * Silveira and Kardar (1999) R. d. Silveira and M. Kardar, Phys. Rev. E 59, 1355 (1999).
  * Shukla and Kharwanlang (2011) P. Shukla and R. S. Kharwanlang, Phys. Rev. E 83, 011121 (2011).