# Functional or imperative? On pleasant semantics for differentiable
programming languages

Michael Innes

###### Abstract

In machine learning (ML), researchers and engineers seem to be at odds. System
implementers would prefer models to be declarative, with detailed type
information and semantic restrictions that allow models to be optimised,
rearranged and parallelised. Yet practitioners show an overwhelming preference
for dynamic, imperative languages with mutable state, and much engineering
effort is spent bridging the resulting semantic divide. Is there a fundamental
conflict? This article explores why imperative and functional styles are used,
and how future language designs might get the best of both worlds.

Machine Learning, ICML

##  1 Why functional?

General purpose, numerical programming languages are the right tools for
building models in machine learning, probabilistic programming, and
differentiable programming ($\partial P$). In general, models may be no less
complex than general-purpose programs, and need similar tools to manage that
complexity (Innes et al., 2018a).

ML practitioners tend to use imperative languages – mainly Python, and to a
lesser extent R, Matlab, Julia and a smattering of others – for their work.
Yet ML, and particularly deep learning, frameworks often look quite unlike
normal libraries, imposing an atypical functional style of programming (Colah,
2015).

The term “functional” originated with the use of higher-order functions, and
is now associated with a number of related ideas. Two that matter for our
purposes are:

  1. 1.

Value semantics, particularly for arrays. Mathematical operations create new
arrays, rather than operating in-place.

  2. 2.

Referential transparency, in which a function call can be replaced by its
return value. In other words, implicit dependencies (such as model
configuration objects) can be assumed fixed.

JAX, a Python library for numerics and code transforms (Frostig et al., 2018),
is a clean example of these issues at play. It turns the Python interpreter
into a compiler, using objects that abstractly represent arrays to run a
partial evaluation of user code. In effect the result of dictionary lookups
and dynamic dispatch is cached on first run, which avoids overhead during
model training. But strictly speaking this is in violation of Python’s
semantics: the result of a dictionary lookup, or almost any other Python
operation, might well change at run time. So JAX simply asks users to provide
a pure function and assumes that they have done so, with undefined behaviour
if the rule is violated. This has some limitations and sharp edges that would
be ameliorated with better language-level support.

Julia, a numerical programming language (Bezanson et al., 2017), is revealing
in a different way. Julia’s abstract interpretation is analogous to JAX’s
partial evaluation (a proof-of-concept hybrid, Innes 2020a, showed that they
are two sides of the same coin). But Julia’s compiler must respect the
language’s semantics, without making additional assumptions. Mutable objects,
such as Dicts, stymie inference and must be given explicit type parameters to
allow inference of downstream code. Limited constant propagation is supported
on immutable objects, like Tuples, but not arrays or dictionaries. Julia’s
source-transform based autodiff library, Zygote, disables array mutation in
user functions due to the complexity of managing their derivatives. Deep
learning with Flux, a Julia-based ML stack (Innes et al., 2018b), therefore
looks functional: but in this case the core language’s use of mutation
impedes, rather than enabling, high performance code.

The trend towards functional-style, relatively declarative programming with
compilers like XLA (Sabne, 2020) is a significant shift from the conventional
style of scientific and numerical computing,111And indeed the style Julia was
designed to cater for, albeit with more modern language affordances. in which
arrays are explicitly allocated and modified in loops, with a close
correspondence between program and machine instructions. The change is partly
driven by increasingly diverse hardware, as well as the relatively fixed set
of core numerical kernels used in ML. Kernels individually benefit from hand
optimisation, but compositions of them have less room for improvement, even
for hardware experts. Low-level control is therefore less valuable for users,
and general optimisations like simple loop fusion, memory placement and model-
level distribution and parallelism are best done by compilers, on
representations far more abstract than machine code (Li et al., 2020).

##  2 Why imperative?

The discussion above seems to suggest that Haskell, OCaml or Scala would be
better frontends for numerical compilers. Yet researchers have voted with
their feet (or at least with their keyboards): Python is by far the most
popular language for ML,222Google Trends supports “python machine learning”
being about 10$\times$ more popular than “matlab machine learning”, which in
turn is about 5$\times$ more searched than “julia machine learning”. despite
being poorly suited to the task on paper. Why are dynamic, imperative
languages so popular? Partly this is to do with conventional (though not
inherent) differences between static and dynamic languages – things like
interactivity, ease of use for beginners and the flexibility of dictionary-
oriented programming (Dubois et al., 1996). We don’t delve into those topics,
except to point out that there are perfectly serviceable dynamic, functional
languages like Clojure and Mathematica, so it can’t be a complete
explanation.333There are also some historical accidents that favour Python
specifically. For example, its unusual choice to use reference counting
happens to be ideal for managing large GPU arrays, where Julia’s tracing GC
struggles.

Given that ML frameworks eschew shared mutable state, the last difference of
note between functional and imperative is their conventional syntax choices.
Clearly, nested, structured control flow (as in the loop, if/else statements
and break, return and continue) is a helpful way to express algorithms,
especially the iterative algorithms common in numerical programming. Just as
importantly, the imperative style allows for a notion of change that is easier
to think with, compared with having to decompose and reconstruct everything:
“change column $i$ of matrix $A$” is a more natural description than “take all
columns but the $i^{th}$ from $A$, and concatenate them together with this new
one in the right position to create $B$”, and more so as updates get more
complex.

Syntax like A[i] = x typically does three things: it provides a syntactic
notion of change as above; it efficiently modifies memory in-place; and it may
alter the behaviour of seventeen other threads of execution that reference
$A$. But the last of these is most often an anti-feature, and the second is
the job of the compiler. It really is only the syntax that we want.

Despite being superficial in a sense, syntax is an important part of the user
experience, and a perfectly good reason for users to prefer imperative over
functional languages. But now we are surely asking to have our cake and eat it
too: we don’t like mutable state, yet we want to write programs with explicit
change. In fact, there is no contradiction. To explain why we’ll have to
challenge some misconceptions about “imperative” syntax.

##  3 The procedural is pure

When discussing imperative programming, it’s common to conflate mutable local
variables with mutable data structures. In fact they are hardly related.
Here’s an example with the former, but not the latter:

fn pow(x, n) { r = 1 while n > 0 { r *= x n -= 1 } return r}

LLVM, a compiler for those most imperative of languages (Lattner & Adve,
2004), lowers such functions to an intermediate representation (IR) similar to
the following:

fn pow(x, n) { fn block1() { r = 1 block2(n, r) } fn block2(n, r) { if n > 0 {
block3(n, r) else { block4(r) } } fn block3(n, r) { n_ = n - 1 r_ = r * x
block2(n_, r_) } fn block4(r) { return r } block1()}

LLVM and similar compilers call this form “static single assignment”, or SSA,
because variables are syntactically immutable – in other words, they behave
like bindings in functional languages. We could equally call this a set of
mutually tail-recursive pure functions (Appel, 1998). This form is excellent
for a wide variety of analysis, optimisation and transformation (eg autodiff
on SSA, Innes 2020b). All local control flow, however complex, can be lowered
to SSA (Cytron et al., 1989).

The main point is that “mutable” locals are a corollary of control flow
structures, rather than of other kinds of mutation. And control flow is just
syntax sugar for simple functions. If syntax sugar does not fundamentally
change our programming paradigm, and the latter example is a functional
program, the former must be too. At the least, we can agree that mutable
locals need not violate referential transparency, because there’s no external
way to tell they are being used.444Shared mutable state can also be emulated
in a functional language. But a more invasive, global transformation would be
needed to capture all of its “features”, like race conditions and
nondeterminism, so it represents a more fundamental difference.

It is right to say that variables like $r$ are changing. But we use the word
change in the sense of changing clothes, without assuming any physical
alteration.555change derives from words meaning “barter” or “substitute”;
“cause to turn or pass from one state to another”. A change might compile down
to a real modification (of registers, memory or whatever) but this is only
coincidence. Change and mutation are distinct; C just happens to conflate
them.

##  4 The more things change…

The exploration above suggests that clearer notions of change may help us get
the benefits of imperative-style programming, without losing those of the
functional style. But simply adding control flow and locals to a functional
language is not enough. For one thing, it is clunky to write

xs = append(xs, x)

every time we want to change $xs$ (in this case, by appending a new element to
the list). It has the benefit of being explicit about what is changing, but at
the cost of some verbosity, especially when multiple inputs are affected.

We propose a design that takes inspiration from the C-like address operator &,
writing this as

append(&xs, x)

which is still explicit, but less verbose. We refer to & as the “swap
operator”. All function calls have a hidden return value, a list of swapped
arguments, so that for example

result = foo(&a, b, &c)

is roughly equivalent to

[result, rs] = foo(a, b, c)a = rs[1]c = rs[3]

We can write a function to swap two variables as follows:

fn switch(&x, &y) { [x, y] = [y, x] return}a = 1b = 2switch(&a, &b) # now a =
2, b = 1

where the &a in the function signature indicates that the final value of a
should be returned to the caller.

The use of & is like the address operator in C, C++, Rust and Swift (Ritchie
et al., 1988). But whereas the change implied is similar, those languages
conflate that change with mutation of the stack, with all the knock-on effects
that implies. In our implementation no references, pointers or mutable state
are created; the & is syntax sugar for a functional program, which means that
program transformations and analyses are not impeded.

A last sprinkle help us recover full convenience. Syntax like A[i] = x and
foo.bar = baz can lower to function calls (eg set(&A, x, i)). We can update
nested structures like append(&foo.xs, x), which roughly becomes:

xs = foo.xsappend(&xs, x)foo.xs = xs

This looks imperative. But only variables are changing – not values.

Note that in most languages with mutable locals, closures capture variables by
reference. In this more functional setting they should of course capture by
value instead. This is more intuitive in many cases anyway (eg a closure that
captures a loop counter by reference can lead to confusing bugs).

A side effect of this approach is to avoid a leaky abstraction most imperative
language have: types that happen to be implemented as heap-allocated
references, like dictionaries or objects, are mutable, whereas basic types
like numbers or strings are not. More generally, possible modifications are
dictated by the underlying representation in memory, and the choice to make an
interface in-place or not may be dictated by implementation concerns. When
change and mutation are distinct, any value is just as changeable as any
other.

With these relatively minor tweaks – structured control flow, mutable locals
and a bit of syntax sugar – we can write numerical programs almost exactly as
we would in Python or Julia. But because our data structures are values, we
have the guarantees needed for aggressive compilation and transformation by
default. Such optimisation can be done robustly by the main compiler, without
special cases or sharp edges, and while playing nicely with existing language
tooling: error messages, debuggers and profilers.

##  5 Turning change into mutation

One reason that change, in the sense we have discussed, and mutation are
conflated is that doing things in-place is a good optimisation. Where
performance matters, it is crucial that this optimisation be predictable to
the programmer. If it isn’t, algorithms may not only have unnecessary
overheads but also higher computational complexity (because mutations are
usually constant time, whereas copies are linear). Having to make copies is
therefore a significant downside of the functional style.

It is not fatal, however. Reference counting (and its compile-time analogues,
as in Rust or Lobster, Matsakis & Klock 2014; van Oortmerssen ) provides a
solution: any time the count is $1$, no-one else can see our mutations, so it
is safe to update the old data in place.

Koka (Leijen, 2016) has an implementation of this idea, and it turns out to
benefit from the functional pattern of deconstructing and rebuilding data. For
example, say we have a type Foo with named fields xs, ys and zs. We want to
append a new item, x, to the xs array. In a functional style we might write
the update something like this, using pattern matching to destructure foo:

Foo(xs, ys, zs) = fooxs2 = append(xs, x)foo2 = Foo(xs2, ys, zs)

After the first line has executed, foo is not used, and can therefore be
released. Assuming foo is unique, releasing it will in turn release its hold
on xs. xs is now only referenced from one place (the current stack frame) and
can safely be updated by mutation.

Consider in contrast the version we wrote above, using an imperative-style
update syntax to change foo:

xs = foo.xsappend(&xs, x)foo.xs = xs

Unlike the last example, foo is live after the first line, because we need it
to run the last line. That in turns means xs is referenced in two places (the
current stack frame, and foo) and can’t be modified in place.

Fortunately, the swap syntax we introduced for nested updates, append(&foo.xs,
x), can also account for this case. Rather than expanding directly to the
above, it expands to code that decomposes and rebuilds foo more like the
functional-style example. This can work recursively, meaning that updates to
deeply-nested structures will still reliably happen in place where possible.

##  6 Related work

There are many interesting efforts to bridge the gap between functional and
imperative, including within numerical and differentiable programming.

Dex (Maclaurin et al., 2019) aims to do numerical programming and code
transformation in a Haskell-like setting, notably using constructs for
structured mutation inspired by algebraic effects.

Koka (Leijen, 2016) is a recent functional language with support for side
effects based on algebraic effect handlers. Notably, effects imply that order
matters, so Koka uses a C-like syntax for sequences of statements.

Clojure (Hickey, 2008) is related for somewhat tangential reasons. A key
benefit of languages like Python is the ease of “dictionary-oriented
programming”, where the built-in data types are flexible enough that users
typically don’t have to build their own. Clojure showed that this was
practical in a functional setting with its persistent data structures, which
will be an important component of any system resembling a “functional Python”.

Rust (Matsakis & Klock, 2014) takes a different tack to eliminating shared
mutable state, not by avoiding mutation but by restricting sharing through
static analysis. This approach is probably inappropriate for our use case, but
Rust libraries have thought carefully about how to express more advanced
techniques, such as multiple threads working on an output array, in a
structured and safe way. Those APIs could inspire functional approaches that
can reliably compile to efficient code.

##  7 Conclusion

The functional and imperative styles are less at odds than they may first
appear. At least in machine learning, the functional approach is valuable for
its semantic benefits, while imperative style brings syntactic and expressive
convenience. And though the benefits of FP are real, there is no need to make
users do SSA transformation by hand. Perhaps surprisingly, it is possible to
construct a very imperative looking language that nonetheless retains the most
important functional properties. Such designs are promising for probabilistic
and differentiable programming, offering high flexibility to users expressing
models, without sacrificing the ability to do advanced analysis, optimisation
and transformation of the resulting code.

Thanks to Zenna Tavares for encouraging this paper, and to Maximilian Schleich
for reviewing a draft.

## References

  * Appel (1998) Appel, A. W.  SSA is functional programming.  _ACM Sigplan Notices_ , 33(4):17–20, 1998.
  * Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.  Julia: A fresh approach to numerical computing.  _SIAM review_ , 59(1):65–98, 2017.
  * Colah (2015) Colah, C.  Neural networks, types, and functional programming.  https://colah.github.io/posts/2015-09-NN-Types-FP/, 2015.  Accessed: 2022-05-24.
  * Cytron et al. (1989) Cytron, R., Ferrante, J., Rosen, B. K., Wegman, M. N., and Zadeck, F. K.  An efficient method of computing static single assignment form.  In _Proceedings of the 16th ACM SIGPLAN-SIGACT symposium on Principles of programming languages_ , pp. 25–35, 1989.
  * Dubois et al. (1996) Dubois, P. F., Hinsen, K., and Hugunin, J.  Numerical python.  _Computers in Physics_ , 10(3):262–267, 1996\.
  * Frostig et al. (2018) Frostig, R., Johnson, M. J., and Leary, C.  Compiling machine learning programs via high-level tracing.  _Systems for Machine Learning_ , pp. 23–24, 2018.
  * Hickey (2008) Hickey, R.  The Clojure programming language.  In _Proceedings of the 2008 symposium on Dynamic languages_ , pp. 1–1, 2008.
  * Innes (2020a) Innes, M.  Partial evaluation & abstract interpretation.  https://mikeinnes.io/2020/07/29/mjolnir, 2020a.  Accessed: 2022-05-24.
  * Innes (2020b) Innes, M.  Sense & sensitivities: The path to general-purpose algorithmic differentiation.  _Proceedings of Machine Learning and Systems_ , 2:58–69, 2020b.
  * Innes et al. (2018a) Innes, M., Karpinski, S., Shah, V., Barber, D., Saito Stenetorp, P., Besard, T., Bradbury, J., Churavy, V., Danisch, S., Edelman, A., et al.  On machine learning and programming languages.  Association for Computing Machinery (ACM), 2018a.
  * Innes et al. (2018b) Innes, M., Saba, E., Fischer, K., Gandhi, D., Rudilosso, M. C., Joy, N. M., Karmali, T., Pal, A., and Shah, V.  Fashionable modelling with Flux.  _arXiv preprint arXiv:1811.01457_ , 2018b.
  * Lattner & Adve (2004) Lattner, C. and Adve, V.  Llvm: A compilation framework for lifelong program analysis & transformation.  In _International Symposium on Code Generation and Optimization, 2004\. CGO 2004._ , pp. 75–86. IEEE, 2004.
  * Leijen (2016) Leijen, D.  Algebraic effects for functional programming.  Technical report, Technical Report. MSR-TR-2016-29. Microsoft Research technical report, 2016.
  * Li et al. (2020) Li, M., Liu, Y., Liu, X., Sun, Q., You, X., Yang, H., Luan, Z., Gan, L., Yang, G., and Qian, D.  The deep learning compiler: A comprehensive survey.  _IEEE Transactions on Parallel and Distributed Systems_ , 32(3):708–727, 2020.
  * Maclaurin et al. (2019) Maclaurin, D., Radul, A., Johnson, M. J., and Vytiniotis, D.  Dex: array programming with typed indices.  2019\.
  * Matsakis & Klock (2014) Matsakis, N. D. and Klock, F. S.  The rust language.  _ACM SIGAda Ada Letters_ , 34(3):103–104, 2014\.
  * Ritchie et al. (1988) Ritchie, D. M., Kernighan, B. W., and Lesk, M. E.  _The C programming language_.  Prentice Hall Englewood Cliffs, 1988.
  * Sabne (2020) Sabne, A.  XLA: Compiling machine learning for peak performance.  2020\.
  * (19) van Oortmerssen, W.  Memory management in Lobster.  https://aardappel.github.io/lobster/memory_management.html.  Accessed: 2022-05-24.


# Ergodicity breaking with long range cavity induced quasiperiodic
interactions

Piotr Kubala<EMAIL_ADDRESS>Institute of Theoretical
Physics, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland
Piotr Sierant  The Abdus Salam International Center for Theoretical Physics,
Strada Costiera 11, 34151, Trieste, Italy  Institute of Theoretical Physics,
Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland  Giovanna
Morigi  Theoretical Physics, Saarland University, Campus E2.6, D–66123
Saarbrücken, Germany  Jakub Zakrzewski  Institute of Theoretical Physics,
Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland  Mark Kac
Complex Systems Research Center, Jagiellonian University, Łojasiewicza 11,
30-348 Kraków, Poland<EMAIL_ADDRESS>
###### Abstract

Many-body localization (MBL) behavior is analyzed in an extended Bose-Hubbard
model with quasiperiodic infinite-range interactions. No additional disorder
is present. Examining level statistics and entanglement entropy of eigenstates
we show that a significant fraction of eigenstates of the system is localized
in the presence of strong interactions. In spite of this, our results suggest
that the system becomes ergodic in the standard thermodynamic limit in which
the energy of the system is extensive. At the same time, the MBL regime seems
to be stable if one allows for a super-extensive scaling of the energy. We
show that our findings can be experimentally verified by studies of time
dynamics in many-body cavity quantum electrodynamics setups. The “quench
spectroscopy” is a particularly effective tool that allows us to
systematically study energy dependence of time dynamics and to investigate a
mobility edge in our system.

##  I Introduction

Many-body localization (MBL) is a robust way of ergodicity breaking in
interacting many-body systems. The seminal early works Gornyi  _et al._
(2005); Basko  _et al._ (2006); Oganesyan and Huse (2007) show that such a
behavior is associated with the presence of strong disorder in the system.
Hundreds of studies in the last 15 years addressed various aspects of MBL so
we refer the reader to reviews of the problem Nandkishore and Huse (2015);
Luitz and Lev (2017); Alet and Laflorencie (2018); Abanin  _et al._ (2019);
Gopalakrishnan and Parameswaran (2020) Despite the numerous works, even the
issue of the existence of MBL in thermodynamic limit is not settled yet. On
one side there is a proof of existence MBL in certain one-dimensional (1D)
spin systems Imbrie (2016a, b) under assumptions regarding properties of
energy levels. On the other side recent numerical studies put into question
Šuntajs  _et al._ (2019) the very existence of MBL as a stable dynamical phase
of matter or point out difficulties in its unambiguous verification Panda  _et
al._ (2020); Sierant  _et al._ (2020a). Subsequent works provided further
arguments against the stability of MBL phase Kiefer-Emmanouilidis  _et al._
(2020); Sels and Polkovnikov (2020), or supported the scenario of stable MBL
phase Luitz and Lev (2020) indicating that earlier predictions Luitz  _et al._
(2015) of critical disorder strengths may underestimate the real value Sierant
_et al._ (2020b) (see also Doggen  _et al._ (2018); Chanda  _et al._ (2020)).
Those works concentrate mainly on the model of spin$-1/2$ systems with short
range interactions and a uniform random disorder present.

Instead, the experimental attempts to observe signatures of MBL often employ a
quasiperiodic potential Choi  _et al._ (2016); Lüschen  _et al._ (2017, 2018)
realized by superimposing two incommensurate optical lattices, the primary one
holding the atoms (allowing for the tight binding approximate description) and
the secondary one producing the disorder by weakly perturbing local minima of
the primary lattice. The quasiperiodic potential leads to MBL similarly to the
uniform random disorder. It is, however, claimed that the universality class
of the transition from delocalized, extended phase to MBL phase is dependent
on the disorder type Khemani  _et al._ (2017); Maksymov  _et al._ (2020).
Similarly, the range of interactions plays a role as well as the dimension of
the problem. We shall restrict ourselves to 1D case, still here, the existence
of MBL seem to be dependent on the interaction range. In particular for power-
law decaying interactions $V\propto 1/r^{\alpha}$ where $r$ is the distance
between interacting particles it is claimed that MBL exists for $\alpha>1$
Burin (2015a, b). This condition is naturally broken in a system with cavity
mediated all-to-all interactions Sierant  _et al._ (2019). The results
presented by two of us in this work, while valid for a finite system may
nevertheless be fragile in the thermodynamic limit as pointed out by Maksymov
and Burin (2020).

The aim of this work is to investigate ergodicity breaking in the presence of
all-to-all cavity mediated interactions in the quasi-periodic disorder
resulting from the mismatch between the cavity mode wavelength and the
wavelength of the optical lattice parallel to the cavity axis. Hence, contrary
to Sierant  _et al._ (2019) we do not assume any additional random or quasi-
periodic on-site disorder. We shall show that while an MBL regime may be seen
for a finite size system, the thermodynamic limit is subtle and the results
depend on the details of the model implementation. The paper is organized as
follows. The model and its Hamiltonian are discussed in the next Section
followed by the presentation of spectral properties of the system. Then we
discuss the time dynamics and consider different possible quenches in the
system. As we shall see some care must be taken in order to excite the system
in the controlled way.

##  II Hamiltonian

Figure 1: Quasi-periodic infinite-range interactions can be realised by a
quantum gas of bosons in an optical lattice (solid, grey line) and strongly
coupling to a standing-wave mode of a high-finesse optical resonator (dashed
blue line). The dipolar transitions couple dispersively with the cavity mode
and with a transverse laser beam. As a result the multiply scattered photons
effectively mediate atom-atom interactions with the periodicity of the cavity
mode wavelength. We assume that the ratio between the periodicity of the
optical lattice and of the cavity mode is incommensurate.

We consider $N$ bosons in a one-dimensional lattice with $K$ sites and open
boundary conditions. Their dynamics is described by the extended Bose-Hubbard
model with infinite-range interactions, whose Hamiltonian reads:

$\hat{H}=\hat{H}_{\text{BH}}+\hat{H}_{\text{cav}}\,.$ (1)

Here, the first term $\hat{H}_{\text{BH}}$ describes the standard Bose-Hubbard
dynamics with nearest-neighbour hopping at rate $J>0$ and repulsive onsite
interaction with amplitude $U$:

$\hat{H}_{\text{BH}}=-J\sum_{i=1}^{K-1}(\hat{b}_{i}^{\dagger}\hat{b}_{i+1}+\text{c.c.})+\frac{U}{2}\sum_{i}^{L}\hat{n}_{i}(\hat{n}_{i}-1)\,,$
(2)

where $\hat{b}_{i}$ and $\hat{b}_{i}^{\dagger}$ are the annihilation and
creation operators for a boson on site $i$, and
$\hat{n}_{i}=\hat{b}_{i}^{\dagger}\hat{b}_{i}$ is the number of particle at
site $i$. The infinite-range interactions is given by the Hamiltonian Habibian
_et al._ (2013a):

$\hat{H}_{\text{cav}}=-\frac{U_{1}}{K}\quantity(\sum_{i=1}^{K}{\mathcal{Z}}_{\beta}(i,\phi_{0})\hat{n}_{i})^{2}\,.$
(3)

where $U_{1}$ scales the interaction amplitude, while $\mathcal{Z}_{\beta}$ is
a function of the site $i$ and of the real parameters $\beta$ and $\phi_{0}$,
that is bounded to unity, $|\mathcal{Z}_{\beta}|\leq 1$. This function has
periodicity $\beta$ in units of the lattice periodicity. In the Appendix we
provide its specific dependence on the parameters of the setup of Fig. 1.
$\phi_{0}$ is an arbitrary phase factor corresponding to different
realisations of the effective “disorder”. In this work we choose $\beta$ to be
an irrational number, and in particular $\beta=(\sqrt{5}-1)/2$. As a result
Hamiltonian (1) is quasiperiodic.

In the rest of this manuscript we assume that $U_{1}>0$. For this choice, the
interaction term favours quasiperiodic density distributions which maximize
the expectation value of Eq. (3) and the quantum ground state exhibits the
features of a Bose glass Habibian  _et al._ (2013a, b). In these works it was
shown that Hamiltonian (1) can be realised in a cavity quantum electrodynamics
setup by overlapping the cavity mode to an optical lattice, tightly confining
the atoms along the cavity axis. The setup is illustrated in Fig. 1. In this
configuration parameter $\beta$ is proportional to the ratio between the
cavity wave number $k$ and the optical lattice wave number $k_{l}$, namely,
$\beta=k/(2k_{l})$. It is known Doggen and Mirlin (2019) that the value of
$\beta$ affects the distribution of on-site energies between neighboring sites
so the results quantitatively depend on the $\beta$ choice.

We finally remark that the single-particle ground state and spectrum of these
dynamics have been studied in Refs. Rojan  _et al._ (2016); Major  _et al._
(2018). Here, it was shown that the ground state can be localized for
sufficiently large $U_{1}$. For these values the spectrum exhibits mobility
edges.

In the rest of this manuscript we analyse the dynamics for $J=U$ and report
the interaction strength $U_{1}$ in units of $J$.

##  III Spectral properties

The localisation properties of our system are studied using exact
diagonalisation (ED) of the Hamiltonian of a system with a unit filling for
different number of sites $K$ ranging from 7 to 9. These modest system sizes
are determined by necessity of including the possibility to have several
bosons at single sites and correspond to dimensions of Hilbert space between
1716 and 24310. All quantities obtained are also averaged over 500-30000
random $\phi_{0}$ phases - corresponding to independent realizations of quasi-
periodic disorder.

Figure 2: The density of states (DOS) for $K=7,9$ and $U_{1}=5,25$ as
indicated in the figure.

We start our analysis with the density of states (DOS) and its dependence on
the system size, $K$ and cavity mediated interaction strength $U_{1}$, shown
in Fig. 2. We rescale eigenenergies $E$ to a unit interval as
$\epsilon=(E-E_{\text{min}})/(E_{\text{max}}-E_{\text{min}})$ with
$E_{\text{min}}$, $E_{\text{max}}$ denoting a minimal and maximal energy from
each disorder realization, respectively. For moderate value of all-to-all
interactions, $U_{1}=5$, the maximal DOS occurs near the middle of the
spectrum. For the higher value $U_{1}=25$ the maximum of DOS is shifted to
higher energies, with low and middle energy states laying in the tail of DOS.
When the system size is increased the shift of DOS to slightly lower energies
is observed. The strong asymmetry of DOS reflects itself in localization
properties of eigenstates in different parts of the spectrum as we show below.

Figure 3: The dependence of the mean gap ratio $\bar{r}$ on the all-to-all
interaction strength $U_{1}$ for different energy ranges: a)
$\epsilon\in[0.3,0.4]$, b) $\epsilon\in[0.6,0.7]$.

To study the localization properties of eigenstates we utilize the mean gap
ratio Oganesyan and Huse (2007). It is defined as

$\bar{r}=\left<\min\quantity{\frac{\epsilon_{i+1}-\epsilon_{i}}{\epsilon_{i}-\epsilon_{i-1}},\frac{\epsilon_{i}-\epsilon_{i-1}}{\epsilon_{i+1}-\epsilon_{i}}}\right>,$
(4)

where the averaging is performed in a band of eigenenergies around a given
rescaled energy $\epsilon$ and over disorder realizations. For ergodic system,
with level statistics well described by Gaussian Orthogonal Ensemble of random
matrices, $\bar{r}\approx 0.53$ while $\bar{r}\approx 0.39$ indicates
localisation and Poissonian level statistics Atas  _et al._ (2013). The
dependence $\bar{r}$ on the interaction strength $U_{1}$ for different parts
of spectrum and system sizes is shown in Fig.3. The $\bar{r}$ plots show a
couple of interesting properties of our system. Firstly, for a fixed system
size and energy range, the system undergoes a crossover between ergodic and
localized regimes as we increase $U_{1}$, meaning that stronger interactions
actually prohibit thermalization. Thus, the interaction induced randomness to
some extent plays a role of a (quasi)random on-site potential in the
interacting system. A similar phenomenon occurs for bosons interacting by
random contact interactions Sierant  _et al._ (2017, 2017); Sierant and
Zakrzewski (2018). Here, we keep the uniform on-site interactions fixed at low
$U=1$ value and the randomness comes from infinite-ranged cavity mediated
interactions only. A second observation is that the higher the energy is, the
stronger $U_{1}$ is necessary to induce the crossover to localized phase.
Therefore, for a fixed $U_{1}$, we observe a so-called mobility edge, namely
the lower energetically part of the spectrum is localized, while the higher is
thermal (delocalized). However, Fig. 3 b) shows that the value of $U_{1}$ at
which the system enters into the localized regime increases with system size
$K$, implying that the mobility edge shifts to lower energies with increasing
$K$. This suggests that the mobility edge may not survive in the thermodynamic
limit and in this limit all eigenstates are extended.

Figure 4: The dependence of $\bar{r}$ on $V=U_{1}/K$ for various energy
ranges. a) $\epsilon\in[0.3,0.4]$, b) $\epsilon\in[0.6,0.7]$.

Observe, however, that we can choose a different scaling the coupling strength
in $\hat{H}_{\text{cav}}^{(0)}$ term, (3). The factor in front of the sum,
$U_{1}/K$, has been chosen, in accordance with the customary approach Habibian
_et al._ (2013a); Dogra  _et al._ (2016) that ensures the appropriate
extensive scaling of the Hamiltonian with increasing $K$. In fact, it may be
justified by assuming that the increase of the system size $K$, for a fixed
wavenumber $k$, implies the increase in the cavity length and thus the
appropriate scaling of the cavity mode volume. However, another approach is
possible - to keep the cavity length fixed and change the lattice spacing.
Then, instead of the coefficient $U_{1}/K$ in front of the all-to-all
interactions term, we would have a constant $V$ that remains fixed when $K$ is
increased. Such a scaling implies a nonstandard thermodynamic limit as the
energy becomes a super extensive quantity. A similar procedure was adopted in
Gopalakrishnan and Huse (2019) to study MBL in presence of power-law
interactions. We show the mean gap ratio $\bar{r}$ as a function of $V$ in
Fig. 4. For a sufficiently large $V$ a crossover to localized regime is
observed. However, contrary to fixed $U_{1}$ case, the resulting
$\overline{r}(W)$ curves deviate from the ergodic value $\overline{r}\approx
0.53$ at a system size independent value of $V$ and the crossover becomes
steeper with increasing system size. This suggests that the observed crossover
becomes a sharp transition to an MBL phase at a fixed critical value of $V$ in
the considered non-standard thermodynamic limit. We note that a similar
behavior at the ergodic-MBL crossover was observed for disordered $U(1)$
lattice gauge theory Giudici  _et al._ (2020), in which the elimination of
gauge field leads to long-range interactions.

Let us note that identifying the transition for $V=const.$ limit implies the
delocalized limit for $U_{1}=const.$ for arbitrary $U_{1}$. On the other hand
all experiments are performed for a finite $K$ value that may reach several
tens or hundreds but necessarily remains finite so our observations for the
standard $U_{1}=const$ scenario are experimentally relevant.

Figure 5: The $\epsilon$ dependence of von Neumann biparite entanglement
entropy $S{/S_{0}}$ of eigenstates normalized by a random Gaussian state value
$S_{0}$. Subsequent panels correspond to different values of $V=U_{1}/K$: a)
$V=0.5$, b) $V=1$, c) $V=2$, d) $V=3$.

The crucial difference between extended and localized states can be seen in
their entanglement properties. For extended states quantum entanglement
between subsystems which we call $A$ and $B$ resulting from splitting of the
entire system into two, follows a volume law, while for localized states area
law is expected. To quantify the entanglement, we calculate the von Neumann
biparite entanglement entropy defined as

$S=-\Tr_{B}[\hat{\rho}_{B}\ln\hat{\rho}_{B}],$ (5)

where $\hat{\rho}_{B}=\Tr_{A}[\hat{\rho}_{AB}]$ is a (usually mixed) state of
the subsystem $B$ obtained by tracing the state
$\hat{\rho}_{AB}=\ket{\psi}\bra{\psi}$ over degrees of freedom of
complementary subsystem $A$. If subsystems $A$ and $B$ are not entangled,
$\hat{\rho}_{B}$ is a pure state and $S=0$. Fig. 5 presents the entanglement
entropy of eigenstates ($\hat{\rho}_{AB}=\ket{\psi_{n}}\bra{\psi_{n}}$, where
$\hat{H}\ket{\psi}=E_{n}\psi_{n}$) averaged over eigenstates with energies
corresponding to the rescaled energy $\epsilon$ for different values of $V$
and $K$, normalized by the entropy $S_{0}$ of random Gaussian states Vidmar
and Rigol (2017). In this case, $A$ and $B$ subsystems are two equal $K/2$
parts.

Again we observe two main features. Firstly, for low lying energy states the
normalized entropy is close to 0 and independent of system size $K$,
indicating localized states. However for higher energies $\epsilon$, the
entanglement entropy $S$ approaches $S_{0}$ as expected for thermal states.
Interestingly, for very high energies we observe the decrease in $S$, but this
appears at the edge of the spectrum where the density of states is again very
low. Secondly, increasing $V$ shifts the the thermal region to higher
energies, which is consistent with level statistics. What is important, for
constant $V$ and changing $K$ we again observe the crossing of curves
corresponding to different system sizes. Interestingly, the crossing is
visible on both sides of $S$ maximum, suggesting that the localized
eigenstates, fulfilling the area law, found at the highest energies $\epsilon$
can also survive in the non-standard thermodynamic limit.

Figure 6: Population imbalance $I(t)$ during time evolution of a density-wave
like state for $K=8$ and a several of $U_{1}$ values: 25, 15, 10, 5, 1
respectively, starting from the top.

##  IV Time dynamics

In the previous Section we analysed spectral properties discussing various
aspects of the crossover between ergodic and MBL regimes as well as the
presence of the mobility edge. Here, we address more experimentally relevant
issue - whether the localization properties may be observed in the time
dynamics starting from a judiciously prepared initial state. Following the
usual experimental strategy Choi  _et al._ (2016); Lüschen  _et al._ (2017) we
might consider the population imbalance, defined as

$I(t)=\left<\frac{1}{K}\sum_{i=1}^{K}(-1)^{i}\hat{n}_{i}(t)\right>.$ (6)

An initial state with an imprinted density pattern, for which the imbalance
$I(t)$ is non-zero at $t=0$ can be used to demonstrate ergodicity breaking.
For unit density assumed by us such a state can be chosen as
$\ket{2020\dots}$. The system evolves governed by $\hat{H}$ and the $I(t)$ is
monitored to reveal the system properties. The initial unit imbalance decays
to 0 for systems obeying the eigenstate thermalization hypothesis Deutsch
(1991); Srednicki (1994). In such a case the information about the initial
configuration is lost during the course of time evolution. In contrast, a non-
zero long-time value of $I$ suggests that the information is partially
preserved and is indicative of the localized phase. The evolution of the
population imbalance for a several values of $U_{1}$ is shown in Fig. 6. The
results are obtained using Chebyshev propagation Tal-Ezer and Kosloff (1984);
Fehske and Schneider (2008) relying on sparse matrix operations, which allows
us to reach system sizes up to $K=12$ with reasonable computational time and
resources. Surprisingly, although the level statistics (Figs. 3,4) indicate
localization in some ranges of $U_{1}$, for all parameters studied here, the
imbalance for the initial density wave state rapidly decays to zero suggesting
thermalization. It is easily understood after inspecting the energy of the
density wave state $\ket{2020...}$ (energy is conserved during propagation)
which, for all reasonable choices of the parameters of the problem,
corresponds to regions identified as ergodic by the mean gap ratio.

Figure 7: The dependence of average energy $\epsilon$ of the quenched state
w.r.t. the final Hamiltonian on several quenched parameters for $K=8$. On both
panels in each row the same energy values are presented, but with different
types of errors represented by vertical lines. In the left column they
represent the quantum error
$\sqrt{\expectationvalue{\hat{H}_{\text{final}}^{2}}-\expectationvalue{\hat{H}_{\text{final}}}^{2}}$,
while in the right column – the standard deviation of the random phase
distribution of mean energies.

One can imagine a preparation of a different initial product state, e.g. a
uniformly occupied state $\ket{1111\dots}$ and use the transport distance
Rispoli  _et al._ (2019) (see also below) to look for signatures of
localization. Our attempts in that direction failed (except for highly unusual
states which did not allow us to study the breaking of ergodicity in our
system in a systematic manner). The reason for that is simple. To reach the
MBL regime in our system we need a sufficiently strong disorder which appears
solely in the cavity mediated interaction term. For large $U_{1}$ this term
dominates the energy of product states. Thus, the occupation numbers of an
initial product state must be adapted in such a way that the all-to-all
interaction term is possibly small in order to reach localized regime.

To inspect the low energy regime in a more systematic manner we employ an
alternative approach introduced by Naldesi  _et al._ (2016) and termed the
quantum quench spectroscopy (QQS). In the first step of QQS the system is
prepared in the ground state of some initial Hamiltonian
$\hat{H}_{\text{ini}}$. An abrupt change of one of its parameters to the
desired $\hat{H}$ follows. In such a procedure the initial state is not
changed but now becomes a non-stationary initial superposition of eigenstates
of $\hat{H}$. One may hope that an appropriate choice of the initial
Hamiltonian $\hat{H}_{\text{ini}}$ may allow to probe different parts of the
spectrum of the final Hamiltonian $\hat{H}$. It was shown in Naldesi  _et al._
(2016) that this procedure, when ramping the disorder strength, may allow one
to probe the mobility edge in the system of interacting spinless fermions in a
quasiperiodic potential.

It is important, however to realize two possible limitations. As the initial
state is not an eigenstate of $\hat{H}$ its energy suffers from quantum
uncertainty

$\Delta
E=\sqrt{\expectationvalue{\hat{H}^{2}}-\expectationvalue{\hat{H}}^{2}},$ (7)

the fate shared, in fact, also by the initial product states. Moreover, as a
target mean energy depends on a random $\phi_{0}$ phase in (3), it is also a
random variable, whose dispersion can be characterized by a sample standard
deviation. This second uncertainty is known to decrease with an increasing
system size Naldesi  _et al._ (2016), so it is not a major concern in the
thermodynamic limit. On the other hand it may be an obstacle for relatively
small system sizes treated by ED in this work.

Comparison of mean energies obtained and both uncertainties is shown in Fig. 7
for different possible quenche scenarios. We fix $\beta=(\sqrt{5}-1)/2\approx
0.618$ and $U_{1}=25$ as the final values. In the first approach we quench the
disorder amplitude from some initial $(U_{1})_{0}$ value to $U_{1}$ \- the
path similar to Naldesi  _et al._ (2016). The results presented in the upper
row of Fig. 7 show low energy selectivity. For small difference between the
initial and final disorder amplitude low excitations are reached only. The
rapid jump in relative energy $\epsilon$ appears for low initial interaction
strength. When $(U_{1})_{0}$ has a different sign than $U_{1}$ a relatively
large $\epsilon$ values may be reached but with quite significant quantum
uncertainty. Similarly, the region of small $(U_{1})_{0}0$ where a rapid
change of final energies is observed suffers from that uncertainty.

Situation is entirely different if, instead of the interaction strength, the
ratio of cavity and laser wavenumbers is changed, i.e. $\beta$ as shown in the
middle row of Fig. 7. The quantum uncertainties are small and a broad range of
final energies may be achieved, however $\phi_{0}$ averaging uncertainty is
very large for the system sizes under consideration. Moreover, the maximal
energy achieved is around $\epsilon\approx 0.70$, which is below the mobility
edge estimated for $K=8$ from gap ratio as $\epsilon\approx 0.75-0.80$ not
enabling us to probe the crossover between ergodic and MBL regimes.

Finally, the bottom row in Fig. 7 shows yet another possible quench. For each
(random) $\phi_{0}$ we choose a different $\phi_{0}(0)$ and define the quench
as a rapid change of the phase from $\phi_{0}(0)$ to $\phi_{0}$. We denote the
magnitude of the phase change by $\Delta\phi_{0}$. As observed the quantum and
statistical uncertainties are relatively small. Moreover, the final energy
approaches 1 for $\Delta\phi_{0}=\pi/2$. Thus, we perform further the QQS
simulations using this type of quench.

Figure 8: a) Time evolution of transport distance for $K=8$ and a several
$\Delta\phi_{0}$ values. b) A similar evolution of fluctuations $F$. c)
Maximal transport pseudo-velocity dependence on initial state energy for a
several $K$. d) Fluctuations pseudo-velocity energy dependence. e) Average
value of $\Delta x$ for $t\in[100,500]$ energy dependence. f) Average value of
$\Delta x$ for the same time interval energy dependence.

After a change of the cavity wavenumber by $\Delta\phi_{0}$ we consider time
evolution of the system. Instead of the Imbalance we calculate, following
Rispoli  _et al._ (2019) the transport distance, defined as

$\Delta
x=2\absolutevalue{\sum_{d=1}^{K-1}d\times\overline{\expectationvalue{G_{c}^{(2)}(i,i+d)}_{i}}},$
(8)

where
$G_{c}^{(2)}(i,i+d)=\expectationvalue{\hat{n}_{i}\hat{n}_{i+d}}-\expectationvalue{\hat{n}_{i}}\expectationvalue{\hat{n}_{i+d}}$
and averaging over sites $\expectationvalue{\dots}_{i}$ is done before
averaging over realizations $\overline{(\dots)}$. Note that the measurement of
$\Delta x$ requires access to two-site correlations available in state of the
art experiments Rispoli  _et al._ (2019). The second quantity we calculate is
an edge fluctuation defined as

$F=\frac{\quantity(\expectationvalue{\hat{n}_{1}^{2}}-\expectationvalue{\hat{n}_{1}}^{2})+\quantity(\expectationvalue{\hat{n}_{K}^{2}}-\expectationvalue{\hat{n}_{K}}^{2})}{2}.$
(9)

It requires a site dependent resolution at the edges of the system but does
not involve two-point correlations. Both quantities were analyzed previously
with the aim of identifying the mobility edge in a Bose-Hubbard model in a
tilted lattice Yao and Zakrzewski (2020). Time evolution of both quantities
for $K=8$, $U_{1}=25$ and different initial state energies are depicted in the
first row of Fig. 8. Initial states were prepared in quantum quench of
$\phi_{0}$, with $\Delta\phi_{0}=0.50,1.00,1.15,1.35$ and $1.50$, where two
first values yield $\epsilon\approx 0.21$ and $0.66$ in the localized regime,
third gives $\epsilon\approx 0.78$ in the vicinity of mobility edge, fourth
gives $\epsilon\approx 0.89$ which is well in extended regime and fifth –
$\epsilon\approx 0.94$, which is near the high end of spectrum, where the
states localize again. Observe that $\Delta x(t)$ evolution starts with an
abrupt, ballistic-like transient behaviour, but after a time of the order of
inverse tunneling energy $1/J=1$ one can observe a slower, subdiffusive
growth. When the Heisenberg time scale (which is proportional to the inverse
of the mean level spacing) is reached, the transport distance value saturates.
For energies corresponding to the localized part of the spectrum, $\Delta x$
saturates at low values below a single site unit distance while for those in
the thermal domain $\Delta x$ is of the order of system size $K=8$. Here, in
fact, the saturation is reached well before the Heisenberg time when the
system size is reached. In the vicinity of the crossover between ergodic and
MBL regimes the saturation of $\Delta x$ occurs only later and at an
intermediate value (of the order of unity). Notice, that $\Delta x$ for the
highest energy considered $0.94$ saturates at lower value then for $0.89$
confirming the presence of the localized regime at the high end of the
spectrum. We observe a similar situation for the edge fluctuation, $F$. The
initial transient behavior is followed by a typical, logarithmic growth
followed by a saturation of $F(t)$.

It was proposed in Yao and Zakrzewski (2020) to detect the mobility edge using
the notion of pseudovelocity. The authors computed the ratio

$V_{O}=\frac{O(t_{\text{max}})-O(t_{\text{min}})}{t_{\text{max}}-t_{\text{min}}}$
(10)

where $O(t)=\Delta x$ or $F$ for fixed times $t_{\text{max}}$ and
$t_{\text{min}}$ for evolution of states with different energies and plotted
them against $\epsilon$. They noticed a peak of this value near mobility edge,
providing an experimental way for its observation. However, in our case, the
mobility edge is relatively high in the rescaled energy $\epsilon$ hindering
the observation of a dip of pseudovelocity after the initial rise. Also, for
large $K$ and maximal simulation time $t=10^{3}$ some observables are not yet
saturated. In effect the pseudovelocities obtained had large fluctuations. In
order to mitigate that issue, we modified the approach – instead of fixing
$t_{\text{max}}$ and $t_{\text{min}}$, we scanned $\Delta x(t)$ and $F(t)$
dependencies with $\Delta t=30$ window starting from $t=10$, computed a linear
fit slope in this window and chose its maximal value as a pseudovelocity. This
allowed us to make $V_{O}(\epsilon)$ dependencies much smoother (Fig. 8c, d).
For both $V_{x}$ and $V_{F}$, on the localized side of spectrum the values are
near 0 and start to visibly rise around the transition energy $\epsilon\approx
0.75$. They also correctly recognize the shift of transition energy to lower
values for larger system sizes. Experimental accessibility of $\Delta x$ is
lower than $F$, because it requires measuring all the correlations, while the
latter one concerns the sites on the edges only. On the other hand $V_{F}$ is
more affected by fluctuations than $V_{x}$. Instead of pseudovelocities one
may, however, consider simply the average values of observables in a fixed
time interval as shown in the bottom row of Fig. 8. Their behavior is much
smoother, providing a similar observation – a rapid increase of both the mean
transport distance and mean fluctuation at the crossover from localized to
extended regime.

Figure 9: a)-c) Time evolution of von Neumann biparite entanglment entropy for
$\Delta\phi_{0}$-quenched states for, respectively, $K=8,10,12$. d) Maximal
value (for $t=10^{3}$) of $S(t)/S_{0}$ as a function of energy $\epsilon$ for
those 3 system sizes.

We may also consider time evolution of the entanglement entropy using the same
quenched initial states. For initially separable state the hallmark of MBL is
the logarithmic entanglement entropy growth Žnidarič  _et al._ (2008);
Bardarson  _et al._ (2012); Serbyn  _et al._ (2013). This has been typically
tested for separable, product-like initial states. It is interesting,
therefore, to observe that entropy growth also for our initial entangled
states - see Fig. 9. For localized initial states $S$ grows very slowly. In
the vicinity of the crossover we observe a delayed saturation, while in the
extended regime $S$ saturates at the value of the order of the average
entanglement entropy of a random Gaussian state $S_{0}$. However, $S_{0}$ is
never reached for system sizes and times under consideration. One may notice,
that, once again, the saturation value for highest energy $\epsilon=0.94$ is
lower that for $\epsilon=0.89$, which is consistent with previous results. We
also plot, in Fig. 9(d), the dependence of maximal $S_{\text{max}}$ against
initial state energy. In this case, the visible rise of $S_{\text{max}}$
starts around $\epsilon\approx 0.6$, below the transition energy.
Nevertheless, it provides additional evidence supporting the mobility edge, at
least for system sizes studied in this paper.

##  V Summary

We considered signatures of many body localization in the cavity QED model
combined with the optical lattice whose wavelength is incommensurate with the
cavity mode wavelength. In such a case the effective long-range interactions
between atoms (bosons) in the system become dependent on the
incommensurability ratio resulting in the effective quasi-random potential.
Importantly, we do not assume any other type of disorder to be present.

We analyse both spectral properties of the system (mean gap ratio as well as
entanglement entropy of eigenstates) as well as the time dynamics from
specially prepared intitial states. For finite (small, amenable to exact
diagonalization) systems we observe clear signatures of MBL with mean gap
ratio reaching the Poisson level characteristic for a localized regime.
However, the interaction strength leading to localization seems to increase
with the system size. In effect, in the standard thermodynamic limit, we
predict that the system will remain delocalized for almost all the states. One
may define, however, a nonstandard thermodynamic limit, in which the
interaction strength is not scaled with system size. While this makes the
energy of the system a super-extensive quantity, it also gives rise to a clear
mobility edge between localized and delocalized regimes.

These findings, based on spectral statistics, were confirmed and verified by
time dynamics in which we addressed quantities that are within the reach of
current experiments. Interestingly, meaningful studies required quantum quench
spectroscopy with initial states obtained by a parameter quench. We have shown
that a great care should be taken in the choice of the quenched parameter, in
particular for necessarily small system sizes considered. We have also shown
that in the localized regime the entanglement entropy growth for such initial
weakly entangled states follows a logarithmic law, a hall mark of MBL observed
typically for initial product states.

Our results indicate that many-body cavity quantum electrodynamics setups are
well suited to studies of quantum ergodicity and its breaking. The setup
considered in this work is considerably simpler than the system investigated
in Sierant  _et al._ (2019). Here, the cavity mediated interactions act
effectively as a source of randomness leading to a non-trivial interplay
between ergodic and MBL eigenstates at different energies.

## Acknowledgements

The support by National Science Centre (Poland) under project OPUS
2019/35/B/ST2/00034 (J.Z.) is acknowledged. P.K. would like to acknowledge the
support of Ministry of Science and Higher Education (Poland) grant no.
0108/DIA/2020/49. P.S. acknowledges the support of Foundation for Polish
Science (FNP) through scholarship START. G. M. acknowledges the support by the
German Research Foundation (Project-ID 429529648 – TRR 306 QuCoLiMa "Quantum
Cooperativity of Light and Matter) and by theGerman Ministry of Education and
Research (BMBF), QuantERA project NAQUAS. Project NAQUAS has received funding
from the QuantERA ERA-NET Co-fund in Quantum Technologies implemented within
theEuropean Union?s Horizon2020 program. Numerical simulations were carried
out with the support of the Interdisciplinary Center for Mathematical and
Computational Modeling (ICM) at the University of Warsaw under grant no.
GB76-1. For linear algebra computations, Armadillo C++ language bindings
Sanderson and Curtin (2016, 2018) with Intel® Math Kernel Library as BLAS and
LAPACK backend were used.

##  Appendix A Details on the Bose-Hubbard model

We consider the setup of Fig. 1 in the regime where the cavity mode dynamics
is characterized by a faster time scale than the atomic motion. In this limit
one can apply a coarse graining to the time evolution Habibian  _et al._
(2013a); Sierant  _et al._ (2019). Eliminating the cavity degrees of freedom
leads to an effective infinite-range interaction term Schütz  _et al._ (2015).
By assuming that the atoms are tightly bound in the lowest band of an optical
lattice with wave number $2k_{l}$, the interaction term can be cast in the
form of Eq. (3), with Habibian  _et al._ (2013a); Dogra  _et al._ (2016);
Sierant  _et al._ (2019)

$\displaystyle\mathcal{Z}_{\beta}(i,\phi_{0})=\int\differential{x}\cos(kx+\phi_{0})w_{i}(x)w_{j}(x)\,,$
(11)

where $w_{i}(x)$ is lowest band Wannier function centred at $i^{\text{th}}$
site. The other quantities are the wave vector of the cavity mode $k$, the
interaction strength $U_{1}$, and an arbitrary phase $\alpha_{0}$. For a
sufficiently deep optical lattice potential, we approximate
$w_{i}(x)\approx\delta(x-\pi/k_{l}(i+1/2))$ and obtain:

$\mathcal{Z}_{\beta}(i,\phi_{0})\approx\cos(2\pi\beta i+\phi_{0})\,,$ (12)

with $\beta=k/(2k_{l})$.

This approach assumes that the cavity dynamics adiabatically follows the
atomic motion. In particular the coarse-grained cavity field operator is now
proportional to the atomic observables

$\hat{a}\propto\sum_{i=1}^{K}\mathcal{Z}_{\beta}(i,\phi_{0})\hat{n}_{i}\,.$
(13)

This quantity can be revealed by homodyne detection of the field at the cavity
output Black  _et al._ (2003); Baumann  _et al._ (2011).

## References

  * Gornyi  _et al._ (2005) I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko  _et al._ (2006) D. Basko, I. Aleiner, and B. Altschuler, Ann. Phys. (NY) 321, 1126 (2006).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015), https://doi.org/10.1146/annurev-conmatphys-031214-014726 .
  * Luitz and Lev (2017) D. J. Luitz and Y. B. Lev, Annalen der Physik 529, 1600350 (2017).
  * Alet and Laflorencie (2018) F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018), quantum simulation / Simulation quantique.
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Gopalakrishnan and Parameswaran (2020) S. Gopalakrishnan and S. Parameswaran, Physics Reports 862, 1 (2020), dynamics and transport at the threshold of many-body localization.
  * Imbrie (2016a) J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016a).
  * Imbrie (2016b) J. Z. Imbrie, Journal of Statistical Physics 163, 998 (2016b).
  * Šuntajs  _et al._ (2019) J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, arXiv e-prints , arXiv:1905.06345 (2019), arXiv:1905.06345 [cond-mat.str-el] .
  * Panda  _et al._ (2020) R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Europhysics Letters 128, 67003 (2020).
  * Sierant  _et al._ (2020a) P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. Lett. 124, 186601 (2020a).
  * Kiefer-Emmanouilidis  _et al._ (2020) M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Phys. Rev. Lett. 124, 243601 (2020).
  * Sels and Polkovnikov (2020) D. Sels and A. Polkovnikov, “Dynamical obstruction to localization in a disordered spin chain,” (2020), arXiv:2009.04501 [quant-ph] .
  * Luitz and Lev (2020) D. J. Luitz and Y. B. Lev, Phys. Rev. B 102, 100202 (2020).
  * Luitz  _et al._ (2015) D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
  * Sierant  _et al._ (2020b) P. Sierant, M. Lewenstein, and J. Zakrzewski, “Polynomially filtered exact diagonalization approach to many-body localization,” (2020b), arXiv:2005.09534 [cond-mat.dis-nn] .
  * Doggen  _et al._ (2018) E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Phys. Rev. B 98, 174202 (2018).
  * Chanda  _et al._ (2020) T. Chanda, P. Sierant, and J. Zakrzewski, Phys. Rev. B 101, 035148 (2020).
  * Choi  _et al._ (2016) J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016).
  * Lüschen  _et al._ (2017) H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Phys. Rev. Lett. 119, 260401 (2017).
  * Lüschen  _et al._ (2018) H. P. Lüschen, S. Scherg, T. Kohlert, M. Schreiber, P. Bordia, X. Li, S. Das Sarma, and I. Bloch, Phys. Rev. Lett. 120, 160404 (2018).
  * Khemani  _et al._ (2017) V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev. Lett. 119, 075702 (2017).
  * Maksymov  _et al._ (2020) A. Maksymov, P. Sierant, and J. Zakrzewski, “Many-body localization in one dimensional optical lattice with speckle disorder,” (2020), arXiv:2008.00219 [cond-mat.dis-nn] .
  * Burin (2015a) A. L. Burin, Phys. Rev. B 91, 094202 (2015a).
  * Burin (2015b) A. L. Burin, Phys. Rev. B 92, 104428 (2015b).
  * Sierant  _et al._ (2019) P. Sierant, K. Biedroń, G. Morigi, and J. Zakrzewski, SciPost Phys. 7, 8 (2019).
  * Maksymov and Burin (2020) A. O. Maksymov and A. L. Burin, Phys. Rev. B 101, 024201 (2020).
  * Habibian  _et al._ (2013a) H. Habibian, A. Winter, S. Paganelli, H. Rieger, and G. Morigi, Phys. Rev. Lett. 110, 075304 (2013a).
  * Habibian  _et al._ (2013b) H. Habibian, A. Winter, S. Paganelli, H. Rieger, and G. Morigi, Phys. Rev. A 88, 043618 (2013b).
  * Doggen and Mirlin (2019) E. V. H. Doggen and A. D. Mirlin, Phys. Rev. B 100, 104203 (2019).
  * Rojan  _et al._ (2016) K. Rojan, R. Kraus, T. Fogarty, H. Habibian, A. Minguzzi, and G. Morigi, Phys. Rev. A 94, 013839 (2016).
  * Major  _et al._ (2018) J. Major, G. Morigi, and J. Zakrzewski, Phys. Rev. A 98, 053633 (2018).
  * Atas  _et al._ (2013) Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Phys. Rev. Lett. 110, 084101 (2013).
  * Sierant  _et al._ (2017) P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev. A 95, 021601 (2017).
  * Sierant  _et al._ (2017) P. Sierant, D. Delande, and J. Zakrzewski, Acta Phys. Polon. A 132, 1707 (2017).
  * Sierant and Zakrzewski (2018) P. Sierant and J. Zakrzewski, New Journal of Physics 20, 043032 (2018).
  * Dogra  _et al._ (2016) N. Dogra, F. Brennecke, S. D. Huber, and T. Donner, Phys. Rev. A 94, 023632 (2016).
  * Gopalakrishnan and Huse (2019) S. Gopalakrishnan and D. A. Huse, Phys. Rev. B 99, 134305 (2019).
  * Giudici  _et al._ (2020) G. Giudici, F. M. Surace, J. E. Ebot, A. Scardicchio, and M. Dalmonte, Phys. Rev. Research 2, 032034 (2020).
  * Vidmar and Rigol (2017) L. Vidmar and M. Rigol, Phys. Rev. Lett. 119, 220603 (2017).
  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Tal-Ezer and Kosloff (1984) H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984), https://doi.org/10.1063/1.448136 .
  * Fehske and Schneider (2008) H. Fehske and R. Schneider,  _Computational many-particle physics_ (Springer, Germany, 2008).
  * Rispoli  _et al._ (2019) M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai, J. Léonard, and M. Greiner, Nature 573, 385 (2019).
  * Naldesi  _et al._ (2016) P. Naldesi, E. Ercolessi, and T. Roscilde, SciPost Phys. 1, 010 (2016).
  * Yao and Zakrzewski (2020) R. Yao and J. Zakrzewski, “Many-body localization of bosons in optical lattice: Dynamics in disorder-free potentials,” (2020), arXiv:arXiv:2007.04745 .
  * Žnidarič  _et al._ (2008) M. Žnidarič, T. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
  * Bardarson  _et al._ (2012) J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * Serbyn  _et al._ (2013) M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 110, 260601 (2013).
  * Sanderson and Curtin (2016) C. Sanderson and R. Curtin, Journal of Open Source Software 1, 26 (2016).
  * Sanderson and Curtin (2018) C. Sanderson and R. Curtin, Lecture Notes in Computer Science (LNCS) 10931, 422 (2018).
  * Schütz  _et al._ (2015) S. Schütz, S. B. Jäger, and G. Morigi, Phys. Rev. A 92, 063808 (2015).
  * Black  _et al._ (2003) A. T. Black, H. W. Chan, and V. Vuletić, Phys. Rev. Lett. 91, 203001 (2003).
  * Baumann  _et al._ (2011) K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger, Phys. Rev. Lett. 107, 140402 (2011).