# Simplified algorithms for adaptive experiment design in parameter estimation

Robert D. McMichael<EMAIL_ADDRESS>Sean M. Blakley  National Institute
of Standards and Technology, Gaithersburg, Maryland 20899, USA

###### Abstract

In experiments to estimate parameters of a parametric model, Bayesian
experiment design allows measurement settings to be chosen based on utility,
which is the predicted improvement of parameter distributions due to modeled
measurement results. In this paper we compare information-theory-based utility
with three alternative utility algorithms. Tests of these utility alternatives
in simulated adaptive measurements demonstrate large improvements in
computational speed with slight impacts on measurement efficiency.

##  I Introduction

Bayesian experiment design is a powerful method for adaptive measurements to
estimate the parameters of a nonlinear model function. Measurements to
estimate parameters are commonplace in the physical sciences, where the
traditional approach is to automate measurements using a sequence of pre-
selected settings (a static design), followed by least-squares fitting of a
model function to the data. This measure-then-fit approach is simple and
effective, and entirely appropriate for many applications, but when
measurement resources are limited, efficiency becomes more important.

As a more efficient alternative to static designs, sequential experiment
designs use accumulated data to inform measurement setting decisions. Whether
by choosing optimal designs, or simply by avoiding wasteful efforts,
sequential design generally requires fewer measurements to achieve a given
level of precision. However, sequential design adds two statistical tasks that
add computational cost to the measurement cycle. Each cycle or epoch includes
a setting-decision/design task and a data analysis task in addition to setting
adjustment and data collection.

Optimal Bayesian experiment design has roots in information theory and
decision theory[1] and is well-suited for nonlinear model functions. Chaloner
and Verdinelli[2] published a review of early work that has become a
touchstone in the field, and in recent decades, the availability of
computational power has sparked a resurgence of interest in Bayesian optimal
design. In ref. 3, Overstall et al. provide a brief review of design methods,
and computational methods have been reviewed by Ryan et al.[4]. More recently,
publication of software has provided direct access to Bayesian design methods
methods, potentially reducing the effort required for implementation.[5, 6, 7,
8, 9]

The Bayesian optimal design method centers around locating maxima of the
utility, $U(d)$, which expresses the goals of the measurement as a function of
candidate setting designs $d$. Several authors have identified computation of
the utility as a particularly difficult part of Bayesian experiment design
that prohibits its use.[3, 4, 10, 11] Importantly, for a sequential design to
be preferable to a static design, the cost of implementing and running the
design processes must not exceed the value of saved resources. For
practitioners, utility calculations must be easy and fast.

In this paper, we test the performance of several utility algorithms on two
physics-related parameter estimation problems. The utility algorithms range in
complexity from a numerical evaluation of the information entropy change to a
minimal working example. Surprisingly, we find that these methods achieve
similar results on a per-measurement basis, while computation time varies
widely. The simplest utility algorithm tested here is orders of magnitude
faster than the information-theoretic version with only a slight decrease in
measurement efficiency.

Section II provides background, including an overview of the sequential design
approach, Bayesian inference methods and the information-theoretic utility
function. Section III.1 presents the numerical utility algorithms used in
testing, with results presented in IV, followed by discussion in section V.

##  II Background

In this section, we provide an overview of the statistical methods. Bayesian
experiment design provides a design method based in decision theory,[1] and is
described in a review by Chaloner and Verdinelli. [2]. The measurement run
comprises a sequence of measurement cycles, or epochs, each including a design
step to select measurement settings, a measurement using those settings and
yielding new data, and analysis to incorporate the new data. The approach is
simulation-based Bayesian experiment design[12] follows the examples provided
by Granade et al.[13] and Huan and Marzouk[10]. Below, we describe the data
analysis process followed by the design process, which is the main focus of
the paper.

The design and inference processes rely on a measurement model

$y=f(\theta,d)+\eta,$ (1)

where $y$ is the measurement output, $f(\theta,d)$ is a deterministic,
nonlinear function of parameters $\theta$ and settings $d$ modeling the mean
measurement value. The model function $f(\theta,d)$ may be nonlinear and non-
monotonic, but is assumed to be a well-behaved and accurate model of
measurement results. Random variable $\eta$ has a zero-mean distribution
$P_{\eta}(\eta)$ and models measurement noise. The form of (1) accommodates
intrinsically probabilistic models, such as quantum mechanics. Finally,
$\theta$ is a vector of model parameters
$\theta=\left\\{\theta_{1},\theta_{2},\ldots\right\\}$ treated as random
variables with joint distribution $P(\theta)$. The object of the measurement
run is to provide precise estimates of the parameters, or equivalently, to
narrow $P(\theta)$ around the true parameter values.

The analysis in each design-measure-analyze epoch uses Bayesian inference to
calculate the influence of new data on the parameter distribution. After $n$
measurements, we wish to infer the posterior parameter distribution
$P_{n}(\theta)\equiv P(\theta|{\bm{y}}_{n},{\bm{d}}_{n})$ given accumulated
measurement data ${\bm{y}}_{n}\equiv\\{y_{1},y_{2},\ldots y_{n}\\}$ and
corresponding settings ${\bm{d}}_{n}\equiv\\{d_{1},d_{2},\ldots d_{n}\\}$.
Bold font indicates accumulated data, while non-bold indicates data from a
single epoch. The inference can be performed iteratively using Bayes rule,

$P(\theta|{\bm{y}}_{n},{\bm{d}}_{n})\propto
P(y_{n}|\theta,d_{n})\,P(\theta|{\bm{y}}_{n-1},{\bm{d}}_{n-1})$ (2)

beginning with prior parameter distribution $P_{0}(\theta)$. The likelihood
$P(y_{n}|\theta,d_{n})$, of obtaining result $y_{n}$ is the same as the
probability that the noise value satisfies (1), or

$P(y_{n}|\theta,d_{n})=P_{\eta}(y_{n}-f(\theta,d_{n})).$ (3)

In this paper we will assume that measurement noise is normally distributed
with standard deviation $\sigma_{\eta}$. The likelihood of result $y_{n}$ is
then given by

$P(y_{n}|\theta,d_{n})=\frac{1}{\sqrt{2\pi}\sigma_{\eta}}e^{[y_{n}-f(\theta,d_{n})]^{2}/2\sigma_{\eta}^{2}},$
(4)

providing a quantitative answer to the question “How well do the parameters
explain the data?”

The distribution of parameters $P(\theta)$ encapsulates the state of belief
regarding different parameter values. For computational work, methods known
variously as particle filters, swarm filters or sequential Monte Carlo methods
have gained prominence as representations of the probability distribution.[14,
15, 16, 17] For a model with $D$ parameters, the distribution is represented
by $N_{\rm p}$ “particles” each with a vector $\theta$ corresponding to a
point in $D$-dimensional parameter space and a corresponding weight $w$ such
that the distribution is approximated by

$P(\theta)\approx\sum_{k=1}^{N_{\rm p}}\delta(\theta-\theta_{k})w_{k}$ (5)

with $\sum_{i}w_{i}=1$. Here $\delta(\cdot)$ is the Dirac delta function. In
inference calculations, unnormalized weights $W_{i}$ are generated by

$W_{i,n}=w_{i,n-1}P(y_{n}|\theta_{i},d_{n})$ (6)

yielding new weights

$w_{i,n}=W_{i,n}/\sum_{i}W_{i,n}.$ (7)

Resampling is a necessary maintenance task for particle filters that
effectively reassigns computer memory from low-probability regions of
parameter space to high-probability regions. Resampling is initiated whenever
the effective number of particles,

$N_{\rm eff}=1/\sum w_{i}^{2},$ (8)

falls lower than a threshold, typically set at $0.5N_{\rm p}$. The covariance
matrix ${\bm{C}}$ of the distribution is calculated for a later step, then
$N_{\rm p}$ samples $\tilde{\theta}_{j}$ are drawn with replacement from the
$\\{\theta_{i}\\}$ set with probability $w_{i}$. This selection process will
tend to miss low-weight particles, while possibly choosing high-weight
particles multiple times. The particles are then given small random
displacements $\delta\theta$ to separate degenerate particles. The resampled
set of particles $\theta_{j}$ are given by

$\theta_{j}=\tilde{\theta}_{j}+\delta\theta_{j}$ (9)

where $\delta\theta_{j}$ is a sample from a multivariate normal distribution
with covariance matrix $\alpha{\bm{C}}$. We use $\alpha=0.01$.

We turn next to the question of selecting a design $d_{n+1}$ for the next
measurement from a finite number of setting combinations. The information-
theoretic approach introduced by Lindley[1] and followed by many others
chooses information entropy as a measure of the parameter distribution quality
and it follows that the predicted change in information entropy following a
measurement is the utility, $U(d^{\prime})$. The entropy ${\cal H}$ of a
distribution $P(\cdot)$ is given by

${\cal H}_{x}[P]=-\int P(x)\log\left(P(x)\right)\,dx.$ (10)

The Kullback-Leibler divergence describing the change in information entropy
in going from $P_{n}(\theta)$ to $P_{n+1}(\theta)$ is

$KL=-\int
P_{n+1}(\theta)\log\left[\frac{P_{n}(\theta)}{P_{n+1}(\theta)}\right]d\theta,$
(11)

where $P_{n}(\theta)$ is the parameter distribution after $n$ epochs, and
$P_{n+1}(\theta)$ is the projected distribution after getting a possible
measurement result $y^{\prime}$ using a candidate setting design $d^{\prime}$.
The “prime” notation indicates speculative values as opposed to values $d$ and
$y$ which are measurement records. Using Bayes rule, the resulting parameter
distribution would be given by

$\displaystyle P_{n+1}(\theta)$ $\displaystyle\equiv$ $\displaystyle
P(\theta|y^{\prime},d^{\prime},{\bm{y}}_{n},{\bm{d}}_{n})$ $\displaystyle=$
$\displaystyle\frac{P(y^{\prime}|\theta,d^{\prime})}{P(y^{\prime})}P(\theta|{\bm{y}}_{n},{\bm{d}}_{n})=\frac{P(y^{\prime}|\theta,d^{\prime})}{P(y^{\prime})}P_{n}(\theta),$

and the Kullback-Leibler divergence would then be written as

$KL(d^{\prime},y^{\prime})=-\int\frac{P(y^{\prime}|\theta,d^{\prime})}{P(y^{\prime})}P_{n}(\theta)\log\left[\frac{P(y^{\prime})}{P(y^{\prime}|\theta,d^{\prime})}\right]d\theta.$
(13)

Averaging over the possible $y^{\prime}$ measurement values yields the utility
$U(d^{\prime})$

$\displaystyle U(d^{\prime})$ $\displaystyle=$ $\displaystyle-\int
P(y^{\prime}|d^{\prime})\log\left[P(y^{\prime}|d^{\prime})\right]dy^{\prime}$
$\displaystyle+\int P_{n}(\theta)\left\\{\int
P(y^{\prime}|\theta,d^{\prime})\log\left[P(y^{\prime}|\theta,d^{\prime})\right]dy^{\prime}\right\\}d\theta$

with

$P(y^{\prime}|d^{\prime})=\int
P(y^{\prime}|\theta,d^{\prime})P_{n}(\theta)\,d\theta.$ (15)

Equation (II) is the conventional, information-theoretic utility of a
candidate setting, and it is the starting point for the investigations in this
paper. This is the utility expression we approximate and/or emulate with
simplified expressions.

The utility (II) has an intuitive interpretation. In the second term of (II)
the distribution $P(y^{\prime}|\theta,d^{\prime})$ appearing in the bracketed
integral is the forecast distribution of measured values given a set of
parameters $\theta$ and design $d^{\prime}$. But with $\theta$ and
$d^{\prime}$ fixed, the arguments of the model function are all fixed and the
distribution of $y^{\prime}$ values is the noise distribution, offset by the
model function value. Therefore, the second term may be interpreted as the
(negative) entropy of the measurement noise, averaged over parameters.

The first term in (II) is also an entropy of the predicted measurement value
distribution, but here the distribution includes noise and also the effects of
the parameter distribution as written explicitly in (15). In total, the
utility is the entropy of possible measurement outputs produced by parameter
uncertainty and noise, discounted by the entropy of noise.

The intuitive interpretation of (II) is simply that useful settings are those
where the spread or dispersion of the parameter distribution produces large
variations in measurement predictions relative to the measurement noise. In
the following section we present utility variations that all share this
interpretation, but differ in their metrics for dispersion.

##  III Methods

###  III.1 Utility Expressions

This subsection describes utility algorithms based on, or at least inspired
by, (II), with the goal of practical application in laboratory settings. To be
useful, design algorithms must produce benefits of measurement efficiency that
outweigh the cost of calculation. To meet that goal, we place a heavy emphasis
on simplicity and speed.

In pursuit of simplicity and speed, we also permit ourselves to invent trial
utility functions without rigorous derivation. To justify this approach, we
argue that the penalty for less-than-perfect design choices is a mild decrease
in measurement efficiency, without invalidating the measurement results.
Acceptable measurements are frequently made using completely arbitrary
settings, after all. So, while we relax requirements for rigor in the
formulation of utility algorithms, we are careful to monitor precision and
accuracy of the parameter estimation results in the results presented below.

Qualitatively, the utility given by (II) suggests that the best measurements
will be made where the dispersion of the parameter distribution has a large
effect on measurement outcomes relative to the dispersion of measurement
noise. In (II), the dispersion measure is information entropy, made so by the
conventional choice of information entropy to gauge improvements in
$P(\theta)$. With the algorithms below, we explore alternative dispersion
measures in utility functions.

####  III.1.1 KLD utility algorithm

This calculation attempts to approximate Eq. (II). For the first term, $N_{\rm
s}$ samples are drawn from the prior $P_{n}(\theta)$ and from the noise
distribution $P(\eta)$ to simulate measurement outcomes using the measurement
model, (1). The simulated outcomes are samples from
$P(y^{\prime}|\theta,d^{\prime})$. The differential entropy is estimated from
these samples using either Ebrahimi’s method[18] or Vasicek’s method[19]. For
the examples in this paper, the noise is normally distributed with standard
deviation $\sigma_{\eta}$, independent of $\theta$ and $d$, so the second term
in (II) can be determined from the known properties of the noise distribution:
${\cal H}_{\eta}=\frac{1}{2}\log(2\pi e\sigma_{\eta}^{2})$.

Pseudocode for $U^{\rm KLD}(d^{\prime})$ is provided in Algorithm 1 of the
appendix.

####  III.1.2 Variance algorithm

For the variance algorithm, we choose the logarithm of variance as a measure
of dispersion[20], yielding the utility

$U^{\rm
var}(d^{\prime})=\frac{1}{2}\log\left[v_{\theta}(d^{\prime})+v_{\eta}\right]-\frac{1}{2}\log(v_{\eta}),$
(16)

which simplifies to

$U^{\rm
var}(d^{\prime})=\frac{1}{2}\log\left[1+\frac{v_{\theta}(d^{\prime})}{v_{\eta}}\right].$
(17)

Here, $v_{\eta}$ is the variance of the noise and $v_{\theta}(d^{\prime})$ is
the variance of a distribution of noise-free model values,
$P(y^{\prime}|d^{\prime},\eta=0)$. Samples of this distribution are provided
by the model function $f(\theta_{i}|d^{\prime})$ using parameter samples
$\\{\theta_{1}\ldots\theta_{N_{s}}\\}$. The first term in (16) represents the
combined effects of parameters and noise, and the second term isolates the
contribution due to noise alone. The variance utility also corresponds to the
Kullback-Leibler utility (II) in the special case where all distributions are
normal. For a normal distribution with standard deviation $\sigma$, the
entropy is

${\cal H}=\frac{1}{2}\log(2\pi e\sigma^{2}),$ (18)

where $\log(e)=1$.

A feature of the variance algorithm is that the variances of the parameter
distribution and noise distributions are simply additive, so that effects of
measurement noise and effects of parameter dispersion are separated. In
contrast, the contributions to entropy are not easily distinguished. Another
advantage is that variance computes faster than entropy.

Pseudocode for $U^{\rm var}(d)$ is provided in Algorithm 2 of the appendix.

####  III.1.3 Pseudo-utility algorithm

Recognizing that standard deviation is not always a good proxy for entropy,
the pseudo-utility attempts to recapture some of the properties of entropy as
a dispersion measure while preserving the separation of parameter effects and
noise effects provided by variance. With the same form as the variance
utility, the pseudo-utility is defined as

$U^{\rm psu}(d^{\prime})=\frac{1}{2}\log\left[1+\frac{v_{\cal
H}(d^{\prime})}{v_{\eta}}\right],$ (19)

where $v_{\cal H}$ is an effective variance derived from the entropy of the
distribution of noise-free model values. As in the variance algorithm, model
values are generated by the model function $f(\theta_{i}|d^{\prime})$ for
parameter samples $\\{\theta_{1}\ldots\theta_{N_{s}}\\}$ and $\eta=0.$ The
effective variance $v_{\cal H}$ requires the entropy[18, 19] of the noise-free
model value distribution ${\cal H}_{y^{\prime}}(d^{\prime},\eta=0)$ which is
then then transformed into a variance, using (18), which yields

$v_{\cal H}=(2\pi e)^{-1}\exp[2{\cal H}_{y^{\prime}}(d^{\prime},\eta=0)].$
(20)

The pseudo-utility is pieced together without a strong claim for validity, but
it has been used in previous work.[21]

Pseudocode for $U^{\rm psu}(d)$ is provided in Algorithm 3 of the appendix.

####  III.1.4 Max-Min utility algorithm

The max-min utility is designed for simplicity and fast computation using the
range statistic as a measure of dispersion. A relatively small set of samples
$\theta_{j}$, $\\{\theta_{1}\ldots\theta_{N_{s}}\\}$ are drawn from
$P(\theta)$ and a corresponding set of $N_{s}$ model values
$\\{y(d^{\prime})\\}\equiv\\{y_{1}(d^{\prime})\ldots y_{N_{s}}(d^{\prime})\\}$
is calculated using $y_{j}(d^{\prime})=f(\theta_{j},d^{\prime})$. The range
$t(d^{\prime})$ is defined as the difference between maximum and minimum
values in the set,

$t(d^{\prime})={\rm max}(\\{y(d^{\prime})\\})-{\rm min}(\\{y(d^{\prime})\\}).$
(21)

Again following the form of the variance utility, the max-min utility is
defined as

$U^{\rm
mm}(d^{\prime})=\frac{1}{2}\log\left[1+\frac{t(d^{\prime})^{2}}{v_{\eta}}\right].$
(22)

####  III.1.5 Random design

As its name suggests, the random design chooses randomly from among the
candidate settings, uninfluenced by measurement data. Because the settings are
chosen with uniform probability, the measurement results are expected to
exhibit the same overall performance as repeated sweeps of the setting value.
However, the random approach avoids artifacts that might be generated by
periodic repetition of the settings.

###  III.2 Sample Reuse

Improved methods have been proposed by several authors to compute the double
integral in the utility expression (II) with sufficient precision [10, 11, 22,
23, 24]. In order to serve its purpose in experiment design, an approximate
utility function need only exhibit a maximum at at design that is close to the
design that maximizes the true utility. The utility values are ultimately
discarded, while the design determines the next measurement. To simplify
calculation of (II), Huan and Marzouk[10] proposed using the same set of
parameter samples both for calculation of (15) and for the averaging over
parameters in (II). Where possible, we take this approach one step further by
using one set of random draws to calculate utility for all candidate settings.
Some modifications would be needed if the noise distribution independent on
$\theta$ or $d$.

Fig. 1 illustrates how sample reuse affects the utility function. Panels a)
and b) plot samples from $P(y^{\prime}|x)$ for a Lorentzian model function

$y^{\prime}=\frac{1}{[(x-x_{0})/0.1]^{2}+1}+\eta$ (23)

with the peak center parameter $x_{0}$ and added noise $\eta$ both normally
distributed. Panel a) plots simulated measurement outcomes where $N_{s}$
parameter samples and $N_{s}$ noise samples are drawn for each of $N_{d}$
candidate designs, requiring $2N_{\rm s}\cdot N_{\rm d}$ samples $[N_{\rm
s}=100;N_{\rm d}=200]$. Fig. 1b) represents the same distribution, but single
sets of $N_{\rm s}$ $x_{0}$ samples and $\eta$ samples are reused for all
values of design variable $x$, requiring only $2N_{s}$ samples. Fig. 1c) shows
the effect of these contrasting approaches on the utility function using the
KLD algorithm described below. Using new samples for each setting value
ensures that the utility values are statistically independent, but introduces
sampling noise that makes the maximum difficult to locate.[4] By reusing
samples, the sampling noise is converted to a systematic error, and the point-
to-point variation is suppressed. Both formulations will yield near-optimal
setting values, but reusing samples greatly reduces computation time. Unless
otherwise indicated we use $N_{\rm s}=1000$ samples.

Figure 1: Illustration of the effects of parameter and noise sample reuse on
utility function estimation. The model function is a Lorentzian given by (23).
Panels a) and b) illustrate distributions of predicted measurement values
$P(y^{\prime}|x)$ using sets of $N_{\rm s}=100$ samples each from normal
distributions of center $x_{0}$ and added noise $\eta$. Standard deviations
are indicated by bars. In a), a fresh set of samples is drawn for each value
of setting $x$, while in b) a single set of samples is reused for all $x$. c)
Corresponding estimates of utility, both yielding maxima near $x=2.6$.

##  IV Results and discussion

###  IV.1 Lorentzian peak

In this subsection, the simulated experiment is a Lorentzian absorption
profile, with a single “unknown” peak-center parameter
$\theta\equiv\\{x_{0}\\}$ and setting $d\equiv{x}$

$f_{\rm L}(x_{0},x)=b+\frac{a}{[(x-x_{0})/\Delta]^{2}+1}.$ (24)

Constants are background $b$, amplitude $a$, and half-width $\Delta$. The
maximum Fisher information for the Lorentzian is found at
$x=x_{0}\pm\Delta/\sqrt{3}$, and the Cramer-Rao bound is given by

$\sigma_{x_{0}}\geq\frac{8}{3\sqrt{3}}\frac{\Delta}{a}\frac{\sigma_{\eta}}{\sqrt{n}}.$
(25)

The Lorentzian function describes the frequency response of a damped harmonic
oscillator, which appears frequently in many branches of science and
engineering. Previous work demonstrated the effectiveness of sequential
Bayesian experiment design on a similar system, estimating 5 parameters from
measurements on triplets of similar peaks. Compared to a swept-setting
approach, the Bayesian method required 40-fold fewer measurements to locate
and measure to a comparable uncertainty.[25]

Figure 2: Simulated measurement runs to determine the center of a Lorentzian.
a) The first 300 simulated data points of a selected run. The Lorentzian model
function is plotted in red with the center parameter $x_{0}=2.6$, and the
standard deviation of the noise is equal to the peak height. The adaptive
design initially calls for settings over a broad range, but later repeats
measurements on the sides of the peak after $\approx 100$ epochs. b) Standard
deviation and c) entropy of $P(x_{0})$ calculated using various utility
algorithms. Mean (solid lines) and 5 % to 95 % credibility intervals (shaded
areas) are calculated from 400 runs. The grey areas result from coincidence of
the credibility intervals. d) Root-mean-square (RMS) error. Variance, KLD and
pseudo algorithms use $N_{s}=1000$ parameter samples; max-min uses $N_{\rm
s}=2$.

Fig. 2 presents results of simulated measurement runs using the sequential
Bayesian design methods along with random setting selection for comparison.
The simulations use constants $b=50000$, $a=-1000$, and $\Delta=0.1$.
Measurement data is simulated using true value $\theta^{*}=\\{x_{0}=2.6\\}$
and measurement noise is normally distributed with $\sigma_{\eta}=1000$. The
prior $P_{0}(x_{0})$ is normal with mean 3.0 and standard deviation 0.5. Fig.
2a) shows data from a single run with the model “true” curve. After
approximately 90 epochs in this run, the design algorithm focuses on settings
close to the Fisher information maxima on the slopes of the peak/dip.

Panels b-d) of Fig. 2 support some of the main conclusions of this paper. In
summaries of 400 measurement runs, panels b), c) and d) respectively trace the
standard deviation and entropy of the parameter distribution, and the RMS
error of the mean.

A striking feature of Fig. 2 is that virtually identical results are obtained
regardless of which utility algorithm is used. The sequential designs all
follow similar approaches to the Cramer-Rao limit (25) after a few hundred
measurements. Compared to uninformed Random setting selection, all adaptive
designs achieve $\sigma_{x_{0}}$ values that are a factor of 4 lower,
corresponding to an order-of-magnitude savings in the number of measurements
needed.

Figure 3: For each of the utility algorithms, setting values from 50 runs are
superposed, demonstrating similar overall patterns in design choices. Epochs
$>50$, are sub-sampled for uniform appearance on the logarithmic scale.

Similarities in the the results are underscored by similarities in the design
settings, as shown in Fig. 3. Each run begins by forming a new particle filter
using independent samples from the prior distribution, and each data
simulation uses independent noise values, so different design choices are
expected in different runs. However, patterns in the design choices have
strong similarities across the various utility algorithms, with an initial
focus on $x$ values close to the center of the prior distribution at
$x_{0}=3.0$, a broadening of the setting choices followed by rapid convergence
to the settings for Fisher information maxima.

Figure 4: Tests of the max-min algorithm for varying numbers of parameter
sample $N_{\rm s}$ including a) standard deviation and b) entropy of the
$x_{0}$ parameter distribution. Performance approximates the KLD algorithm for
small $N_{\rm s}$.

Figure 4 highlights behavior of the max-min utility algorithm and its
dependence on the number of parameter samples $N_{\rm s}$ drawn for design
calculations. The max-min algorithm uses the range of model function outputs
over the parameter samples as a crude measure of dispersion due to parameter
uncertainty. The results show that the minimum number of parameter samples
($N_{\rm s}=2$) performs as well or better than larger sample numbers, and
that its performance approaches that of the full KLD calculation.

Figure 5: Comparison of design computation times among the utility algorithms
in 1000-epoch runs. The runs simulate measurements to determine the center,
$x_{0}$, of a Lorentzian peak using 200 candidate settings. Overall bar length
indicates total design computation, with sub-tasks shown as segments. Model
function and entropy calculations consume the majority of the total time. The
max-min algorithm with $N_{\rm s}=2$ executes very quickly.

Figure 5 shows the computation time required for the design algorithms and
selected sub-tasks. Each of the design calculations requires computation of
the model function using $N_{d}=200$ candidate designs for each of $N_{\rm
s}=1000$ parameter samples, followed by $N_{d}$ calculations of a statistic
over $N_{s}$ samples. The exception is the max-min algorithm using only
$N_{\rm s}=2$ parameter samples, which is unsurprisingly much faster.

###  IV.2 Decaying sinusoid

This subsection presents simulations of measurements to determine four
parameters of a decaying sinusoid model function:

$f(\\{h,c,\omega_{0},T_{2}\\},\tau)=h+c\;\sin(\omega_{0}\tau)\;e^{-(\tau/T_{2})^{2}},$
(26)

where $\theta\equiv\\{h,c,\omega_{0},T_{2}\\}$ are “unknown” parameters
corresponding respectively to a background level, contrast, mean angular
frequency, and dephasing with true values
$h^{*}=0.8,c^{*}=0.13,\omega_{0}^{*}=9.4$ and $T_{2}^{*}=10$. The lone setting
is $d\equiv\tau$. For each measurement, the system is initialized, the
oscillations are started at $t=0$ and a single measurement is allowed after
the a selected delay time $\tau$. For optimal measurements of $\omega_{0}$,
the maximum slope $df/d\omega_{0}$ occurs for $\tau\approx T_{2}/\sqrt{2}$ and
$\sin(\omega_{0}\tau)\approx 0$. For $n$ repeated measurements at this
setting,

$\sigma_{\omega_{0}}=\frac{\sqrt{2e}}{cT_{2}}\frac{\sigma_{y}}{\sqrt{n}}.$
(27)

For optimal measurements of $T_{2}$, the maximum slope $df/dT_{2}$ occurs for
$\tau\approx T_{2}$ and $\sin(\omega_{0}\tau)\approx\pm 1$.

The decaying sinusoid describes a Ramsey measurement, which is the canonical
method for measuring energy differences $\Delta E=\hbar\omega_{0}$ between
quantum states[26], and it is closely related to the problem of measuring
quantum phase.

In (27), the frequency uncertainty is limited by measurement noise, but even a
single noise-free measurement would yield ambiguous results since many values
of $\omega_{0}$ can produce the same value of $y$. Interest in Ramsey
measurements is driven by the possibility that the energy difference $\Delta
E$ can be measured (i.e., disambiguated) with precision $\Delta E\propto
1/\Delta t$ in a time $\Delta t$. [27, 28] This behavior is termed Heisenberg
scaling in analogy with the Heisenberg uncertainty principle, $\Delta E\Delta
t\geq\hbar/2$. When readout fidelity is high, but not perfect, Heisenberg
scaling can be approached using adaptive[29, 30, 31, 32, 33] and by optimized
non-adaptive[34, 35, 36] measurements.

Figure 6: Simulated measurement run to determine four parameters of a decaying
sinusoid. a) The model function is plotted in red as a function of setting
$\tau$. The standard deviation of the noise is equal to the oscillation
amplitude at $\tau=0$. The adaptive design initially chooses low $\tau$
settings, but later favors $\tau\approx T_{2}$. Data from epochs later than 50
have been subsampled for clarity. b) Standard deviation and c) entropy of
$P(\omega_{0})$ of simulated measurements using various utility algorithms.
Mean (solid lines) and 5 % to 95 % credibility intervals (shaded areas) are
calculated from 400 runs. d) RMS error.

Figure 6 presents results of Ramsey sequence simulations. Fig. 6a) plots the
“true” model function (26) and simulated data. Allowed settings cover the
range $0.1\leq\tau\leq 20$ in steps of 0.01. The simulated noise is normally
distributed, with $\sigma_{\eta}=0.13$. Fig. 6b) compares the standard
deviation of $\omega_{0}$ for various utility algorithms as a function of the
number of measurement epochs. As with the results of the Lorentzian peak
example, the utility algorithms achieve similar performance on a per-
measurement basis.

Figure 7: Setting values plotted for 50 runs each using the indicated utility
algorithms. In later epochs, the algorithms choose small $\tau$ and also
$\tau$ distributed across ranges, for example the range indicated by the
dashed rectangle in panel a). The chosen $\tau$ values tend to fall at local
extrema of the model function and local extrema of slope. These correlations
are more clearly visible in Fig. 6a)

The setting values generated by the algorithms are plotted in Fig. 7. Strong
similarities between the designs are visible, despite differences in the
algorithms. All the design algorithms begin with low-$\tau$ measurements and
increase $\tau$ gradually. After 100 epochs, the algorithms select
concentrated measurements at low $\tau$ and distributed measurements in the
region roughly $5\lesssim\tau\lesssim 15$, for example within the boxed region
in Fig. 7a).

Within the diffuse region of settings, there are strong correlations between
$\tau$ values and the sinusoid phase. In Fig. 6a), close examination shows
that the low-$\tau$ settings are correlated with maxima and minima of the
model function. For $8\lesssim\tau\lesssim 11$, $\tau$ values are correlated
with zero crossings, suggesting optimal settings for determining $\omega_{0}$.
Also, $\tau$ above 11 are correlated with maxima and minima suggesting optimal
settings for determining $T_{2}$.

In contrast to the Lorentzian model results in Fig. 2, where there are large
performance differences between the adaptive algorithms and random parameter
selection, the random design applied to the Ramsey model achieves similar
performance to adaptive designs after approximately 400 epochs. We attribute
the relatively good performance of random setting selection here to the
feature-rich structure of the Ramsey model function and the larger number of
parameters. Taken together, these factors create a situation where virtually
any selection of $\tau$ is somewhat informative. The broad distributions of
settings chosen by adaptive designs also suggest that utility is not
concentrated on a few settings in the decaying sinusoid measurement.

Figure 8: Comparison of design task computation among the utility algorithms
in 1000-epoch runs using the decaying sinusoid model with four parameters and
1991 setting values. Within the total time, model function and entropy
calculations are the most costly parts of the algorithms. The max-min
algorithm with $N_{\rm s}=2$ executes very quickly.

Figure 8 shows the computation time required for the design algorithms and
selected sub-tasks. Each of the design calculations requires computation of
the model function using $N_{d}=1991$ candidate designs for each of $N_{\rm
s}=1000$ parameter samples, followed by $N_{d}$ calculations of a statistic
over $N_{s}$ samples. The exception is the max-min algorithm using only
$N_{\rm s}=2$ parameter samples.

##  V Discussion

This paper compares results of simulated measurement runs using two nonlinear
experiment models and four utility algorithms using different statistical
functions, but all treating the utility as a measure of the dispersion of
forecast measurement results relative to the dispersion of measurement noise.
In the regime where signal contrast was on par with measurement noise, the
five utility algorithms generated similar design choices and similar
reductions in parameter variance, entropy, and RMS error. While the KLD
utility has a solid theoretical basis, the KLD algorithm does not produce
dramatically better designs than simplified algorithms.

We point to the max-min algorithm with $N_{\rm s}=2$ as an important outcome
of this study. This algorithm pushes the utility-as-dispersion idea to its
simplest extreme, calculating the model function for each candidate design and
only two parameter samples. This simplified utility algorithm calculates much
faster than the other tested algorithms with only a slight decrease in
efficiency.

An interesting point is that the max-min algorithm shows poorer performance
with larger numbers of parameter samples. Fig. 4 shows that using $N_{\rm
s}>2$ makes less efficient design decisions, despite more thorough sampling of
parameters. We explain the relative inefficiency of $N_{\rm s}>2$ by noting
that the max-min algorithm disregards representative parameter values from the
center of a distribution when $N_{\rm s}>2$, and instead takes direction from
the atypical extreme values.

By reusing parameter samples in utility calculations, we make utility
estimates correlated across designs. The benefits of this technique are first,
that the utility estimates become smooth functions of the setting value with
easily located maxima, and second, that fewer random parameter samples are
required. Although these benefits come at the cost of a systematic error, the
effects of the systematic error are transient, because new parameter samples
are drawn for each epoch. The $N_{\rm s}=2$ max-min algorithm especially
benefits from suppressed sampling noise, because it lacks the natural noise
suppression provided by large sample numbers.

In this paper, we have demonstrated speed and quality of the max-min $N_{\rm
s}=2$ algorithm on CPU hardware. We speculate that this algorithm will be well
suited for implementation on FPGAs or other dedicated hardware with fixed-
point arithmetic, because the max-min statistic requires only subtraction and
absolute value operations. However, practical implementation of a dedicated
system would still face challenges in implementing a parameter probability
density representation and mathematical operations to support model function
and likelihood calculations.

##  Appendix A Pseudocode

Algorithm 1 Kullback-Leibler algorithm

Parameter distribution $P(\theta)$

Noise distribution $P(\eta)$

for $j\leftarrow 1\ldots N_{s}$ do

Sample $\theta_{j}$ from $P(\theta)$ parameter distribution

Sample $\eta_{j}$ from $P(\eta)$ from noise distribution

for all candidate settings $d_{i}$ do

Simulate measurement $y_{i,j}\leftarrow{\rm model}(\theta_{j},d_{i})+\eta_{j}$

$\triangleright$ Reusing $\theta_{j}$ and $\eta_{j}$

end for

end for

for all candidate settings $d_{i}$ do

Estimate entropy from samples $h_{y,i}={\cal H}(y_{i,1}\ldots y_{i,N_{s}})$

Estimate entropy of noise distribution $h_{\eta,i}\leftarrow{\cal H}[P(\eta)]$

Calculate utility $U^{\rm KLD}_{i}\leftarrow h_{y,i}-h_{\eta,i}$

end for

Find max utility $i_{\rm best}\leftarrow$ argmax$(U^{\rm KLD}_{i})$

corresponding setting $d_{i_{\rm best}}$

Algorithm 2 Variance utility algorithm

Parameter distribution $P(\theta)$

for $N_{s}$ parameter samples $\theta_{j}$ drawn from $P(\theta)$ do

for all candidate designs, $d_{i}$ do

Evaluate model $y_{i,j}\leftarrow f(\theta_{j},d_{i})$

$\triangleright$ reusing $\theta_{j}$ samples

end for

end for

for all candidate designs $d_{i}$ do

Variance $v_{\theta,i}\leftarrow{\rm Var}(y_{i,1}\ldots y_{i,N_{s}})$ over
parameters

Variance of noise $v_{\eta,i}\leftarrow{\rm Var}[P(\eta)]$

Calculate utility $U^{\rm var}_{i}\leftarrow\log[1+v_{\theta,i}/v_{\eta,i}]/2$

end for

Find max utility $i_{\rm best}\leftarrow{\rm Argmax}(U^{\rm var}_{i})$

corresponding setting $d_{i_{\rm best}}$

Algorithm 3 Pseudo-utility algorithm

Parameter distribution $P(\theta)$

for $N_{s}$ parameter samples $\theta_{j}$ drawn from $P(\theta)$ do

for all candidate designs, $d_{i}$ do

Evaluate model $y_{i,j}\leftarrow f(\theta_{j},d_{i})$

$\triangleright$ reusing $\theta_{j}$ samples

end for

end for

for all candidate designs $d_{i}$ do

Estimate entropy from samples $h_{y,i}={\cal H}(y_{i,1}\ldots y_{i,N_{s}})$

Calculate effective variance $v_{\theta,i}^{\rm
eff}\leftarrow\exp[2h_{y,i}]/(2\pi e)$

Variance of noise distribution $v_{\eta,i}\leftarrow{\rm Var}[P(\eta)]$

Calculate utility $U^{\rm PU}_{i}\leftarrow\log[1+v_{\theta,i}^{\rm
eff}/v_{\eta,i}]/2$

end for

Find max utility $i_{\rm best}\leftarrow{\rm Argmax}(U^{\rm PU}_{i})$

corresponding setting $d_{i_{\rm best}}$

Algorithm 4 Max-min utility

Parameter distribution $P(\theta)$

for $N_{s}$ parameter samples $\theta_{j}$ drawn from $P(\theta)$ do

for all candidate designs, $d_{i}$ do

Evaluate model $y_{i,j}\leftarrow f(\theta_{j},d_{i})$

$\triangleright$ reusing $\theta_{j}$ samples

end for

end for

for all candidate designs $d_{i}$ do

Find maximum model value $y_{i}^{\rm max}={\rm max}(y_{i,1}\ldots
y_{i,N_{s}})$

Find minimum model value $y_{i}^{\rm min}={\rm min}(y_{i,1}\ldots
y_{i,N_{s}})$

Calculate range $t_{i}\leftarrow y_{i}^{\rm max}-y_{i}^{\rm min})$

Variance of noise $v_{\eta,i}\leftarrow{\rm Var}[P(\eta)]$

Calculate utility $U^{\rm mm}_{i}\leftarrow\log[1+t_{i}^{2}/v_{\eta,i}]/2$

end for

Find max utility $i_{\rm best}\leftarrow{\rm Argmax}(U^{\rm mm}_{i})$

corresponding setting $d_{i_{\rm best}}$

## References

  * Lindley [1956] D. V. Lindley, Annals of Mathematical Statistics 27, 986 (1956).
  * Chaloner and Verdinelli [1995] K. Chaloner and I. Verdinelli, Statistical Science 10, 273 (1995).
  * Overstall, McGree, and Drovandi [2018] A. M. Overstall, J. M. McGree, and C. C. Drovandi, Statistics and Computing 28, 343 (2018).
  * Ryan  _et al._ [2016] E. G. Ryan, C. C. Drovandi, J. M. McGree, and A. N. Pettitt, International Statistical Review 84, 128 (2016).
  * Overstall, Woods, and Adamou [2020] A. M. Overstall, D. C. Woods, and M. Adamou, Journal of Statistical Software 95 (2020), 10.18637/jss.v095.i13.
  * Olofsson  _et al._ [2019] S. Olofsson, L. Hebing, S. Niedenführ, M. P. Deisenroth, and R. Misener, Computers & Chemical Engineering 125, 54 (2019).
  * McMichael, Blakley, and Dushenko [2021] R. D. McMichael, S. M. Blakley, and S. Dushenko, Journal of Research of the National Institute of Standards and Technology 126, 126002 (2021).
  * Ha [2018] F. Ha, ACS Central Science , 12 (2018).
  * Granade  _et al._ [2017] C. Granade, C. Ferrie, I. Hincks, S. Casagrande, T. Alexander, J. Gross, M. Kononenko, and Y. Sanders, Quantum 1, 5 (2017), arXiv: 1610.00336.
  * Huan and Marzouk [2013] X. Huan and Y. M. Marzouk, Journal of Computational Physics 232, 288 (2013).
  * Long  _et al._ [2013] Q. Long, M. Scavino, R. Tempone, and S. Wang, Computer Methods in Applied Mechanics and Engineering 259, 24 (2013).
  * Müller, Sansó, and De Iorio [2004] P. Müller, B. Sansó, and M. De Iorio, Journal of the American Statistical Association 99, 788 (2004).
  * Granade  _et al._ [2012] C. E. Granade, C. Ferrie, N. Wiebe, and D. G. Cory, New Journal of Physics 14, 103013 (2012).
  * Gordon, Salmond, and Smith [1993] N. Gordon, D. Salmond, and A. Smith, IEE Proceedings F Radar and Signal Processing 140, 107 (1993).
  * Arulampalam  _et al._ [2002] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, IEEE Transactions on Signal Processing 50, 174 (2002).
  * Carpenter, Clifford, and Fearnhead [1999] J. Carpenter, P. Clifford, and P. Fearnhead, IEE Proceedings - Radar, Sonar and Navigation 146, 2 (1999).
  * Elfring, Torta, and van de Molengraft [2021] J. Elfring, E. Torta, and R. van de Molengraft, Sensors 21, 438 (2021).
  * Ebrahimi, Pflughoeft, and Soofi [1994] N. Ebrahimi, K. Pflughoeft, and E. S. Soofi, Statistics & Probability Letters 20, 225 (1994).
  * Vasicek [1976] O. Vasicek, Journal of the Royal Statistical Society. Series B (Methodological) 38, 54 (1976), publisher: [Royal Statistical Society, Wiley].
  * Solonen, Haario, and Laine [2012] A. Solonen, H. Haario, and M. Laine, Journal of Computational and Graphical Statistics 21, 234 (2012).
  * McMichael, Dushenko, and Blakley [2021] R. D. McMichael, S. Dushenko, and S. M. Blakley, Journal of Applied Physics 130, 144401 (2021).
  * Drovandi and Tran [2018] C. C. Drovandi and M.-N. Tran, Bayesian Analysis 13 (2018), 10.1214/16-BA1045.
  * Beck  _et al._ [2020] J. Beck, B. M. Dia, L. F. R. Espath, and R. Tempone, International Journal for Numerical Methods in Engineering 121, 3402 (2020), arXiv: 1811.11469.
  * Ryan [2003] K. J. Ryan, Journal of Computational and Graphical Statistics 12, 585 (2003).
  * Dushenko, Ambal, and McMichael [2020] S. Dushenko, K. Ambal, and R. D. McMichael, Physical Review Applied 14, 054036 (2020).
  * Degen, Reinhard, and Cappellaro [2017] C. Degen, F. Reinhard, and P. Cappellaro, Reviews of Modern Physics 89, 035002 (2017).
  * Kitaev [1995] A. Y. Kitaev, arXiv:quant-ph/9511026 (1995), arXiv: quant-ph/9511026.
  * Kitaev [1996] A. Y. Kitaev, in  _Electr. Coll. Comput. Complex._ (1996) pp. TR96–003.
  * Higgins  _et al._ [2007] B. L. Higgins, D. W. Berry, S. D. Bartlett, H. M. Wiseman, and G. J. Pryde, Nature 450, 393 (2007).
  * Berry  _et al._ [2009] D. W. Berry, B. L. Higgins, S. D. Bartlett, M. W. Mitchell, G. J. Pryde, and H. M. Wiseman, Physical Review A 80, 052114 (2009).
  * Said, Berry, and Twamley [2011] R. S. Said, D. W. Berry, and J. Twamley, Physical Review B 83, 125410 (2011).
  * Cappellaro [2012] P. Cappellaro, Physical Review A 85, 030301 (2012).
  * Bonato  _et al._ [2016] C. Bonato, M. S. Blok, H. T. Dinani, D. W. Berry, M. L. Markham, D. J. Twitchen, and R. Hanson, Nature Nanotechnology 11, 247 (2016).
  * Higgins  _et al._ [2009] B. L. Higgins, D. W. Berry, S. D. Bartlett, M. W. Mitchell, H. M. Wiseman, and G. J. Pryde, New Journal of Physics 11, 073023 (2009).
  * Nusran, Momeen, and Dutt [2012] N. M. Nusran, M. U. Momeen, and M. V. G. Dutt, Nature Nanotechnology 7, 109 (2012).
  * Waldherr  _et al._ [2012] G. Waldherr, J. Beck, P. Neumann, R. S. Said, M. Nitsche, M. L. Markham, D. J. Twitchen, J. Twamley, F. Jelezko, and J. Wrachtrup, Nature Nanotechnology 7, 105 (2012).


#  Entanglement-enabled symmetry-breaking orders

Cheng-Ju Lin  Perimeter Institute for Theoretical Physics, Waterloo, Ontario
N2L 2Y5, Canada  Liujun Zou  Perimeter Institute for Theoretical Physics,
Waterloo, Ontario N2L 2Y5, Canada

###### Abstract

A spontaneous symmetry-breaking order is conventionally described by a tensor-
product wavefunction of some few-body clusters. We discuss a type of symmetry-
breaking orders, dubbed entanglement-enabled symmetry-breaking orders, which
cannot be realized by any tensor-product state. Given a symmetry breaking
pattern, we propose a criterion to diagnose if the symmetry-breaking order is
entanglement-enabled, by examining the compatibility between the symmetries
and the tensor-product description. For concreteness, we present an infinite
family of exactly solvable gapped models on one-dimensional lattices with
nearest-neighbor interactions, whose ground states exhibit entanglement-
enabled symmetry-breaking orders from a discrete symmetry breaking. In
addition, these ground states have gapless edge modes protected by the
unbroken symmetries. We also propose a construction to realize entanglement-
enabled symmetry-breaking orders with spontaneously broken continuous
symmetries. Under the unbroken symmetries, some of our examples can be viewed
as symmetry-protected topological states that are beyond the conventional
classifications.

Introduction. Spontaneous symmetry breaking is a ubiquitous phenomenon and
plays fundamental roles in numerous physical systems, ranging from
astronomical objects like neutron stars to quantum mechanical objects like
atoms forming a crystal. It is also the basis of many important modern
technologies, such as hard drives, maglev trains, and spintronics.

In quantum physics, the standard recipe in describing a spontaneous symmetry
breaking order is by representing it as a tensor-product wavefunction of some
few-body clusters. For example, as depicted in Fig 1, a quantum ferromagnetic
order is often described as a tensor product of spin-ups or spin-downs; a
valence-bond-solid order is described as a tensor product of spin-singlets.

Figure 1: “Classical” symmetry-breaking orders. (a) A quantum ferromagnetic
order. (b) A valence-bond-solid order.

Nevertheless, remarkably, there exist symmetry-breaking orders that cannot be
described by any tensor-product state, fundamentally requiring entanglement.
We dub such an order entanglement-enabled symmetry-breaking order (EESBO),
which will be defined more precisely later. A simple scenario of EESBOs is
where the remaining symmetry has some Lieb-Schultz-Mattis (LSM) constraints
that mandate the states to be long-range entangled. But it is more striking if
the symmetries do not have such constraints but are still incompatible with
any tensor-product wavefunctions. This scenario is the focus of this work.

Here we point out the significance of EESBO. The usual description of the
symmetry-breaking orders suggests that all symmetry-breaking orders are
essentially “classical” Beekman  _et al._ (2019), since they allow a
description by a tensor-product wavefunction. On the other hand, EESBO is
intrinsically quantum, indicating the incompleteness of the above description.
For instance, some enigmatic hidden orders Mydosh and Oppeneer (2011) in some
quantum materials may turn out to be EESBO, evading the usual “classical”
symmetry-breaking order or mean-field-like descriptions. As another example,
recent development in the systematic understanding of the magnon topological
band structures is based on “classical” symmetry-breaking orders Lu and Lu
(2018); Corticelli  _et al._ (2022); Liu  _et al._ (2022); McClarty (2022);
Karaki  _et al._ (2022). However, the complete understanding must incorporate
the magnons arising from EESBOs. Furthermore, the properties of the symmetry
defects associated with EESBOs are largely unexplored and could pave a way to
the future technologies, much like the roles of skyrmions in spintronics Fert
_et al._ (2017); Zhou (2018); Marrows and Zeissler (2021).

Definition of EESBO. To define EESBOs precisely, we need to clarify how to
specify the symmetry setting of a physical system, and the notion of tensor-
product state of few-body clusters. We will see that EESBOs arise from some
incompatibility between the symmetries and tensor-product wavefunctions.

For concreteness, we consider a quantum system on a lattice, described by a
tensor product of local Hilbert spaces. Changing the structures of the local
Hilbert spaces, i.e., their dimensions and operator contents, changes the
physical system. Below we often use “spin-$S$ system” to mean a bosonic system
with $2S+1$ dimensional local Hilbert spaces.

To describe spontaneous symmetry breaking, we first have to specify the
symmetries of the Hamiltonian111We only consider 0-form, ordinary symmetries..
In particular, we need the symmetry setting specified by i) the group $G_{0}$
formed by the symmetry actions on local operators, ii) how states in each
local Hilbert space transform under the symmetry, and iii) the locations of
the degrees of freedom (DOF)222In this definition of a symmetry setting, if
there is a $U(1)$ symmetry, the filling factor is allowed to be tuned after
specifying these three pieces of data, unless it is pinned by other
symmetries.. For example, a spin system on a kagome or honeycomb lattice
described by the usual Heisenberg Hamiltonian has
$G_{0}=SO(3)\times\mathbb{Z}_{2}^{T}\times p6m$, where $SO(3)$,
$\mathbb{Z}_{2}^{T}$ and $p6m$ describe the spin rotation, time reversal and
lattice symmetries, respectively. We stress that the group $G_{0}$ does not
fully specify this symmetry setting without specifying how each spin
transforms under $G_{0}$, e.g., the magnitude of the spin, whether it is a
Kramers doublet, etc. We also need to specify the locations of the spins,
i.e., whether they live on a kagome or honeycomb lattice. However, for
convenience, we just use the group $G_{0}$ to denote the symmetry setting,
keeping other information implicit.

The ground states of the system may spontaneously break the symmetry $G_{0}$
to its subgroup $G$,333How states transform under $G$ is determined by how
they transform under $G_{0}$, and the locations of the DOF are often unchanged
after spontaneous symmetry breaking. So the symmetry setting after spontaneous
symmetry breaking is fully specified. and we denote this symmetry-breaking
order by “$G_{0}\rightarrow G$”. Given $G_{0}$ and $G$, we call the symmetry-
breaking order “classical” if its $G$-symmetric ground states $|\psi\rangle$
can be represented by a local product state, i.e.,
$|\psi\rangle=\bigotimes_{i}|\psi_{i}\rangle_{\Lambda_{i}}$, where
$|\psi_{i}\rangle_{\Lambda_{i}}$ is supported only in a local region
$\Lambda_{i}$ for all $i$. Here by local we mean that the “size” of
$\Lambda_{i}$ does not scale with the system size, and we will call the state
“ultralocal” if each $\Lambda_{i}$ is just one site. Note that although a
“classical” $G_{0}\rightarrow G$ symmetry-breaking order can be realized by
local product states, generically it can also be realized by other states,
which may have (even long-range) entanglement, since $G$ alone does not fully
determine the ground states. If all possible $G$-symmetric ground states
cannot be realized by local product states but can be realized by some
entangled states, we call the $G_{0}\rightarrow G$ symmetry-breaking order an
EESBO.

Mechanisms for EESBOs. What are the mechanisms for EESBOs? A simple scenario
is that there are LSM-type constraints444As commented before, in our
definition the $U(1)$ filling factor is not fixed unless it is pinned by other
symmetries, so there is no LSM constraint purely associated with the filling
factor. that mandate all $G$-symmetric ground states to be long-range
entangled Lieb  _et al._ (1961); Oshikawa (2000); Hastings (2004); Po  _et
al._ (2017), so the corresponding symmetry-breaking order must be
entanglement-enabled. For instance, consider a kagome lattice qubit system
with $G_{0}=SO(3)\times\mathbb{Z}_{2}^{T}\times p6m$, i.e., $G_{0}$ includes
an $SO(3)$ spin rotation symmetry where the qubits carry spin-1/2, a
$\mathbb{Z}_{2}^{T}$ time reversal symmetry where the qubits transform as
Kramers doublets, and a $p6m$ lattice symmetry that moves the locations of the
qubits. It is proposed that $G_{0}\rightarrow G=SO(3)\times p6m$ for some
Hamiltonian Gong  _et al._ (2014) (i.e., $\mathbb{Z}_{2}^{T}$ is spontaneously
broken), which has an LSM constraint mandating all $G$-symmetric ground states
to be long-range entangled.

Although interesting by itself, the entanglement-enabled nature of such
symmetry-breaking orders is expected due to the LSM constraints. It may be
more remarkable that EESBOs can arise even if $G$-symmetric states have no LSM
constraint, i.e., $G$-symmetric states can be short-range entangled but cannot
be local product states. Given the original and remaining symmetries $G_{0}$
and $G$ without the LSM constraints, how do we diagnose if the symmetry-
breaking order is entanglement-enabled? Assuming the $G_{0}\rightarrow G$
symmetry-breaking order is realizable, we can have a useful general criterion
to show its entanglement-enabled nature: $G$ contains a lattice symmetry which
constrains the structure of any local product state, and an on-site symmetry
which is incompatible with this constraint. This symmetry $G$ therefore
forbids the possibility of local-product-state realization, forcing the
symmetry-breaking order to be entanglement-enabled. In practice, the lattice
symmetry will often constrain any local product state to be ultralocal, which
is then incompatible with the on-site symmetry.

We note that models exhibiting short-range EESBOs which can be diagnosed in
this fashion have appeared before, see e.g., Ref. Affleck  _et al._ (1991).
However, the motivations to study such models appear to not be based on
EESBOs, and the obstruction of the local-product-state realization of the
order was not shown explicitly therein. It is therefore desirable to
demonstrate the application of this condition for EESBOs in a few classes of
simpler examples, which we show below.

An infinite family of EESBOs. We first show an infinite family of EESBOs that
can be constructed by “dimers”. Each of this infinite family of EESBOs lives
on a one dimensional lattice, where each site hosts an $n^{2}$ dimensional
local Hilbert space, labeled by an integer $n\geqslant 3$. We represent a
basis of the local Hilbert space at site $j$ by $|\alpha,\beta\rangle_{j}$,
where $\alpha,\beta=0,1,\cdots,n-1$ are defined modulo $n$. We consider the
generalized Pauli operators $\mu^{x,z}_{j}$ and $\nu^{x,z}_{j}$, where

$\displaystyle\begin{split}&\mu_{j}^{z}|\alpha,\beta\rangle_{j}=e^{\frac{2\pi
i}{n}\alpha}|\alpha,\beta\rangle_{j}~{},\
\mu_{j}^{x}|\alpha,\beta\rangle_{j}=|\alpha+1,\beta\rangle_{j}~{},\\\
&\nu_{j}^{z}|\alpha,\beta\rangle_{j}=e^{\frac{2\pi
i}{n}\beta}|\alpha,\beta\rangle_{j}~{},\
\nu_{j}^{x}|\alpha,\beta\rangle_{j}=|\alpha,\beta+1\rangle_{j}~{}.\end{split}$

While it is convenient to view variables $\alpha$ and $\beta$ as fictitious
DOF, it is crucial to keep in mind that the local Hilbert space is formed by
$\alpha$ and $\beta$ together; in particular, one cannot discuss the
entanglement between the $\alpha$ and $\beta$ DOF at the same site.

The symmetry group $G_{0}$ of the Hamiltonian contains a
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$ internal symmetry, where
$\mathbb{Z}_{n}^{x}$ is implemented by $\prod_{j}\mu^{x}_{j}\nu^{x}_{j}$, and
$\mathbb{Z}_{n}^{z}$ is implemented by
$\prod_{j\in\text{even}}(\mu^{z}_{j}\nu^{z}_{j})\prod_{j\in\text{odd}}(\mu^{z}_{j}\nu^{z}_{j})^{-1}$.
$G_{0}$ also contains a lattice translation symmetry $T$:
$\bigotimes_{j}|\alpha_{j}\beta_{j}\rangle_{j}\rightarrow\bigotimes_{j}|\alpha_{j}\beta_{j}\rangle_{j+1}$
and a reflection symmetry with respect to the site $j\\!=\\!0$, $\sigma$:
$\bigotimes_{j}|\alpha_{j}\beta_{j}\rangle_{j}\rightarrow\bigotimes_{j}|\beta_{j}\alpha_{j}\rangle_{-j}$.

We comment on some noteworthy aspects of the symmetry setting. First, while
the symmetry $G_{0}$ can be defined on an infinite 1d lattice, defining it on
a finite lattice with $L$ sites requires a periodic boundary condition and $L$
to be even. Second, all operators transform in a linear representation of the
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$ symmetry. But the states at each
site transform in its projective representations, and pairs of adjacent sites
together form a linear representation Jiang  _et al._ (2021). Third, $T$ and
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$ do not commute, and our setting
has no LSM constraint.

We consider the case where $T$ is broken to $T^{2}$ with other symmetries
intact, namely, $G_{0}\rightarrow G$, where $G$ is generated by
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$, $\sigma$ and $T^{2}$. Below we
show the obstruction in realizing such a symmetry breaking order as a local
product state using the strategy laid out previously. We first show that the
lattice symmetries $T^{2}$ and $\sigma$ enforces the local-product state to be
ultralocal. Assume sites $2j_{1}$ and $2j_{2}$ (or sites $2j_{1}+1$ and
$2j_{2}+1$) are entangled, then $T^{2}$ requires extensively many sites to be
entangled, contradicting the assumption of local product state. If sites
$2j_{1}$ and $2j_{2}+1$ are entangled, then $\sigma$ requires sites $2j_{1}$
and $4j_{1}\\!-\\!2j_{2}\\!-\\!1$ to be entangled. $T^{2}$ further requires
sites $4j_{1}\\!-\\!2j_{2}\\!-\\!1$ and $6j_{1}-\\!4j_{2}\\!-\\!2$ to be
entangled. By iteration, we see that extensively many sites would be
entangled, resulting in a contradiction again. So any local product state must
be ultralocal. However, because the states at each site transform in
projective representations of the $\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$
symmetry (recall $n\geqslant 3$), no ultralocal product state is invariant
under $\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$.

Figure 2: The symmetry-breaking ground states of the Hamiltonian Eq. (1) with
periodic boundary condition. With open boundary condition, the ground states
have dangling edge degrees of freedom, depicted by the same figure but with
the dashed lines removed.

Despite the obstruction of a local-product-state realization, can this
symmetry breaking order be realized at all? Indeed in the following, we
present an explicit Hamiltonian and the ground states, which establish the
existence of EESBOs for this $G_{0}\rightarrow G$. Consider the Hamiltonian on
a 1d lattice with even number of sites $L$ and a periodic boundary condition
($j+L\equiv j$),

$H=\sum_{j=1}^{L}(I-P_{j,j+1})~{},$ (1)

where $I$ is the identity operator, $P_{j,j+1}$ is the projector onto the
linear space
$\text{span}\\{|D^{(\alpha)}_{a,b}\rangle,|D^{(\beta)}_{a,b}\rangle,a,b=1\dots
n\\}$ and

$\displaystyle|D^{(\alpha)}_{a,b}\rangle$ $\displaystyle=$
$\displaystyle\sum_{d=1}^{n}|\alpha_{j}\\!=\\!d,\beta_{j}\\!=\\!a,\alpha_{j+1}\\!=\\!d,\beta_{j+1}\\!=\\!b\rangle$
$\displaystyle|D^{(\beta)}_{a,b}\rangle$ $\displaystyle=$
$\displaystyle\sum_{d=1}^{n}|\alpha_{j}\\!=\\!a,\beta_{j}\\!=\\!d,\alpha_{j+1}\\!=\\!b,\beta_{j+1}\\!=\\!d\rangle$

are the maximum entangled states (“dimers”) on the degrees of freedom $\alpha$
or $\beta$. It is easy to verify that the two (unnormalized) $G$-symmetric
states

$\displaystyle|\psi_{A}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha_{j},\beta_{j}\\}}\prod_{j=1}^{L/2}\delta_{\alpha_{2j}\alpha_{2j+1}}\delta_{\beta_{2j-1}\beta_{2j}}|\alpha_{1}\beta_{1}\dots\alpha_{L}\beta_{L}\rangle~{},$
$\displaystyle|\psi_{B}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha_{j},\beta_{j}\\}}\prod_{j=1}^{L/2}\delta_{\alpha_{2j-1}\alpha_{2j}}\delta_{\beta_{2j}\beta_{2j+1}}|\alpha_{1}\beta_{1}\dots\alpha_{L}\beta_{L}\rangle~{}$

are the ground states, i.e., $H|\psi_{A(B)}\rangle=0$, which are depicted in
Fig. 2. In SM sup , we show that these two states are indeed the only two
ground states of $H$, and that the model has a spectral gap in the
thermodynamic limit. Note that these two states are orthogonal only in the
thermodynamic limit but are linearly independent at any finite $L$.

In fact, these symmetry-breaking ground states are
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$ symmetry-protected topological
(SPT) states. This can be shown by examining the ground states with open
boundary condition. In this case, the dimension of the ground state subspace
becomes $2n^{2}$, where the degeneracy comes from the edge states of the
dangling $\alpha$ or $\beta$ DOF, as depicted in Fig. 2 with the dashed lines
removed. These edge DOF transform projectively under
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$, which is a hallmark of 1d
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}$ SPTs Pollmann  _et al._ (2012,
2010); Chen  _et al._ (2011). More interestingly, under the remaining
$\mathbb{Z}_{n}^{x}\times\mathbb{Z}_{n}^{z}\times\sigma$ symmetry, these
ground states can be viewed as SPTs beyond the conventional group-cohomology-
based classification Jiang and Ran (2017); Thorngren and Else (2018), which is
possible because the states at each site transform projectively under the on-
site symmetries Jiang  _et al._ (2021). To capture such states, a more refined
classification like the ones in Refs. Else and Thorngren (2020a); Jiang  _et
al._ (2021) is needed. We leave the derivation of these statements to SM sup .

Examples with spontaneously broken continuous symmetries. The above discussion
firmly establishes that EESBOs can exist even without LSM constraints
associated with the remaining symmetry $G$. Below we present examples of
entanglement-enabled spontaneously broken continuous symmetries, which result
in richer structures and illustrate some important features of EESBOs.

Consider a $d$-dimensional spin-$S$ system with
$G_{0}=SO(3)\times\mathbb{Z}^{d}$, where the $\mathbb{Z}^{d}$ translation
symmetric lattice has one site per unit cell, and a spin-$S$ moment lives on
each site. Consider the case where $d>1$ and
$G=\mathbb{Z}_{2}\times\mathbb{Z}^{d}$, i.e., $SO(3)$ is broken to
$\mathbb{Z}_{2}$ while the translation is intact555The Mermin-Wagner theorem
forbid such a symmetry-breaking order for $d=1$.. It is easy to verify that
the translation symmetry forces all local product states to be ultralocal. If
$S\geqslant 1$, ultralocal product states can realize such a symmetry-breaking
order without any obstruction sup . On the other hand, if $S=1/2$, all
translation symmetric ultralocal product states have a remaining $U(1)$ on-
site symmetry, so $SO(3)$ can never be broken down to $\mathbb{Z}_{2}$. We
therefore see that for $S=1/2$, the $G_{0}\rightarrow G$ symmetry-breaking
order must be entanglement-enabled. Note that there is no LSM constraint
associated with the $G$ symmetry. A representative wavefunction of this EESBO
can be taken as a stack of spin-$1/2$ chains (see Fig. 3), where each chain is
described by a matrix-product state $|\psi\rangle=\sum_{\\{s_{j}\\}}{\rm
tr}(A^{[s_{1}]}\dots A^{[s_{L}]})|s_{1}\dots s_{L}\rangle$ with matrices

$\displaystyle A^{[s_{z}=\frac{1}{2}]}=\left(\begin{array}[]{cc}1+a&0\\\
0&1-a\end{array}\right),A^{[s_{z}=-\frac{1}{2}]}=\left(\begin{array}[]{cc}0&b\\\
c&0\end{array}\right)$ (6)

Each chain has a translation symmetry along itself, and different chains are
arranged so that the entire system has $\mathbb{Z}^{d}$ translation symmetry.
In SM sup , we show that for generic nonzero real-valued $a,b$ and $c$, the
remaining symmetry of this system is indeed $G$, and the system can be viewed
as a $G$-SPT beyond the conventional group-cohomology-based classification in
Refs. Cheng  _et al._ (2016); Jiang and Ran (2017); Thorngren and Else (2018)
and in some sense even beyond the more refined classification in Refs. Else
and Thorngren (2020a); Jiang  _et al._ (2021), which is possible because the
$\mathbb{Z}_{2}$ symmetry acts on the states at each site as a
$\mathbb{Z}_{4}$ symmetry. Furthermore, in these symmetry-breaking orders, for
all $S\geqslant 1/2$, there are generically two gapless Goldstone modes, where
one of them has a linear dispersion and the other has a quadratic dispersion.

Figure 3: Stacks of spin-1/2 chains that realize the
$G_{0}=SO(3)\times\mathbb{Z}^{d}\rightarrow
G=\mathbb{Z}_{2}\times\mathbb{Z}^{d}$ EESBOs in 2 (left) and 3 (right)
dimensions, where each chain is described by the matrix-product state in Eq.
(6).

This example highlights some intriguing points. First, in the previous
examples, after fixing $G$, $G_{0}$ is immaterial in determining whether the
$G_{0}\\!\rightarrow\\!G$ symmetry-breaking order is entanglement-enabled.
However, the nature of $G_{0}$ is important here. For instance, if the $SO(3)$
spin rotation symmetry is replaced by $O(2)$, $G_{0}=O(2)\times\mathbb{Z}^{d}$
still has $G=\mathbb{Z}_{2}\times\mathbb{Z}^{d}$ as its subgroup, but the
$G_{0}\rightarrow G$ symmetry-breaking order is “classical” for all $S$.
Second, it is interesting to see that fixing $G_{0}$ but enlarging
$G=\mathbb{Z}_{2}\times\mathbb{Z}^{d}$ to $G=U(1)\times\mathbb{Z}^{d}$ will
make the would-be $G_{0}\rightarrow G$ EESBO for $S=1/2$ “classical”. Third,
this example shows that whether a symmetry-breaking order is entanglement-
enabled can depend on the precise representation of DOF under the on-site
symmetry, rather than only the equivalence class of projective
representations.

In SM sup , we present two more examples of EESBOs with spontaneously broken
continuous symmetries (but these by no means exhaust all such EESBOs). One
example is a spin-1 lattice system with $G_{0}=PSU(3)\times\mathbb{Z}^{d}$ and
$G=SO(3)\times\mathbb{Z}^{d}$, which may be realized in some generalizations
of the bilinear-biquadratic model (see, e.g., Ref. Papanicolaou (1988) for a
discussion of the bilinear-biquadratic model). The previous two examples can
be viewed as SPTs under their remaining symmetries that are beyond the
conventional classifications, but this state is a within-group-cohomology SPT
under its remaining symmetry. The other example is a honeycomb lattice
spin-$1/2$ system with $G_{0}=SO(3)\times\mathbb{Z}_{2}^{T}\times p6m$ and
$G=\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}^{T}\times p6m$,
i.e., an entanglement-enabled spin-nematic order. Ref. Ye  _et al._ (2021)
proposed that the vicinity of this order may host a quantum spin liquid that
is beyond the usual parton description.

Discussion. We have discussed various examples of EESBOs, which are
intrinsically quantum symmetry-breaking orders that elude the standard
“classical” descriptions. While one mechanism of EESBO is due to the LSM
constraints of the remaining symmetry, we present a criterion for EESBO in the
cases without LSM constraints. Our work paves the way to discovering more
interesting aspects and insightful results about EESBOs. We finish this paper
by discussing some future directions.

For our infinite family of one dimensional EESBOs, it is interesting to
generalize them to higher dimensions. Since these EESBOs have symmetry-
protected gapless boundaries, it is natural to expect the transitions out of
these phases to be gapless SPT states Keselman and Berg (2015); Scaffidi  _et
al._ (2017); Li  _et al._ (2022), which warrant a further study. For the
EESBOs with spontaneously broken continuous symmetries, finding their parent
Hamiltonians is another important problem.

It might also be rewarding to connect the intrinsic quantumness of EESBOs to
the complexity of (path-integral) quantum Monte Carlo calculations from the
perspective of quantum-classical mapping. The Hamiltonian Eq. (1) is not
stoquastic in either the $\mu^{z}$/$\nu^{z}$-diagonal basis or the
$\mu^{x}$/$\nu^{x}$-diagonal basis, though the ground state wavefunctions are
nonnegative in both basis. What would a corresponding “classical statistical
mechanics model” look like for EESBO, if such a sign problem of its
Hamiltonian is curable? Could the intrinsic quantumness of EESBOs suggest that
it is not curable by some (symmetry-preserving) local basis transformation,
signifying the sign-problem being symmetry protected or even intrinsic Ringel
and Kovrizhin (2017); Smith  _et al._ (2020); Ellison  _et al._ (2021); Lin
_et al._ ?

Although much of the universal low-energy physics of EESBOs can be captured by
the usual effective field theory approach, it may require new methods to
understand other important properties of an EESBO due to the absence of a
“classical” mean-field description. For example, in EESBOs with spontaneously
broken continuous symmetry, it is desirable to understand the band structures
of the Goldstone modes. The conventional method is to analyze the effective
theory of these Goldstone modes based on a “classical” mean-field description,
which is unavailable in EESBOs. To proceed, one may need techniques like the
parton approach and/or tensor networks. It is also interesting to understand
the entanglement-enabled nature of these symmetry-breaking orders from a field
theoretic perspective, similar to the description of LSM constraints using
topological field theories Cheng  _et al._ (2016); Else and Thorngren (2020a);
Ye  _et al._ (2021).

A more formal challenge is to classify all EESBOs. From our examples, we see
that EESBOs often arise when the local DOF carry specific representations
under the symmetry. This is analogous to the physics of obstructed insulators,
fragile topology and generalized LSM constraints studied in the context of
symmetric states Po  _et al._ (2018); Else  _et al._ (2019); Latimer and Wang
(2021); Lu (2017); Else and Thorngren (2020b); Jiang  _et al._ (2021), because
in all these cases certain many-body entanglement pattern is incompatible with
the properties of local DOF. It is interesting to see if this analogy can be
made sharper and useful. For example, one can ask how the EESBO ground states
with spontaneously broken continuous symmetries constrain the possible band
topology of the Goldstone modes. To this end, anomaly-based methods developed
in Refs. Zou  _et al._ (2021); Ye  _et al._ (2021) may be helpful.

Finally, it is important to identify more models and experimental systems
exhibiting EESBOs. It is also worth exploring whether the corresponding
symmetry defects, whose properties are not well understood so far, have
potential technological applications.

###### Acknowledgements.

We thank Yin-Chen He, Timothy Hsieh, Olexei Motrunich, Yang Qi and Chong Wang
for valuable discussions and feedback. Research at Perimeter Institute is
supported in part by the Government of Canada through the Department of
Innovation, Science and Economic Development Canada and by the Province of
Ontario through the Ministry of Colleges and Universities.

## References

  * Beekman  _et al._ (2019) Aron J. Beekman, Louk Rademaker, and Jasper van Wezel, “An Introduction to Spontaneous Symmetry Breaking,” SciPost Phys. Lect. Notes , 11 (2019).
  * Mydosh and Oppeneer (2011) J. A. Mydosh and P. M. Oppeneer, “Colloquium: Hidden order, superconductivity, and magnetism: The unsolved case of URu2Si2,” Rev. Mod. Phys. 83, 1301–1322 (2011), arXiv:1107.0258 [cond-mat.str-el] .
  * Lu and Lu (2018) Fuyan Lu and Yuan-Ming Lu, “Magnon band topology in spin-orbital coupled magnets: classification and application to $\alpha$-RuCl3,” arXiv e-prints , arXiv:1807.05232 (2018), arXiv:1807.05232 [cond-mat.str-el] .
  * Corticelli  _et al._ (2022) A. Corticelli, R. Moessner, and P. A. McClarty, “Spin-space groups and magnon band topology,” Phys. Rev. B 105, 064430 (2022), arXiv:2103.05656 [cond-mat.str-el] .
  * Liu  _et al._ (2022) Pengfei Liu, Jiayu Li, Jingzhi Han, Xiangang Wan, and Qihang Liu, “Spin-Group Symmetry in Magnetic Materials with Negligible Spin-Orbit Coupling,” Phys. Rev. X 12, 021016 (2022), arXiv:2103.15723 [cond-mat.mtrl-sci] .
  * McClarty (2022) Paul A. McClarty, “Topological magnons: A review,” Annual Review of Condensed Matter Physics 13, 171–190 (2022).
  * Karaki  _et al._ (2022) Mohammed J. Karaki, Xu Yang, Archibald J. Williams, Mohamed Nawwar, Vicky Doan-Nguyen, Joshua E. Goldberger, and Yuan-Ming Lu, “An efficient material search for room temperature topological magnons,” arXiv e-prints , arXiv:2206.06248 (2022), arXiv:2206.06248 [cond-mat.mtrl-sci] .
  * Fert  _et al._ (2017) Albert Fert, Nicolas Reyren, and Vincent Cros, “Magnetic skyrmions: advances in physics and potential applications,” Nature Reviews Materials 2, 17031 (2017), arXiv:1712.07236 [cond-mat.mtrl-sci] .
  * Zhou (2018) Yan Zhou, “Magnetic skyrmions: intriguing physics and new spintronic device concepts,” National Science Review 6, 210–212 (2018).
  * Marrows and Zeissler (2021) C. H. Marrows and K. Zeissler, “Perspective on skyrmion spintronics,” Applied Physics Letters 119, 250502 (2021).
  * Lieb  _et al._ (1961) Elliott Lieb, Theodore Schultz, and Daniel Mattis, “Two soluble models of an antiferromagnetic chain,” Annals of Physics 16, 407 – 466 (1961).
  * Oshikawa (2000) Masaki Oshikawa, “Commensurability, Excitation Gap, and Topology in Quantum Many-Particle Systems on a Periodic Lattice,” Phys. Rev. Lett. 84, 1535–1538 (2000), arXiv:cond-mat/9911137 [cond-mat.str-el] .
  * Hastings (2004) M. B. Hastings, “Lieb-Schultz-Mattis in higher dimensions,” Phys. Rev. B 69, 104431 (2004), arXiv:cond-mat/0305505 [cond-mat.str-el] .
  * Po  _et al._ (2017) Hoi Chun Po, Haruki Watanabe, Chao-Ming Jian, and Michael P. Zaletel, “Lattice Homotopy Constraints on Phases of Quantum Magnets,” Phys. Rev. Lett. 119, 127202 (2017), arXiv:1703.06882 [cond-mat.str-el] .
  * Gong  _et al._ (2014) Shou-Shu Gong, Wei Zhu, and D. N. Sheng, “Emergent Chiral Spin Liquid: Fractional Quantum Hall Effect in a Kagome Heisenberg Model,” Scientific Reports 4, 6317 (2014), arXiv:1312.4519 [cond-mat.str-el] .
  * Affleck  _et al._ (1991) Ian Affleck, D.P. Arovas, J.B. Marston, and D.A. Rabson, “Su(2n) quantum antiferromagnets with exact c-breaking ground states,” Nuclear Physics B 366, 467–506 (1991).
  * Jiang  _et al._ (2021) Shenghan Jiang, Meng Cheng, Yang Qi, and Yuan-Ming Lu, “Generalized Lieb-Schultz-Mattis theorem on bosonic symmetry protected topological phases,” SciPost Phys. 11, 24 (2021).
  * (18) Supplemental Material.
  * Pollmann  _et al._ (2012) Frank Pollmann, Erez Berg, Ari M. Turner, and Masaki Oshikawa, “Symmetry protection of topological phases in one-dimensional quantum spin systems,” Phys. Rev. B 85, 075125 (2012), arXiv:0909.4059 [cond-mat.str-el] .
  * Pollmann  _et al._ (2010) Frank Pollmann, Ari M. Turner, Erez Berg, and Masaki Oshikawa, “Entanglement spectrum of a topological phase in one dimension,” Phys. Rev. B 81, 064439 (2010), arXiv:0910.1811 [cond-mat.str-el] .
  * Chen  _et al._ (2011) Xie Chen, Zheng-Cheng Gu, and Xiao-Gang Wen, “Classification of gapped symmetric phases in one-dimensional spin systems,” Phys. Rev. B 83, 035107 (2011), arXiv:1008.3745 [cond-mat.str-el] .
  * Jiang and Ran (2017) Shenghan Jiang and Ying Ran, “Anyon condensation and a generic tensor-network construction for symmetry-protected topological phases,” Phys. Rev. B 95, 125107 (2017), arXiv:1611.07652 [cond-mat.str-el] .
  * Thorngren and Else (2018) Ryan Thorngren and Dominic V. Else, “Gauging spatial symmetries and the classification of topological crystalline phases,” Phys. Rev. X 8, 011040 (2018).
  * Else and Thorngren (2020a) Dominic V. Else and Ryan Thorngren, “Topological theory of Lieb-Schultz-Mattis theorems in quantum spin systems,” Phys. Rev. B 101, 224437 (2020a), arXiv:1907.08204 [cond-mat.str-el] .
  * Cheng  _et al._ (2016) Meng Cheng, Michael Zaletel, Maissam Barkeshli, Ashvin Vishwanath, and Parsa Bonderson, “Translational Symmetry and Microscopic Constraints on Symmetry-Enriched Topological Phases: A View from the Surface,” Physical Review X 6, 041068 (2016), arXiv:1511.02263 [cond-mat.str-el] .
  * Papanicolaou (1988) N. Papanicolaou, “Unusual phases in quantum spin-1 systems,” Nuclear Physics B 305, 367–395 (1988).
  * Ye  _et al._ (2021) Weicheng Ye, Meng Guo, Yin-Chen He, Chong Wang, and Liujun Zou, “Topological characterization of Lieb-Schultz-Mattis constraints and applications to symmetry-enriched quantum criticality,” arXiv e-prints , arXiv:2111.12097 (2021), arXiv:2111.12097 [cond-mat.str-el] .
  * Keselman and Berg (2015) Anna Keselman and Erez Berg, “Gapless symmetry-protected topological phase of fermions in one dimension,” Phys. Rev. B 91, 235309 (2015), arXiv:1502.02037 [cond-mat.mes-hall] .
  * Scaffidi  _et al._ (2017) Thomas Scaffidi, Daniel E. Parker, and Romain Vasseur, “Gapless Symmetry-Protected Topological Order,” Phys. Rev. X 7, 041048 (2017), arXiv:1705.01557 [cond-mat.str-el] .
  * Li  _et al._ (2022) Linhao Li, Masaki Oshikawa, and Yunqin Zheng, “Symmetry Protected Topological Criticality: Decorated Defect Construction, Signatures and Stability,” arXiv e-prints , arXiv:2204.03131 (2022), arXiv:2204.03131 [cond-mat.str-el] .
  * Ringel and Kovrizhin (2017) Zohar Ringel and Dmitry L. Kovrizhin, “Quantized gravitational responses, the sign problem, and quantum complexity,” Sci. Adv. 3, e1701758 (2017).
  * Smith  _et al._ (2020) Adam Smith, Omri Golan, and Zohar Ringel, “Intrinsic sign problems in topological quantum field theories,” Phys. Rev. Research 2, 033515 (2020), arXiv:2005.05343 .
  * Ellison  _et al._ (2021) Tyler D. Ellison, Kohtaro Kato, Zi-Wen Liu, and Timothy H. Hsieh, “Symmetry-protected sign problem and magic in quantum phases of matter,” Quantum 5, 612 (2021).
  * (34) Cheng-Ju Lin, Weicheng Ye, Yijian Zou, Shengqi Sang, and Timothy H. Hsieh, “Probing sign structure using measurement-induced entanglement,” arXiv:2205.05692 .
  * Po  _et al._ (2018) Hoi Chun Po, Haruki Watanabe, and Ashvin Vishwanath, “Fragile Topology and Wannier Obstructions,” Phys. Rev. Lett. 121, 126402 (2018), arXiv:1709.06551 [cond-mat.mes-hall] .
  * Else  _et al._ (2019) Dominic V. Else, Hoi Chun Po, and Haruki Watanabe, “Fragile topological phases in interacting systems,” Phys. Rev. B 99, 125122 (2019), arXiv:1809.02128 [cond-mat.str-el] .
  * Latimer and Wang (2021) Katherine Latimer and Chong Wang, “Correlated fragile topology: A parton approach,” Phys. Rev. B 103, 045128 (2021), arXiv:2007.15605 [cond-mat.str-el] .
  * Lu (2017) Yuan-Ming Lu, “Lieb-Schultz-Mattis theorems for symmetry protected topological phases,” arXiv e-prints , arXiv:1705.04691 (2017), arXiv:1705.04691 [cond-mat.str-el] .
  * Else and Thorngren (2020b) Dominic V. Else and Ryan Thorngren, “Topological theory of Lieb-Schultz-Mattis theorems in quantum spin systems,” Phys. Rev. B 101, 224437 (2020b), arXiv:1907.08204 [cond-mat.str-el] .
  * Zou  _et al._ (2021) Liujun Zou, Yin-Chen He, and Chong Wang, “Stiefel Liquids: Possible Non-Lagrangian Quantum Criticality from Intertwined Orders,” Phys. Rev. X 11, 031043 (2021), arXiv:2101.07805 [cond-mat.str-el] .