††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Order by disorder in frustration-free systems: a quantum Monte Carlo study
of a two-dimensional PXP model

Mingxi Yue  Wilczek Quantum Center, School of Physics and Astronomy, Shanghai
Jiao Tong University, Shanghai 200240, China  Zijian Wang  Wilczek Quantum
Center, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China  Zhiyuan College, Shanghai Jiao Tong University,
Shanghai 200240, China  Bhaskar Mukherjee  Wilczek Quantum Center, School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh,
Pennsylvania 15260, USA  Zi Cai<EMAIL_ADDRESS>Wilczek Quantum Center,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Key Laboratory of Artificial Structures and Quantum Control,
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Shanghai Research Center for Quantum Sciences, Shanghai 201315,
China

###### Abstract

In this study, an order by disorder mechanism has been proposed in a two-
dimensional PXP model, where the extensive degeneracy of the classical ground-
state manifold is due to strict occupation constraints instead of geometrical
frustrations. By performing an unbias large-scale quantum Monte Carlo
simulation, we find that local quantum fluctuations, which usually work
against long-range ordering, lift the macroscopic classical degeneracy and
give rise to a compressible ground state with charge-density-wave long-range
order. A simple trial wavefunction has been proposed to capture the essence of
the ground-state of the two-dimensional PXP model. The finite temperature
properties of this model have also been studied, and we find a thermal phase
transition with an universality class of two-dimensional Ising model.

Introduction – Order by disorder mechanism is one of the most remarkable
phenomena in geometrically frustrated materials, where it is impossible to
simultaneously minimize interacting energy for all bonds in the system. Hence
extensive degeneracy arises in the classical ground-state manifold. Such a
classical ground-state degeneracy may be partially or completely lifted by
subtle effects that normally work against ordering in a process termed “order
by disorder”, where particular configurations are selected out of the
degenerate manifold and long-range order is restoredVillain (1979); Henley
(1997, 1989). For instance, thermal fluctuations can give rises to entropic
differences between configurations, hence degeneracy-breaking free-energy
terms at finite temperatureVillain (1979); Henley (1997), while at zero
temperature, taking into account quantum fluctuations may leads to degenerate-
breaking zero point energy that favors ordered statesHenley (1989); Chubukov
(1992). Up to now, most studies in this field concentrate on the frustrated
systems, a question thus arises: whether frustration is a necessity for us to
observe the order by disorder phenomenon?

In this paper, we attempt to answer this question, focusing for simplicity on
a two-dimensional(2D) PXP model on a square lattice. In spite of its extreme
simplicity, this model does not only provide a prototypical example of order
by disorder phenomena in frustration-free system, but is also closely related
with recent progress with quantum simulators based on Rydberg atomic
systemsEndres et al. (2016); Barredo et al. (2016); Kumar et al. (2018), where
a neutral atom can strongly interact with its neighbors via a Rydberg blockade
mechanism which prevents two nearby atoms from being simultaneously excited
into the excited statesJaksch et al. (2000). In such a highly tunable and
controllable platform, the interplay between constraints resulting from the
Rydberg blockade and the quantum fluctuations could give rise to exotic
quantum mattersLaycock et al. (2011); Lesanovsky (2012); Labuhn et al. (2016);
de Léséleuc et al. (2019); Samajdar et al. (2020) as well as non-equilibrium
quantum many-body dynamicsBernien et al. (2017); Keesling et al. (2019);
Turner et al. (2019); Ho et al. (2019).

By performing a numerically exact quantum Monte Carlo (QMC)
simulationProkof’ev et al. (1998), we study the ground-state of a PXP model in
a square lattice, which can be realized in an interacting Rydberg atomic
system with zero detuningLesanovsky and Katsura (2012). It is shown that in
the absence of laser driving, any classical configuration satisfying the
Rydberg blockade condition is a ground-state of this model, yielding a
macroscopic degeneracy in the ground-state manifold. For each atom in the
lattice, the driving laser field induces coupling between the ground state and
Rydberg state, thus acts as a local quantum fluctuation. Counter-intuitively,
such an off-diagonal term (in the Fock basis), which usually works against
long-range ordering, is found to lift the macroscopic classical degeneracy,
and give rise to a charge-density-wave state with a spontaneous $Z_{2}$
symmetry breaking. This result has provided a different perspective view of
the order by disorder phenomenon, where the extensive degeneracy lifted by
quantum fluctuations comes from occupation constraints instead of geometrical
frustration.

Model and method – In terms of the hard-core boson language, the Hamiltonian
of a PXP model reads

$H=\Omega\sum_{\mathbf{i}}(\tilde{b}_{\mathbf{i}}+\tilde{b}_{\mathbf{i}}^{\dagger})$
(1)

where $\tilde{b}_{\mathbf{i}}$ ($\tilde{b}_{\mathbf{i}}^{\dagger}$) is the
annihilation (creation) operator of the hard-core bosons operating on the
constraint Hilbert space $\mathbb{H}$ that excludes those configurations with
two bosons located on neighboring sites. Similar models have also been
proposed to describe correlated bosonic systems in tiled optical
latticeSachdev et al. (2002); Sengupta et al. (2004); Pielawa et al. (2011).
In terms of spin-$\frac{1}{2}$ operators, Eq.(1) describes a PXP model with
the Hamiltonian
$H_{PXP}=\Omega\sum_{\mathbf{i}}\tilde{\sigma}_{\mathbf{i}}^{x}$. Compare to
the one dimensional version of the related models whose ground state and non-
equilibrium properties have been intensively studiedOlmos et al. (2012a, b);
Turner et al. (2019); Ho et al. (2019); Lin and Motrunich (2019); Mark et al.
(2020), the 2D PXP model is much less studiedMichailidis et al. (2020a, b).

The PXP model can also be realized in strongly interacting Rydberg atomic
system subjected to a coherent laser driving field, whose Hamiltonian reads as
follows:

$H_{R}=\sum_{\mathbf{i}}[\Omega(|g\rangle_{\mathbf{i}}\langle
e|+|e\rangle_{\mathbf{i}}\langle g|)-\delta|e\rangle_{\mathbf{i}}\langle
e|]+V\sum_{\langle\mathbf{ij}\rangle}\hat{n}^{e}_{\mathbf{i}}\hat{n}^{e}_{\mathbf{j}}$
(2)

where $\langle\mathbf{ij}\rangle$ indicates a pair of adjacent lattice sites
in the $L\times L$ square lattice. $|g\rangle_{i}$ ($|e\rangle_{i}$) denote
the internal atomic ground (excited) state of the atom on site $\mathbf{i}$,
and $\hat{n}_{\mathbf{i}}^{e}=|e\rangle_{\mathbf{i}}\langle e|$ is the density
operator of the excited state on site $\mathbf{i}$. $\Omega$ and $\delta$ are
the Rabi frequency and the detuning of the coherent laser driving field, and
$V$ is the strength of the interaction between the Rydberg atoms. Here we only
consider the nearest neighbor interaction, which can be experimentally
realized by properly tuning the Rydberg blockade radius. The ground-state
phase diagram of Eq.(2) has been studied by density matrix renormalization
group (DMRG) method, where different crystalline ground states and
corresponding quantum phase transitions have been foundSamajdar et al. (2020).
Here, to illustrate its relationship with the order by disorder phenomena, we
consider Eq.(2) in the limit $V\rightarrow\infty$ and $\delta\rightarrow 0$,
where Eq.(2) turns to the 2D PXP model.(1).

Throughout this paper, we study the 2D PXP model described by Eq.(1) using a
continuous-time QMC algorithm, which allows us to study both the ground state
and thermal properties up to a significantly large system size. Notice that
local quantum fluctuation in Eq.(1) does not result in sign problem in the
Monte Carlo sampling, nor does the occupation constraint, thus the QMC
simulation is unbias and numerically exact. To sample those (imaginary)time-
space configurations according to their weights in the partition function, we
choose an ergodic set of updates consists of randomly inserting/removing
onsite pair vertices in the configurationsGreitemann and Pollet (2018). The
constraint complicates the simulation by introducing correlations between
adjacent sites. To approach the ground state, in the QMC simulations, we scale
the inverse temperature as $\beta=L$ and take the thermodynamic limit
$L\rightarrow\infty$.

Figure 1: (Color online) Finite size scaling of the CDW order parameter $m$,
density of the bosons $n$ and the compressibility $\kappa$ in the 2D $L\times
L$ PXP model obtained by QMC simulations with $\beta=L$.

Classical limit – We first focus on the case without quantum fluctuations
($\Omega=0$), where any classical configuration in the Fock basis satisfying
the occupation constraint has the same energy $E=0$. In an one-dimensional(1D)
lattice, the number of these states diverges exponentially with the system
size L as $\sim e^{L\ln\alpha}$ with $\alpha=(\sqrt{5}+1)/2$Domb (1960). In a
2D square lattice, the dimension of this constraint Hilbert space $\mathbb{H}$
still diverges exponentially with the number of lattice sites $N=L^{2}$, while
it is difficult to find an analytical expression as simple as that in the 1D
case.

Order by disorder in the ground state – Now we consider the effect of the
local quantum fluctuation in Eq.(1), which breaks the particle number
conservation. We first focus on the ground state properties, the finite size
scaling of the CDW order parameter

$m=\frac{2}{N}\sqrt{\langle[\sum_{\mathbf{i}}(-1)^{i_{x}+i_{y}}\hat{n}_{\mathbf{i}}]^{2}\rangle}$
(3)

as well as the density of the bosons

$n=\frac{1}{N}\sum_{\mathbf{i}}\langle n_{\mathbf{i}}\rangle$ (4)

are plotted in Fig.1, from which we can find a non-vanishing CDW order
parameter $m_{0}=0.403(6)$ at the density $n_{0}=0.221(5)$ in the
thermodynamic limit. These results indicates that the local quantum
fluctuation, which usually works against long-range ordering, select a quantum
state with CDW long-range order from the extensively degenerate classical
ground-state manifold. Notice that the physical quantities of this model
barely depend on the system size $L$, which indicates that a feature of
locality. In another word, the ground state might be close to a product state.
More information of the ground state can be obtained from the compressibility,
which is defined based on the total number fluctuation as:

$\kappa=\frac{\beta}{L^{2}}(\langle\hat{N}^{2}\rangle-\langle\hat{N}\rangle^{2}),$
(5)

where $\hat{N}=\sum_{\mathbf{i}}\hat{n}_{\mathbf{i}}$ is the total particle
number operator of the hard-core bosons. The finite size scaling of $\kappa$
is also plotted in Fig.1, which indicates that the ground state is
compressible. However, due to the explicit breaking of the U(1) symmetry in
Eq.(1), its ground state is not a superfluid.

The order by disorder phenomenon in this frustration-free system and emergent
CDW order can be understood as a consequence of competition between the
quantum fluctuation and occupation constraint. The quantum fluctuations (spin
flip terms $\tilde{\sigma}^{x}_{i}$ in the spin-$\frac{1}{2}$ language) could
further lower the energy of the system by resonating the degenerate classical
configurations, thus are energetically favored in the ground state. However,
two “spin-flip” operators acting on a pair of adjacent sites will inevitably
give rise to the configurations that violate the occupation constraint. As a
compromise, the system favors a ground state with “spin-flips” as many as
possible, while they try to keep away from each other to avoid Rydberg
blockade. One of the simplest example satisfying these requirements is a
product state with a wavefunction: $\Phi_{A}=\prod_{i}|\phi\rangle_{i}$, where
$|\phi\rangle_{i}=|0\rangle_{i}$ if $i\in A$, and
$|\phi\rangle_{i}=\frac{1}{\sqrt{2}}(|0\rangle_{i}-|1\rangle_{i})$ if $i\in B$
with A(B) indicating the A(B) sublattice in the square lattice.
$|0\rangle_{i}$ ($|1\rangle_{i}$) indicates the ground (excited) state of the
Rydberg atoms on site i. The state $\Phi_{A}$ as well as its sublattice
symmetric counterpart $\Phi_{B}=\Phi_{A\rightarrow B}$ spontaneously break the
translational symmetry, thus exhibits a long rang order which can
qualitatively explain the observed CDW order in the ground state of Eq.(1).

Since there is no spontaneous symmetry breaking for any finite system, we
propose a trial wavefunction

$\Phi_{trial}=\frac{1}{\sqrt{2}}(\Phi_{A}+\Phi_{B})$ (6)

which is an equally-weighted superposition of the symmetry breaking states to
restore the symmetry but keep the long-range CDW correlations. In the Fock
basis, the wavefunction $\Phi_{A}$ ( $\Phi_{B}$) indicate an equally-weighted
superposition of all the possible classical configurations in sublattice A(B),
while leave the other sublattice empty. A similar trial wavefunction without
translational symmetry breaking has been proposed for the 1D PXP
modelOvchinnikov et al. (2003).

Figure 2: (Color online) Schematic diagram of the trial wavefunction in Eq.(6) in terms of spin language.  |  $n$ |  m |  $e_{g}$
---|---|---|---
trial |  0.25 |  $0.5$ |  $-0.5\Omega$
Lanczos($6\times 6$) |  0.209 |  0.417 |  $-0.513\Omega$
QMC ($L\rightarrow\infty$) |  0.221 |  0.403 |  $-0.512\Omega$
Table 1: Comparisons of physical quantities (particle density n, CDW order
parameter m, and average energy $e_{g}=E_{g}/N$) of the ground state obtained
by the trial wavefunction, Lanczos method for a $6\times 6$ system and the QMC
simulations.

Such a trial wavefunction in Eq.(6), in spite of its extreme simplicity, can
qualitatively capture most features of the ground state of Eq.(1). To verify
this point numerically, we compare the results predicted by this trial
wavefunction to our numerical results. We first focus on the physical
quantities, for instance the average energy predicted by the trial
wavefunction .(6) is $-0.5\Omega$, which is only $2.4\%$ higher than its exact
value. In addition, the particle density predicted by wavefunction (6) is
$0.25$, which is also not very far from its exact value $0.22$. To directly
compare the trial wavefunction to the exact one, we perform the exact
diagonalization to obtain the exact ground-state wavefuction
$|\Phi_{E}\rangle$. Notice that the occupation constraint significantly
reduces the dimension of the Hamiltonian matrix, thus enables us to obtain
$|\Phi_{E}\rangle$ by Lanczos algorithm for a small system. Take a $6\times 6$
system for an example, the dimension of the constraint Hilbert space
$\mathbb{H}$ is only $0.0035\%$ of that of the full Hilbert space, and the
overlap between the the exact wavefunction and the trial wavefunction.(6)
$\langle\Phi_{trial}|\Phi_{E}\rangle=0.84888$, which is pretty high
considering the facts that the dimension of these two wavefunctions are
2406862, and there are no any tunable parameter in the trial wavefunction.(6)

Figure 3: (Color online) (a) CDW order parameter m and (b) renormalized
correlation length $\xi/L$ as a function of temperature T for different system
size $L$ in the QMC simulation. The insets are the data collapse for (a) the
order parameter and (b) the renormalized correlation length using the critical
exponents $\nu=1$ and $\beta=1/8$.

Finite temperature phase transition – It has been shown that the ground state
of the 2D PXP Eq.(1) spontaneously breaks the $Z_{2}$ translational symmetry
and exhibit the CDW long range order. Considering the fact that there is only
one parameters $\Omega$ in PXP Eq.(1), at zero temperature, this CDW order
will persist irrespective of the strength of the quantum fluctuation $\Omega$.
At finite temperature, one may expect a competition between quantum and
thermal fluctuations, which may result in a finite-temperature phase
transition. To verify this point numerically, we calculate CDW order parameter
as a function of temperature, as shown in Fig.3(a). We also calculate the CDW
correlation length $\xi$ (along x-direction) in the QMC simulation asSandvik
(2010):

$\xi=\frac{1}{q_{1}}\sqrt{\frac{S(\pi,\pi)}{S(\pi+q_{1},\pi)}-1}$ (7)

where
$S(\mathbf{Q})=\frac{1}{N^{2}}\sum_{ij}e^{i\mathbf{Q}\cdot(\mathbf{i}-\mathbf{j})}\langle
n_{\mathbf{i}}n_{\mathbf{j}}\rangle$ is the structure factor of the density-
density correlation. $q_{1}=\frac{2\pi}{L}$. The CDW order parameter $m$ is
defined the same as Eq.(3). The normalized correlation length $\xi/L$ as a
function of temperature for systems with different sizes $L$ has been plotted
in Fig.3 (b), from which we can find a scaling invariant point indicating a
critical point at $T_{c}=0.4\Omega$. Near the critical point, the correlation
length diverges as $\xi\sim|T-T_{c}|^{-\nu}$, and the order parameter scales
as $m\sim|T-T_{c}|^{\beta}$. To explore the universal class of this phase
transition, we study the scaling behavior around the critical point. The data
collapses as shown in the insets of Fig.3 indicate that the critical exponents
$\beta=1/8$ and $\nu=1$, which agrees with those of the 2D Ising universal
class.

Conclusion and outlook – In summary, we observe a order by disorder phenomenon
in a frustration-free system, where the extensive classical ground state
degeneracy is due to occupation constraint instead of geometrical frustration.
It is shown that in a 2D PXP model, the competition between local quantum
fluctuation and occupation constraint gives rise to a compressible ground
state with CDW long-range order, which could be destroyed by thermal
fluctuation via a finite-temperature phase transition with an universality
class of 2D Ising model.

Future developments will include a generalization of our method and result to
systems with different types of CDW orders and quantum fluctuations. For
instance, Samajdar et al have proposed complex CDW orders in 2D interacting
Rydberg atomic systemsSamajdar et al. (2020). Even though most of them are due
to the direct interactions between the Rydberg atoms in the classical limit,
there could be some regimes in the global phase diagram whose CDW orders are
induced by the order-by-disorder mechanism. Different types of quantum
fluctuations, especially those preserving the total particle
number($\tilde{b}_{i}^{\dagger}\tilde{b}_{j}$ for an instance), are of
particular interest, and one may expect a CDW-to-superfluid transition with
increasing doping in this case. In addition, the groundstate phase diagram of
the PXP models with additional terms (e.g. the chemical potential
$-\mu\tilde{n}_{i}$) and the corresponding quantum critical behavior can also
been studiedSamajdar et al. (2020), even though this term with $\mu\neq 0$
lifts the extensively degeneracy of the corresponding classical model
($\Omega=0$) thus is not directly related with the order-by-disorder phenomena
discussed here. Furthermore, our QMC algorithm is ready to be generalized to
frustrated systems (e.g. kagome and triangle lattices), where the interplay
between the frustration, quantum fluctuations and constraint may give rise to
exotic quantum many-body statesSamajdar et al. (2021); Verresen et al. (2020).

Acknowledgments.—This work is supported by the National Key Research and
Development Program of China (Grant No. 2020YFA0309000 and No.
2016YFA0302001), NSFC of China (Grant No.11574200), Shanghai Municipal Science
and Technology Major Project (Grant No.2019SHZDZX01).

## References

  * Villain (1979) J. Villain, Z. Phys. B 33, 31 (1979).
  * Henley (1997) C. L. Henley, J. Appl. Phys. 61, 3962 (1997).
  * Henley (1989) C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Chubukov (1992) A. Chubukov, Phys. Rev. Lett. 69, 832 (1992).
  * Endres et al. (2016) M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Science 354, 1024 (2016).
  * Barredo et al. (2016) D. Barredo, S. D. Leseleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
  * Kumar et al. (2018) A. Kumar, T.-Y. Wu, F. Giraldo, and D. S. Weiss, Nature 561, 83 (2018).
  * Jaksch et al. (2000) D. Jaksch, J. I. Cirac, P. Zoller, S. L. Rolston, R. Côté, and M. D. Lukin, Phys. Rev. Lett. 85, 2208 (2000).
  * Laycock et al. (2011) T. Laycock, B. Olmos, and I. Lesanovsky, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 184017 (2011).
  * Lesanovsky (2012) I. Lesanovsky, Phys. Rev. Lett. 108, 105301 (2012).
  * Labuhn et al. (2016) H. Labuhn, D. Barredo, S. Ravets, S. D. Leseleuc, T. Macri, T. Lahaye, and A. Browaeys, Nature 534, 667 (2016).
  * de Léséleuc et al. (2019) S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Science 365, 775 (2019), ISSN 0036-8075.
  * Samajdar et al. (2020) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Phys. Rev. Lett. 124, 103601 (2020).
  * Bernien et al. (2017) H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, et al., Nature 551, 579 (2017).
  * Keesling et al. (2019) A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, S. Sachdev, et al., Nature 568, 207 (2019).
  * Turner et al. (2019) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papic, Nature Physics 14, 745 (2019).
  * Ho et al. (2019) W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Phys. Rev. Lett. 122, 040603 (2019).
  * Prokof’ev et al. (1998) N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, Phys. Lett. A 238, 253 (1998).
  * Lesanovsky and Katsura (2012) I. Lesanovsky and H. Katsura, Phys. Rev. A 86, 041601 (2012).
  * Sachdev et al. (2002) S. Sachdev, K. Sengupta, and S. M. Girvin, Phys. Rev. B 66, 075128 (2002).
  * Sengupta et al. (2004) K. Sengupta, S. Powell, and S. Sachdev, Phys. Rev. A 69, 053616 (2004).
  * Pielawa et al. (2011) S. Pielawa, T. Kitagawa, E. Berg, and S. Sachdev, Phys. Rev. B 83, 205135 (2011).
  * Olmos et al. (2012a) B. Olmos, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. Lett. 109, 020403 (2012a).
  * Olmos et al. (2012b) B. Olmos, R. González-Férez, I. Lesanovsky, and L. Velázquez, Journal of Physics A: Mathematical and Theoretical 45, 325301 (2012b).
  * Lin and Motrunich (2019) C.-J. Lin and O. I. Motrunich, Phys. Rev. Lett. 122, 173401 (2019).
  * Mark et al. (2020) D. K. Mark, C.-J. Lin, and O. I. Motrunich, Phys. Rev. B 101, 094308 (2020).
  * Michailidis et al. (2020a) A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Phys. Rev. Research 2, 022065 (2020a).
  * Michailidis et al. (2020b) A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Phys. Rev. X 10, 011055 (2020b).
  * Greitemann and Pollet (2018) J. Greitemann and L. Pollet, SciPost Phys. Lect. Notes p. 2 (2018).
  * Domb (1960) C. Domb, Advances in Physics 9, 194 (1960).
  * Ovchinnikov et al. (2003) A. A. Ovchinnikov, D. V. Dmitriev, V. Y. Krivnov, and V. O. Cheranovskii, Phys. Rev. B 68, 214406 (2003).
  * Sandvik (2010) A. W. Sandvik, AIP Conference Proceedings 1297, 135 (2010).
  * Samajdar et al. (2021) R. Samajdar, W. W. Ho, H. Pichler, M. D. Lukin, and S. Sachdev, Proceedings of the National Academy of Sciences 118 (2021), ISSN 0027-8424.
  * Verresen et al. (2020) R. Verresen, M. D. Lukin, and A. Vishwanath, arXiv e-prints arXiv:2011.12310 (2020), eprint 2011.12310.