Entangled quantum cellular automata, physical complexity, and Goldilocks rules

Cellular automata are interacting classical bits that display diverse emergent behaviors, from fractals to random-number generators to Turing-complete computation. We discover that quantum cellular automata (QCA) can exhibit complexity in the sense of the complexity science that describes biology, sociology, and economics. QCA exhibit complexity when evolving under ‘Goldilocks rules’ that we define by balancing activity and stasis. Our Goldilocks rules generate robust dynamical features (entangled breathers), network structure and dynamics consistent with complexity, and persistent entropy fluctuations. Present-day experimental platforms—Rydberg arrays, trapped ions, and superconducting qubits—can implement our Goldilocks protocols, making testable the link between complexity science and quantum computation exposed by our QCA.

Complexity science is primarily phenomenological and so does not admit of an axiomatic definition [16].Instead, a system is regarded as complex if data analysis reveals multiple signals widely associated with complexity [17,18].Examples include natural selection [19], diversity [20] power-law statistics [14], robustness-fragility tradeoffs [21], and complex networks [17,18,22].Physical complexity has been observed at the classical and biological scales.Does it arise in the quantum dynami-cal regime?
We answer affirmatively, showing that QCA can exhibit physical complexity when subject to evolutions that we term Goldilocks rules.A Goldilocks rule changes a qubit's state when, and only when, half the qubit's neighbors occupy the +1 σz eigenstate, |1 .If too many neighbors occupy |1 , or too few neighbors, the central qubit remains static.Goldilocks rules balance evolution with stasis.This tradeoff parallels the tradeoffs in traditional complex systems [21].Hence one should expect Goldilocks rules to produce complexity.They do, we discover-even in some of the simplest QCA constructable: unitary, one-dimensional (1D) nearestneighbor circuits.
We quantify QCA's complexity using quantum generalizations of measures used to quantify the brain's complexity in electroencephalogram/functional-magneticresonance-imaging (EEG/fMRI) [24,25]: mutualinformation networks.For each pair (j, k) of qubits, we calculate the mutual information M jk between them.We regard the qubits as a graph's nodes and the M jk as the links' weights.The nodes and links form a graph, or network.Upon calculating complex-network measures, we discover behaviors traditionally associated with complexity.We complement these complex-network measures with complexity measures based on many-body quantities: persistent entropy fluctuations and emergent phe- ) if and only if its nearest neighbors are 0s [4].This L=19-site chain is initialized with a 1 centered in 0s.A periodic pattern, called a blinker, propagates upward in discrete time steps.(b) We extend C201, for QCA, to rule T1, defined in the main text.The classical bit becomes a quantum bit in an L=19-site chain.The average spin σz j has richer discrete-time dynamics, oscillating between the classical extremes of 0 (black) and 1 (white) and filling the lattice.A truly quantum analog of the classical blinker-an entangled breather-appears in Figs.2(c)-(d).(c) The quantum lattice evolves into a high-entropy entangled state.sj denotes the site-j von Neumann entropy.(d) An analog quantum computer has implemented rule T1, of which the recently named PXP model is a particular case.Optical tweezers were used to trap a Rydberg-atom chain [23].(e) A quantum circuit can realize digital T1 dynamics.One QCA time step requires two layers of controlled-controlled-Hadamard gates.The first layer evolves even-index qubits; and the second layer, odd-index.The dashed line represents the boundary qubits, which remain fixed in |0 's.
nomena.This constellation of observations, in the empirical and phenomenological spirit of complexity science, leads to our conclusion that Goldilocks QCA produce physical complexity.
Furthermore, we show that these Goldilocks rules can be implemented with two classes of quantum computers: First, analog quantum simulators [26] include Rydberg atoms in optical lattices [23,27] and trapped ions [28,29].Second, digital quantum computers include superconducting qubits, used to demonstrate quantum supremacy [30].Quantum simulators [26,31,32] have energized the developing field of nonequilibrium quantum dynamics [33].Our work establishes a guiding question for this field: What are the origins of complexity, that we find complexity in simple, abiological quantum systems?
In summary of our work's significance, we bridge the fields of quantum computation and complexity science via QCA.QCA, we discover, can exhibit physical complexity similar to that observed in biological, social, and ecological systems.We quantify this complexity and prescribe how to define QCA evolution rules that generate it.Today's Rydberg arrays, trapped ions, and superconducting qubits can realize our QCA; hence our bridge between disciplines is not only theoretical, but also experimentally accessible.This paper is organized as follows.Section II introduces our QCA evolution schemes, including Goldilocks rules.
Section III demonstrates our central claim: Goldilocks QCA generate physical complexity despite their simple unitary structure.Sections II-III suffice for grasping the novelty of Goldilocks QCA's complexity generation.Our QCA's mathematical details appear in Sec.IV.In Sec.V, we propose a physical implementation of our QCA in Rydberg systems, an experimentally realized quantum simulator [26].Section VI concludes with this work's significance.

II.
OVERVIEW OF QCA In this section, we first introduce classical elementary cellular automata and compare them with our QCA (Sec.II A).Second, we establish the terminology and notation used to label and construct our QCA (Sec.II B).For generality, and to enable simulation by discrete and analog quantum computers, we define digital and continuous-time QCA.The number of QCA models is vast.Therefore, we pare down the options to a diverse and illustrative, but manageable, set (Sec.II C).We identify the Goldilocks rules, which we contrast with trivial and near-Goldilocks rules (Sec.II D).

A. Classical elementary cellular automata
The best-known classical cellular automata are the 1D elementary cellular automata (ECA).An ECA updates each bit of a length-L string in discrete time steps.Each bit is updated in accordance with a local transition function dependent on the bit's neighborhood.The neighborhood is defined as the bit and its two nearest neighbors.The transition function is encoded in the ECA rule number, C R .The R = 0, 1, . . ., 255.C R specifies the bit's next state, given the neighborhood's current state.Section IV details the map from R to the local transition function.
During each time step, all bits are updated simultaneously.Such an update requires a copy of the current bit string: Without the copy, one updates the j th bit, based on its neighbors, then must update the neighbors based on the earlier j th bit's value.But the earlier value has been overwritten; hence the necessity of the copy.Though simple to define, ECA dynamics generate emergent features.Examples include fractal structures, high-quality random number generation [3], and Turingcompleteness (the ability to simulate any classical computer program) [2].
Our QCA partially resemble, and partially contrast with, ECA.Similarities include one-dimensionality, spatial discreteness, and local updates' affecting only a neighborhood's center site.But we trade bits for qubits, securing new features for, and restrictions on, our QCA.First, time updates can entangle qubits, generating quantum entropy absent from ECA.Second, our QCA can run in discrete or continuous time.Hence digital and analog quantum computers can implement our QCA.Third, a QCA rule number specifies the neighborhood conditions under which an operator evolves a qubit's state.In contrast, an ECA rule number specifies a bit's next state.Fourth, our QCA evolve a qubit under an operator independent of the qubit's state.Our nearest-neighbor QCA therefore parallel the 16 locally invertible (unitary) ECA.Fifth, in digital QCA, not all qubits are updated simultaneously.As explained above, a simultaneous update requires a copy of the current bit string, in classical ECA.So would a simultaneous update require a copy of the current joint quantum state, in digital QCA.The no-cloning theorem debars such a copy [5,51].Finally, we generalize the local update rule to depend on ≥ 2 neighbors, as the ECA have been generalized [57].

B. Defining quantum cellular automata
Section IV will detail our QCA mathematically.Here, we provide the key terminology and notation.σα denotes the α = x, y, z Pauli operator.The σz eigenbasis forms the computational basis: |0 denotes the eigenstate associated with the eigenvalue 1, and |1 denotes the eigenstate associated with −1.
However, our definitions conform to the most general notion of a QCA, as a dynamical quantum system that is causal, unitary, and translation-invariant [6] [58].Furthermore, our family of QCA overlaps with the family defined in [51,59].The overall structure of their QCA definition resembles the structure of our definition.Differences include our QCA's accommodating five-site neighborhoods and their QCA's accommodating nonunitary evolutions.Additionally, our family of analog QCA includes, as a special case, the QCA in [9].A QCA updates each qubit of a length-L chain according to a unitary update rule.Qubit j has a set Ω j of neighbors: the sites, excluding j, within a radius r of j.Ω j corresponds to a Hilbert space on which are defined density operators that form a vector space spanned by projectors P(i) Ωj .The i enumerates the neighbors' possible configurations, and the projectors form a complete set: The neighborhood of site j consists of j and its neighbors.
We define digital time evolution as follows.Given an ECA rule number, we convert into a unitary that acts on a qubit's neighborhood.The gate evolves different neighborhoods in different circuit layers.This layered evolution resembles the Trotterization of a many-body Hamiltonian for simulation on a digital quantum computer.
The unitary on site j's neighborhood has the form Vj := 1 ⊗j ⊗ V ⊗ 1 ⊗(L−1−j) denotes a unitary activation operator V that nontrivially evolves site j.The rule number's binary expansion determines the c i ∈ {0, 1}.The map between rule and {c i } is tedious, so we relegate it to Sec.IV B. The matrix power implies that V ci ∈ { 1, V }.If V = σx and the initial state is a product of |0 's and |1 's, the dynamics are classical.We therefore omit this V from this paper.
To define analog time evolution, we convert the rule number into a many-body Hamiltonian j Ĥj .The single-neighborhood Hamiltonian is defined as ĥj denotes a Hermitian activation operator that acts on site j.The rule number's binary expansion determines the c i ∈ {0, 1}.
To evolve the qubits under an analog QCA rule, we exponentiate the many-body Hamiltonian, forming Û = exp(−iδt Ĥ).This Û evolves the system for a time δt.Thus, evolving for one time unit requires 1/δt applications of Û .Choosing δt 1, we approximate Û with a product of local unitaries to order δt 4 [60] .
Analog and digital QCA tend to yield qualitatively similar outcomes.The relationship between analog and digital QCA is formalized in App.B.

C. Scope
Our digital QCA (1) and analog QCA (2) offer rich seams to explore: Many neighborhood sizes, timeevolution schemes, activation operators, boundary conditions, etc. are possible.One paper cannot report on all the combinations, so we balance breadth with depth.Below, we identify a manageable but representative subset of options.Appendix C demonstrates that conclusions of ours generalize to other initial conditions, boundary conditions, and activation operators.
We illustrate digital QCA with r = 1, or nearestneighbor, rules.This choice facilitates comparisons with ECA.We denote r = 1 digital QCA by T R .The T evokes three-site neighborhoods; and R = 0, 1, . . ., 15 denotes the rule number.We illustrate T R rules with a Hadamard activation operator, V = ÛH := (σ z + σx )/ √ 2, which interchanges the Bloch sphere's xand z-axes.We choose the Hadamard so that the QCA generate entanglement even when the lattice begins in a computational-basis state.During each time step, even-j sites are updated via the unitary (1), followed by odd-j sites [Fig.1(e)].We showcase these digital rules due to their conceptual simplicity, ease of implementation, and similarities with ECA.
We illustrate analog QCA with a time resolution of δt = 0.1 and r = 2, i.e., 5-site neighborhoods, [61].As 6.55 × 10 4 such rules exist, we must pare them down.We do so by focusing on totalistic rules, which update a site conditionally on the total number of neighbors in |1 .Only 32 totalistic r = 2 rules exist.Furthermore, totalism endows left-hand and right-hand neighbors with the same influence, enforcing a physically common inversion symmetry.We denote these rules by F R .The F refers to five-site neighborhoods; and R = 0, 1, . . ., 31 denotes the totalistic-rule number.We choose the activation ĥj = σx j , for three reasons.First, σx j serves as the Hamiltonian in many quantum simulations.Second, σx j transforms the computational basis, and so our typical initial states, nontrivially.Third, σx j can generate entanglement when serving as a Hamiltonian ĥj (as opposed to serving as a unitary Vj ).
In all simulations reported on in the main text, the boundary conditions are fixed: Consider evolving a boundary site, or any of the boundary's r neighbors.The evolution is dictated partially by imaginary sites, lying past the boundary, fixed to |0 's.Again, we exhibit other boundary conditions in App.C, supporting the generality of the main text's conclusions.
Throughout the main text, T R simulations are initialized with the product state |. . .000 1 000 . . . .This choice promotes simplicity; concreteness; and ease of comparison with ECA, which are often initialized similarly.F R simulations are initialized in |. . .000 101 000 . . .: One neighborhood must contain at least two |1 's for the Goldilocks F R rule to evolve the state nontrivially, as we will see below.Furthermore, the Goldilocks F R rule evolves this state into an entangled breather.This breather, first shown in Sec.III A and further analyzed in Sec.III D, parallels the breather of discrete nonlinear wave theory [62,63].

D. Entrée into Goldilocks and non-Goldilocks rules
To gain intuition about Goldilocks rules, consider the extreme digital rules T 0 and T 15 .T 0 always applies the identity operator to the target site.T 15 applies the activation operator to all sites at every time step.Both dynamics are trivial and generate no entanglement.As Goldilocks would say, T 0 is too passive, and T 15 is too active.
Interpolating between the extremes yields rules that increasingly approach what Goldilocks would deem "just right."Consider extending the do-nothing rule T 0 to updating sites whose two neighbors are in |0 .The resulting rule, T 1 , leads to nontrival dynamics.Similarly, consider constraining the update-everything rule T 15 : Each site is updated unless both its neighbors are in |0 's.Nontrivial dynamics emerge from the resulting rule, T 14 .
Optimizing the tradeoff between passivity and activity leads to the Goldilocks rule T 6 .Similar tradeoffs have been observed to generate and maintain complexity [21].Section III shows that tradeoffs generate complexity also in quantum dynamics and computation.
Having introduced the digital Goldilocks rule, we define the analog Goldilocks rule, F 4 .F 4 updates a site if exactly 2 of the site's 4 neighbors are in |1 .Let us construct the single-neighborhood Hamiltonian, Ĥj , for a Hermitian activation ĥj .We sum over the ways of distributing two P (1) 's across four neighbors: Having introduced the Goldilocks rules, we identify a near-Goldilocks rule that has been realized experimentally but whose status as a QCA has not been widely recognized.The PXP model is a Hamiltonian that has been implemented with Rydberg-excited neutral atoms [23] [Fig.1(d)].Such models generate diverse quantum dynamical features, such as many-body scars [64].Studies of PXP have spotlighted approaches conventional in many-body physics, such as energy-level statistics.Any physical complexity generated by the model, and the model's standing within the broader class of QCA, have not been explored.We quantify the complexity and identify the model as a QCA, although PXP is not a focus of this paper: As shown below, PXP (more precisely, a digital version of PXP) is a particular case of the near-Goldilocks rule T 1 .The model therefore generates less complexity than Goldilocks rules, our main focus.Additionally, rule T 14 has been explored theoretically under the name "rule 54" [65,66].We contextualize T 14 and T 1 , demonstrating that they are near-Goldilocks rules that belong to a broader family of QCA, which contains QCA that produce greater complexity.
To demonstrate the the power of this cellularautomaton-and-complexity-driven approach, we arrive at the PXP model not via conventional many-body considerations, but via ECA; and we observe rich many-body physics accessible to QCA beyond PXP: Figure 1(a) features the ECA rule C 201 , known to produce a simple periodic pattern called a blinker.The name derives from the oscillations of black and white squares that represent 1 and 0 bits.The closest digital QCA rule is T 1 , of which the digital PXP model is a particular case, which generates Figs.1(b)-(c).A quantum analog of a bit's being 1 or 0 is σz .T 1 produces not only visually richer dynamics than the ECA, but also quantum entropy.Progressing beyond PXP to other QCA reveals a truly nonclassical analog of a blinker, we show in the next section.

III. MAIN RESULTS: GOLDILOCKS QCA GENERATE COMPLEXITY
The prevalence of tradeoffs in complexity, with Goldilocks rules' regulation of a tradeoff, suggests that Goldilocks rules may generate complexity.This section supports this expectation with numerical evidence.As many readers may be unfamiliar with complexity science, Sec.III A initiates our numerical study of Goldilocks rules with measures familiar from many-body physics: onepoint correlation functions and von Neumann entropies.These metrics equilibrate under non-Goldilocks evolution but not under Goldilocks evolution.Building on this observation of Goldilocks rules' noteworthy dynamics, we calculate four complex-network measures [67][68][69][70] in Sec.III B. We evaluate these measures on networks formed from the mutual information between qubit pairs, in states generated by QCA.We complement these twopoint functions with many-body entropy fluctuations in Sec.III C. Additionally, the Goldilocks rule F 4 exhibits an emergent feature, the nonclassical analog of the ECA blinker.We study the structure and robustness of the entangled breather in Sec.III D. As explained in the introduction, complexity science is phenomenological.A system is regarded as complex if it exhibits several features suggestive of complexity [17].Together, the network measures, entropy fluctuations, and entangled breather form six signatures of complexity.We therefore conclude that Goldilocks QCA generate complexity.

A. One-point measures
As complexity science may be unfamiliar to many readers, our numerical study of Goldilocks rules begins with measures familiar from many-body physics.Figure 2 shows σz j (top row) and the single-site von Neumann entropy s j (bottom row) as functions of time.The entropy of a subsystem A in a reduced state ρ A is defined as s A := −Tr(ρ A log ρ A ).Further definitions and calculations appear in Apps.A 1 and A 2.
Under the QCA evolutions illustrated in the right-hand panels, σz j and s j equilibrate: After a time of the order of or much longer than the propagation time across the system, the plots become fairly uniform across space and show only small fluctuations.The Goldilocks rules, T 6 and F 4 , defy equilibration most dramatically [Figs.2(a We illustrate far-from-Goldilocks rules with F 26 , Goldilocks rules with T 6 and F 4 , and near-Goldilocks rules with T 1 and T 14 throughout the main text.Appendix C widens the survey to rule T 13 , a far-from-Goldilocks, nontotalistic digital rule.
One might wonder whether integrability or finite-size effects prevent the Goldilocks rules from equilibrating the measures in Fig. 2. Finite-size effects cannot take responsibility because the non-Goldilocks rules promote equilibration at the same system sizes used for the Goldilocks rules.Furthermore, App.C demonstrates our results' robustness with respect to changes in boundary conditions.Whether Goldilocks rules are integrable is discussed in Sec.VI but is irrelevant to the present paper, which focuses on complexity.The important point is that Goldilocks rules differ from other QCA even according to many-body measures.We now show that complexity distinguishes Goldilocks rules.

B. Complexity measures evaluated on mutual-information networks
We quantify QCA's complexity using complex networks.A network, or graph, is a collection of nodes, or vertices, connected by links, or edges.A complex network is neither regular nor random.Completely con- nected graphs (in which all nodes link to all others) are not complex; their links are regular and orderly.For the same reason, lattice graphs are not complex networks.Erdös-Rényi graphs is not complex for the opposite reason: Their links are created randomly.
Complex networks lie between the random and regular extremes.Concrete examples include social, airportconnection, neurological, and metabolic networks.Abstract examples include small-world networks, which exhibit clustering (nodes collect into tight-knit groups) and a short average path length (one node can be reached from another via few links, on average).
We study networks that emerge from quantum manybody systems.Networks can alternatively be imposed on quantum systems, such as networked NISQ (noisy intermediate-scale quantum) computers and a quantum Internet [71,72].Yet complex networks can also emerge naturally in quantum states generated by noncomplex systems or models.This emergence has been demonstrated with the Ising and Bose-Hubbard models near quantum critical points [67,73].We progress beyond critical phenomena to dynamical quantum systems and quantum computation.
Networks have been formed from the mutual information, as in brain activity that exhibits a small-world structure [24,25].The mutual information quantifies the correlation between two members of a system.We use the quantum mutual information [74,75] between qubits j and k = j, The 1/2, though nonstandard, conveniently normalizes the quantity as M jk ≤ 1.The mutual information depends on two-site von Neumann entropies that can be measured experimentally with quantum overlapping tomography [76].M jk upper-bounds every same-time, twopoint correlator between j and k [77].Appendix A 2 contains further details.We denote by M the matrix of elements M jk .Figure 3(a) displays mutual-information networks formed by QCA after 500 time units.Open circles represent the qubits, which serve as nodes.The greater the correlation M jk , the closer together nodes j and k lie, and the thicker the line between them.The Goldilocks rules T 6 and F 4 produce visibly clustered networks.The non-Goldilocks rules' networks more closely resemble the network encoded in the random state |R .|R is an L-qubit state of 2 L complex amplitudes.The amplitudes' real parts are drawn from independent, identi-cal Gaussian distributions, as are the imaginary parts.The near-Goldilocks rules T 1 and T 14 generate networks that interpolate between the Goldilocks and random networks.The far-from-Goldilocks rule F 26 produces a network strikingly similar to |R 's in shape and connectivity.Progressing beyond these qualitative features, we now quantify the networks' complexity with four metrics: clustering, average path length, disparity fluctuations, and the probability density over node strengths.

Clustering and average path length
The clustering coefficient C quantifies transitivity, easily intuited in the context of social networks: If Alice is friends with Bob and friends with Charlie, how likely is Bob to be friends with Charlie?The clustering coefficient is proportional to the number of connected triangles in the network, divided by the number of possible connected triangles: We will measure the clustering as a function of time [Fig.3(b)] and system size [Fig.3(d)].Figure 3(b) shows the clustering coefficient averaged over the network, C, as a function of time.Complex networks exhibit large clusterings, as discussed above.The Goldilocks rules (ruby and amber curves) maintain the highest clusterings.The nearest clustering characterizes the near-Goldilocks rule T 1 (turquoise curve).The T 6 clustering (ruby curve) exceeds the T 1 clustering by a factor of 10.5, and the F 4 clustering (amber curve) does by a factor of 2.2.The lowest clustering appears in the random network (black curve), generated as follows: The lattice is initialized in a random state |R , then evolved under one of the five QCA.The clustering is measured, then averaged over the five QCA.
The qualitative trends remain robust as the system size, L, grows: Figure 3(d) shows the clustering, averaged over the network and over late times t ∈ [5, 000, 10, 000] as a function of L. Only under the Goldilocks rules T 6 and F 4 (ruby and amber dots) does C grow with L. Under the other QCA, C decays exponentially, then plateaus.These behaviors approximate the behavior in the random state, wherein C strictly decays.Hence only the Goldilocks rules' complexity, as quantified with clustering, remains robust with respect to system-size growth.
As mentioned above, clustering combines with short average paths in small-worlds networks, which are complex.All our QCA simulations generated networks with short average path lengths ∼ 1.5 [78].Since Goldilocks rules generate high clustering and short average paths, Goldilocks rules generate small-world networks and so complexity.

Disparity fluctuations
Having measured the clustering and average path length, we turn to disparity.The disparity quantifies how heterogeneous a network is, or how much the network resembles a backbone [79]: Technically, Y is the disparity averaged over the network; the disparity is conventionally defined on one site.However, we call Y "the disparity" for conciseness.Two examples illuminate this definition.First, consider a uniform network: M jk = a is constant for all j = k.Along the diagonal, M jj = 0.The uniform network lacks a backbone and so has a small disparity: If L > 1 denotes the number of nodes, Y = 1/(L − 1) → 0 as L → ∞.In contrast, a 1D chain has an adjacency matrix M jk = a(δ j(k+1) + δ j(k−1) ).The chain consists entirely of a backbone, so the disparity is order-one for all system sizes L: Y = (L + 2)/(2L).Disparity has been applied in complex-network theory, e.g., to classify chemical reactions in bacterial metabolism [80].
Figure 3(c) shows fluctuations ∆Y in QCA disparity as a function of time.We define ∆Y as the standard deviation of Y in a rolling time window of L time units.∆Y reflects changes in a network's uniformity.The Goldilocks rule T 6 (ruby line) always produces disparity fluctuations above the pack's.The Goldilocks rule F 4 (amber line) always remains at the top of the pack, though the near-Goldilocks rules (T 1 and T 14 , represented with turquoise and yellow lines) do, too.The far-from-Goldilocks rule F 26 (green line), like the random state |R (black line), produce 1-2 orders of magnitude fewer disparity fluctuations.The |R data are generated analogously to in the clustering analysis.
Similar trends surface in Fig. 3(e), which shows disparity fluctuations averaged over space and over late times, ∆Y, as a function of system size.The ∆Y produced by the Goldilocks rules (ruby and amber dots) exceed the ∆Y produced by the other QCA, at all system sizes.Again, the near-Goldilocks rules (turquoise and yellow) mimic the Goldilocks rules most, and the far-from-Goldilocks rule (green) mimics the random state most (black).High disparity fluctuations imply a lack of equilibration in network structure: The network sometimes exhibits heterogeneity and sometimes exhibits uniformity.Such diversity signals the complexity of, e.g., Internet traffic [81].Here, the diversity signals the complexity of Goldilocks and, to some extent, near-Goldilocks rules.

Probability density over node strengths
Our fourth complex-network measure is the probability density over node strengths.The node strength, g j := Σ k M jk ≥ 0, measures how strongly node j is connected to the rest of the network.We draw this measure from analyses of social and airport-connection networks, real-world complex networks [82,83].These networks contain hubs: Few nodes carry a significant fraction of the network's total node strength.Hubs make the the node-strength distribution broad, e.g., power-law tails.
Figure 3(f) displays the node-strength densities of the networks generated by the QCA.The probability densities were calculated at late times t ∈ [5, 000, 10, 000].The Goldilocks rules (T 6 and F 4 , shown in ruby and amber) produce distributions that are significantly broader than the others: The T 6 points form a lopsided peak with a fat tail.The F 4 points form two peaks that ensure the distribution's width.Hence Goldilocks rules cause a few nodes to be strongly connected, as in known complex networks.The other rules, and the random state, lead to narrowly peaked distributions: Most connection strengths are centered on the mean.
We have investigated four complexity measures on mutual-information networks generated by QCA: clustering, average path length, disparity fluctuations, and the probability density over node strengths.
The Goldilocks rules tend to exhibit the most complexity; the near-Goldilocks rules, the second-most; the far-from-Goldilocks rule, less; and the random state, the least.

C. Bond entropy and its fluctuations
The previous section presented complexity measures evaluated on a network formed from the mutual information M jk .M jk depends on just two sites.We expand our study to larger subsystems, finding further evidence that Goldilocks rules generate complexity.Large subsystems have von Neumann entropies that have not been measured experimentally.However, the second-order Rényi entropy is accessible [84].This section therefore focuses on the second-order Rényi entropy.
The entropy is defined as follows.Consider partitioning the lattice at its center.The lattice's halves have equal entropies, associated with the bond between them and the bond entropy, s bond [85].If one lattice half occupies the quantum state ρ, s bond := − log 2 Tr(ρ 2 ) .To bridge this measure to complexity science, we consider fluctuations in s bond .Entropy fluctuations quantify complexity in biological systems [86].In quantum statistical mechanics [87], entropy fluctuations diverge at certain transitions between ordered and disordered phases [88].Complexity sits at the boundary between order and disorder, as explained in Sec.III B. For these reasons, we interpret high entropy fluctuations generated by QCA as signaling physical complexity.
Figure 4(a) shows s bond as a function of time.The curves generated by the Goldilocks rules T 6 and F 4 (the ruby and amber curves) fluctuate the most and persistently.Figure 4(b) shows the fluctuations ∆s bond plotted against time, for a 19-site system.The fluctuations produced by the Goldilocks rules remain at least an order of magnitude above the fluctuations produced by the non-Goldilocks rules.Entropy fluctuations therefore signal Goldilocks' rules complexity.
This behavior remains robust as the system grows.Let ∆s bond denote the entropy fluctuations averaged over late times.Figure 4(d) shows ∆s bond plotted against the system size.As L grows, the entropy fluctuations shrink most quickly for the random state (black dots) and the far-from-Goldilocks rule F 26 (green dots), less quickly for the near-Goldilocks rules (yellow and blue dots), and most slowly for the Goldilocks rules (ruby and amber dots).Hence Goldilocks rules generate complexity robustly, according to entropy fluctuations.
As an aside, although the Goldilocks rules produce entanglement, according to Figs. 4(a) and 4(c), they produce the least.This behavior might remind one of manybody scars-nonequilibrating features discovered under continuous-time T 1 evolution.Many-body scars mark a few area-law energy eigenstates states amidst a volumelaw sea.Yet observing many-body scars requires special initial conditions.Fine-tuning is not necessary for Goldilocks rules to hinder the entropy's equilibration: Complexity emerges for diverse initializations, activation operators, and boundary conditions, as shown in App. C.

D. Entangled breather
The sixth signal of Goldilocks rules' complexity is a robust emergent phenomenon.The emergent phenomena of life and consciousness signal complexity most famously.Other examples appear in economic markets [89], Internet traffic [90], physiological airways [91], and neural activity [92].Figures 2(c)-(d) display the entangled breather produced when the Goldilocks rule F 4 evolves the initial state |00 . . .0 101 00 . . .0 .Here, we analyze the breather's dynamics and robustness.
Breathers fall under the purview of discrete nonlinear wave theory.There, a breather is defined as a solution that oscillates in time and that is localized spatially [62,63].Let us establish that the phenomenon we observe satisfy these criteria.Imagine measuring the σz of each site j.The probability of obtaining the +1 outcome is P j plotted against time for key sites j.The central site is at j = L/2 , whose notation we simplify to L/2.The central site (j = L/2, purple curve) and its two nearest neighbors (j = L/2±1, amber curve) exhibit P (1) j 's that oscillate throughout the evolution.We have confirmed the first defining property of breathers.
To confirm the second property, we compare the oscillations near the lattice's center to the oscillations at the edges (j = 0, L − 1, red curve).The central oscillations have much greater amplitudes.We quantify this comparison with the fluctuations ∆ P , the temporalstandard deviation across a sliding window of L time units.Averaging the fluctuations over late times produces ∆ P , plotted against j in Fig. 5(b).The fluctuations peak strongly near the lattice's center.Hence the emergent phenomenon is strongly localized and so is a breather.
Two sets of evidence imply the breather's entanglement.One consists of the entanglement measure in Section III C. (Section III B displays results consistent with entanglement, although the mutual information captures classical correlations in addition to entanglement.)Second, the Schmidt-truncation study described later in this section implies the breather's entanglement.
This emergent phenomenon could be delicate; its lifetime could be exponentially suppressed by errors in, e.g., an experiment.The complexity generated by the Goldilocks rule F 4 would be tenuous.We demonstrate, however, that the entangled breather is robust with re- Entropy, and fluctuations in the entropy, of the qubit chain's central bond.Large entropy fluctuations are characteristic of macro-scale complex systems.Goldilocks rules display the largest such entropy entropy fluctuations of all the rules considered here.The black lines and dots, labeled |R , were generated as described in the caption of Fig. 3.The digital TR QCA and analog FR QCA were initialized as described there.(a) The central cut's second-order Rényi bond entropy, s bond , as a function of time for systems of size L=19.(b) The bond entropy averaged over late times t ∈ [5, 000, 10, 000], s bond , as a function of system size.(c) Entropy fluctuations ∆s bond , computed as a rolling standard deviation with a window size of L=19 time units.(d) Time-averaged fluctuations in the entropy, ∆s bond , as a function of system size.spect to perturbations of three types: perturbations of the Hamiltonian, of the entanglement, and of the initalization.Left unperturbed, the entangled breather oscillates for as long as our numerics run.To quantify the breather's robustness, we perturb the breather and measure the center site's neighbors (j = L/2 ± 1).We fit the fluctuations ∆ P (1) L/2±1 with the exponential A+Be −t/τ to extract breather's lifetime, τ .
The first perturbation is of the Hamiltonian.We scale the far-from-Goldilocks F 26 Hamiltonian with a positive positive ε 1 and add the product to the F 4 Hamiltonian.The perturbed dynamics produce an entangled breather whose oscillations decay.Figure 5(c) shows the fluctuations ∆ P (1) L/2±1 and their fits for various ε.The lifetime τ decays with ε as the power law τ ∼ ε −1.6 , as shown in Fig. 5(d).The power law decays more slowly than an exponential, indicating the entangled breather's robustness with respect to Hamiltonian perturbations.
Second, we perturb the QCA via Schmidt truncation: We evolve the system with a matrix-product-state approximation.We reduce the number of singular values used to represent the state, χ, from 2 L/2 = 2 8 until the breather becomes unstable.This truncation time-adaptively caps the amount of entanglement in the state [93].Under Schmidt truncation, the lifetime τ changes sharply at χ ∼ 11.Below the threshold, τ ∝ e 0.6χ .At the threshold, τ ∼ 340.Above, τ exceeds the simulation time of 6, 000 times units.Hence we regard the entangled breather as robust with respect to Schmidt truncation.
Third, we initialized the lattice imperfectly: The |1 's were replaced with copies of e iφ δ |0 + √ 1 − δ 2 |1 , wherein φ ∈ [0, 2π) and 0 < δ ≤ 1. Whenever δ ≤ 1/2, the breather's lifetime exceeds the simulation time, regardless of φ.When δ > 1/2, the breather's oscillations still persist indefinitely, but with a smaller amplitude.The entangled breather is therefore robust also with respect to a wide range of perturbations to the initialization.This emergent phenomenon, with the bond-entropy fluctuations and the four mutual-information measures, points to the Goldilocks rules' generation of complexity.

IV. MODEL DETAILS
In this section, we present formulae for the T R unitaries and the F R Hamiltonians, as functions of rule number.Our rule numbering is inspired by elementary-cellularautomata rule numbering, which we review first.

A. Classical-rule numbering
A classical-ECA neighborhood consists of a central bit and its two nearest neighbors.Classical ECA rules are denoted by C R .The rule number R encodes the local transition function, which prescribes a bit's next state for every possible neighborhood configuration.A bit has two possible values (0 or 1), and ECA neighborhoods consist of three sites.Hence a neighborhood occupies one of 2 3 = 8 configurations.Each configuration can lead to the central bit's being set to 0 or 1; so 2 2 3 = 256 C R rules exist.
The system is evolved through one time step as follows.First, a copy of the bit string is produced.Each site's neighborhood is read off the copy and inputted into the local transition function.The function's output dictates the central site's next state, which is written to the original bit string.The copying allows all sites to be updated independently, simultaneously.
The label C R encodes the rule as follows.Given C R , expand R into 8 binary digits, including leading zeroes: R = 7 n=0 c n 2 n .The index n is a base-10 number that can be translated into binary: n = 2 m=0 k m 2 k = k 2 k 1 k 0 .These three bits represent a possible initial configuration of the neighborhood.If the neighborhood begins in the n th possible configuration, the central site is updated to the c n ∈ {0, 1} in the binary expansion of R. Figure 6 shows the update table for rules R = 30, 60, and 110.Also depicted are the time evolutions of a 1 centered in 0s.

B. Digital quantum cellular automata
T R denotes a digital 3-site-neighborhood QCA rule.To expose the rule's dynamics, convert R from base-10 into four base-2 bits (including leading zeroes): R = c 11 2 3 + c 10 2 2 + c 01 2 1 + c 00 2 0 = 1 m,n=0 c mn 2 2m+n .The m specifies the state of a central site's left-hand neighbor, and n specifies the right-hand neighbor's state.c mn ∈ {0, 1} specifies whether the central bit evolves under a single-qubit unitary V (when c mn = 1) or under the identity 1 (when c mn = 0), given the neighbors' states.The three-site update unitary has the form The superscript on Vj denotes a matrix power.We update the lattice by applying Ûj (V ) to each 3-site neighborhood centered on site j, for every j.Even-j neighborhoods update first; and odd-j neighborhoods, second.One even-site update and one odd-site update, together, form a time step.

C. Analog quantum cellular automata
F R denotes a totalistic analog QCA rule that updates each site j, conditionally on four neighbors, with a Hermitian activation operator ĥj .The value of R dictates the form of the many-body Hamiltonian Ĥ = j Ĥj .We expand R into 5 binary digits: R = 4 q=0 c q 2 q .Site j evolves conditionally on its neighbors through Ĥj = σx N (q) j denotes the projector onto the q-totalistic subspace, in which exactly q neighbors are in |1 's: The sum runs over 4 q permutations of the neighboring |0 's and |1 's.We report on ĥj = σx j for reasons explained in Sec.II C.

V. PHYSICAL IMPLEMENTATION OF QUANTUM CELLULAR AUTOMATA
Having established the QCA's rule numbering, we prescribe a scheme for implementing QCA physically.We focus on analog QCA because they naturally suit our target platform, an analog quantum simulator.In this section, T R refers to an analog 3-site QCA.Appendix B details the relation between analog and digital T R protocols.
We begin with general results, before narrowing to an example platform.First, we show that the longitudinaltransverse quantum Ising model can reproduce QCA dynamics.Many experimental platforms can simulate this Ising model and so can, in principle, simulate QCA.Examples include superconducting qubits, trapped ions, and Rydberg atoms.
We illustrate our experimental scheme with neutral atoms trapped in an optical-tweezer lattice and excited to Rydberg states with principle quantum numbers 1. Excited atoms interact through a van der Waals coupling that induces a Rydberg blockade: Two close-together atoms cannot be excited simultaneously.We exploit this blockade to condition a site's evolution on the site's neighbors.Moreover, Rydberg atoms can be arranged in various geometries subject to short-and long-range interactions [94].We exploit these freedoms to engineer nearest-neighbor rules with a linear atom chain and next-nearest-neighbor rules with a ladder-like lattice [Fig.7 This section is organized as follows.In Sec.V A, we map analog QCA onto the Ising Hamiltonian.Section V B shows how to implement this Hamiltonian with Rydberg atoms.We simulated these Rydberg-atom QCA numerically, obtaining results presented in Sections V C and V D.
A. From QCA to the Ising model Analog QCA involve the Hermitian activation operator ĥj , which we have chosen to be σx j in this paper.According to Eq. ( 2), ĥj evolves site j if the neighboring sites satisfy the conditions imposed by the projectors.Such a state change can be modeled with the Ising Hamiltonian: The projectors are implemented with local energy penalties that shift undesired transitions far off-resonance.
Let h x denote the transverse field; and h z , the longitudinal field.The Ising Hamiltonian has the form The interaction strength J is uniform across the neighbors k ∈ Ω j and vanishes outside.If h x = 0, the resulting Hamiltonian, Ĥ0 , can be represented by a block-diagonal matrix relative to the σz product basis.Each block corresponds to one total-energy value.Consider switching on h x J ≈ h z .Some transitions between Ĥ0 eigenstates are resonant in energy, while others require energy jumps ∼ J, depending on the relationship between J and h z .An initial state in one Ĥ0 block will evolve within the block.By tuning h z and J, we can impose the resonance condition for a set transitions.Now, we show how to choose values for h z and J. Spin j evolves under the influence of the other spins, in an interacting global quantum system.However, we approximate j as evolving in a classical field.h x serves as the j th qubit's intrinsic gap.h z + Jm serves as the effective external field's frequency.The effective external field creates an energy penalty dependent on the magnetization, m = k∈Ωj σz k .The magnetization counts the |0 's in Ω j and subtracts off the number of |1 's.
Under each of several T R rules, the activation operates only if a certain number of neighbors occupy |1 's: m must have one particular value.Examples include T 1 , T 6 , and T 8 .The Ising model can implement these rules if the longitudinal field is h z = −mJ.The central spin's oscillations will be resonant with the effective external field at the required m value.If m does not have the required value, the central spin's oscillations are strongly off-resonant.This off-resonance effects the projectors in Eq. (2).The central spin's maximum probability of flipping, across a Rabi cycle, is h x / h 2 x + (h z + Jm) 2 .Rules that activate at multiple m values would require additional fields and resonances [59].Examples include rule T 14 , which activates when m = 0, 2.
B. From the Ising Model to the Rydberg-atom quantum simulator Rydberg atoms are natural platforms for experimental implementations of the Ising model [96] in the parameter regime described above.Indeed, the Rydberg blockade naturally realizes the energy penalty necessary for implementing a QCA.The blockade has been experimentally observed [23] and recently exploited to simulate lattice gauge theories [97][98][99].Each atom's ground state, |g , is coupled to a highly excited Rydberg state, |r , via the Rabi frequency Ω. Atoms j and k = j experience the isotropic van der Waals interaction V jk = C 6 /a 6 jk .The C 6 > 0 denotes an experimental parameter, and a jk > 0 denotes the interatomic separation.We use the notation V a = V jk for all |j − k| = 1 and V nn = V jk for all |j − k| = 2.
In the next two sections' numerical simulations, we include the next-highest-order terms, V ho , in the van der Waals interaction [dashed lines in Figs.7)(j)-(k)].These V ho , we prove, do not affect the results significantly.
Rydberg atoms evolve under the effective Hamiltonian [23] ĤRyd.= Ω 2 ∆ denotes a detuning in the natural Hamiltonian that is transformed and approximated to yield Eq. ( 12).The detuning affects nearly all the atoms uniformly: The QCA, recall, has boundary sites set to |0 's.We implement these boundary conditions by tailoring ∆ for the boundary sites and keeping ∆ constant elsewhere.We map the Rydberg Hamiltonian onto the Ising Hamiltonian as follows.First, the lattice geometry is chosen such that each site j has 2r nearest neighbors k ∈ Ω j .Second, we set Ω = 2h x and ∆ = 2(h z − 2rJ).Third, we set the interatomic distance a such that V a = 4J for the neighbors k ∈ Ω j .
In the Ising Hamiltonian, V j,k = 0 for all every site k outside Ω j .This condition is well-approximated by the Rydberg Hamiltonian (12), due to the van der Waals interaction's rapid decay with separation.We fix parameters to reflect the physics of 87 Rb atoms with the excited state |r = 70S 1/2 , for which C 6 = 863 GHz µm 6 .The interatomic distance a 5.4 µm, and Ω = 2 MHz V a = 36 MHz [23].

C. Simulating TR rules with the Rydberg Hamiltonian
Here, we specify the parameters needed to simulate rules T 6 , T 1 and T 8 with the Rydberg Hamiltonian.Consider a linear atom chain, with interatomic distance a, in the Rydberg-blockade regime [Fig.7(j)].Appropriately setting the detuning ∆ can enhance the oscillations undergone by a neighborhood's central spin, conditionally on the spin's neighbors.
The analog version of rule T 1 imposes the Hamiltonian Ĥ = j P (0) The Ising Hamiltonian implements the projectors when the resonance condition is imposed between |000 Ωj and |010 Ωj .In each of these states, two neighbors occupy |0 's; two neighbors point upward, so the magnetization is m = 2. Hence the resonance condition implies that h z = −2J, or ∆ = −2V a for the Rydberg Hamiltonian.Rule T 8 effects the Hamiltonian (13), except each P (0) is replaced with a P The resonance must be enhanced when the neighboring sites satisfy m = 0. Hence we set h z = 0 or, equivalently, ∆ = −V a in the Rydberg Hamiltonian.
Numerical simulations comparing analog QCA and their Rydberg-atom implementation for T 1 are shown in Fig. 7(a)-(c).Figure 7(d)-(f) shows the comparison for T 6 .The QCA simulations (in each figure triplet's leftmost column) agree qualitatively with the Rydberg simulations (depicted in the center column).Each figure triplet's rightmost column compares the QCA and Rydberg simulations on the same axes at key sites.The agreement exists despite the small, high-order interaction terms V ho ∼ V a /64 [Fig.7(j)].

D. Simulating FR rules with the Rydberg Hamiltonian
Simulating the F R rules with Rydberg atoms requires that excited next-nearest-neighbor atoms have the same interaction energy as nearest-neighbor atoms.One can meet this condition by forming a zigzag lattice, shown in Fig. 7(k).In this arrangement, the distance a jk for all |j − k| = 1 is equal to the distance a jk for all |j − k| = 2. Hence the interaction energies are equal: V nn = V a .
We discuss, as an example, the implementation of the Goldilocks rule F 4 .According to this rule, a neighborhood's central spin changes state if 2 of its 4 neighbors are in |1 's.Hence, as for T 6 , m = 0, leading to the resonance condition ∆ = −2V a .Figures 7(g)-(i) depict the numerical results, and Fig. 7(k) depicts the atomic arrangement.The plots confirm strong qualitative agreement between the QCA dynamics and the quantum-simulator dynamics.Due to the lattice geometry, higher-order interaction terms are larger for F 4 than for T 1 and T 6 : In the F 4 case, V ho ∼ V a /27.Consequently, the QCA dynamics diverge from the Rydberg dynamics more: The entangled breather remains in the Rydberg simulation, but with about half the oscillation frequency.We have checked that, if the higher-order interaction terms are suppressed, the agreement is significantly improved, supporting the mapping between the Rydberg Hamiltonian (12) and the QCA Hamiltonian.
The Rydberg-atom dynamics agree visually with the analog QCA.We quantify the agreement with the normalized absolute deviation | σz j QCA − σz j Ryd.|/2, averaged over sites j and times t ∈ [0, 20].The results are 3.9% for T 1 , 5.9% for T 6 , and 4.2% for F 4 -all errors of, or slightly above, order 1%.Thus, Rydberg atoms offer promise for quantum simulation of QCA, for the observation of complexity's emergence under Goldilocks rules, and for studying QCA's robustness with respect to physical-implementation details.

VI. CONCLUSIONS
We have bridged the fields of QCA and complexity science, discovering the emergence of physical complexity under Goldilocks rules.Goldilocks rules balance activity and inactivity.This tradeoff produces structured entangled quantum states.These states are persistently dynamic and neither uniform nor random, according to complexity measures.Complexity signatures of Goldilocks rules include high clustering; short average paths; high disparity fluctuations; broad node-strength distributions; persistent entropy fluctuations; and robust emergent dynamical features, entangled breathers.
Non-Goldilocks rules still display interesting dynamical features.For example, the PXP model, a particular case of rule T 1 , has received much attention.As a near-Goldilocks rule, T 1 produces some complexity.However, the confluence of many strong complexity signatures distinguishes the Goldilocks rules.
This work establishes several opportunities for future research.For example, whether Goldilocks rules are integrable merits study.This question is orthogonal to the question of whether Goldilocks rules generate complexity, the present paper's focus.The present study, however, offers hints.First, we searched for nontrivial quantities conserved by the Goldilocks rules.We found only one: T 6 with a Hadamard activation unitary conserves j σz j σz j+1 , whose sum includes the boundaries, be they fixed or periodic.Second, we have performed preliminary studies on QCA's entanglement spectra [100,101], which suggest that Goldilocks rules might lie between integrability and nonintegrability [102].Many more measures have been leveraged to diagnose integrability and chaos [103][104][105][106]; these measures merit scrutiny.Furthermore, we observed Goldilocks QCA hinder equilibration (Sec.III).How does this resistance relate to mechanisms for avoiding thermalization in quantum systems, e.g., many-body localization [107], quantum scars [108], and Hilbert space fragmentation [106,109]?Contrariwise, we observed equilibration under non-Goldilocks QCA.Do they obey the eigenstate thermalization hypothesis [110]?
Finally, we have focused on finite-size systems accessible to classical computing, which exist in the physical world and can be simulated numerically.The thermodynamic limit forms an alternative, idealized setting.Do the complexity features observed in our finite systems persist, contract, or grow richer in the infinite-system-size limit?Preliminary indications based on finite-size scaling in Fig. 3(d) show that at least one of our complexity measures stays large for Goldilocks rules in the thermodynamic limit, namely clustering.Disparity and entropy fluctuations presents a borderline cases that cannot be resolved on classical computers based on our present studies, as their decay is slow and may in fact be leveling off to a plateau.Likewise, given that quantum entangled breathers require about 6-7 sites including their supporting tails, finite-size scaling based on multiple breathers as fundamental emergent objects is inaccessible on classical resources.For example, to scale to 10 breathers requires 60-70 qubits.Thus, understanding complexity in the thermodynamic limit is overall a good case for usage of quantum computers [30], and presents an excellent  Blue circles represent the atoms; red lines, the interactions with strength Va; and gray, dashed lines, the next-nearest-neighbor interactions (under TR rules) and next-next-nearest neighbor interactions (under FR rules), which generate higher-order effects.A local detuning is applied at the boundaries to simulate the QCA boundary conditions.Parameters: L=17.In the Rydberg simulations, Ω = 2 MHz, and Va = 36 MHz.Time is expressed in units of 2π µs.
topic for future research in NISQ-era computing.
Aside from these theoretical opportunities, we have demonstrated that extant digital and analog quantum computers can implement our QCA.For example, Rydberg atoms can be arranged to engineer five-site QCA neighborhoods, large enough to support an entangled breather.can simulate QCA dynamics that support complexity even with three-site neighborhoods.Our work therefore uncovers a direction for experimental quantum computation: to explore the rich features of biological and social complexity that manifest in simple, abiological quantum systems.funding from the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1125565) with support of the Gordon and Betty Moore Foundation (GBMF-2644), for a Barbara Groce Graduate Fellowship, and for an NSF grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory.

Appendix C: Supplemental Material
In this appendix, we provide additional analysis of several digital T R rules to support our claim that Goldilocks rules produce physical complexity for a variety of activation operators, initial states, and boundary conditions.Additionally, we provide an alternative visualization of the data in Fig. 2. We illustrate the class of non-Goldilocks rules with T 13 , a rule not covered in the main text, for the sake of exhibiting a wider variety of QCA.T 13 updates a site j if (i) both neighbors are in |0 's, (ii) both neighbors are in |1 's, or (iii) the right-hand neighbor is in |1 , while the left-hand neighbor is in |0 .This rule is nontotalistic.
For a more general activation operator, we consider the product of a Hadamard and a phase gate: V (υ) = ÛH Ûphase (υ).In the matrix representation with |0 = (1, 0) T and |1 = (0, 1) The main text details the υ = 0 case.Figure 8 shows key complexity metrics: The mutualinformation-network clustering [panel (a)], disparity fluctuations [panel (b)], and bond-entropy fluctuations [panel (c)] are observed under various rules and phase-gate angles υ.The long-time averages are calculated over the interval t ∈ [500, 1, 000].The Goldilocks rule T 6 exhibits the largest clustering, disparity fluctuations, and entropy fluctuations for all υ.Like T 6 , the near-Goldilocks rule T 1 is, in comparison to the other rules, less sensitive to nonzero υ.Conversely, the near-Goldilocks rule T 14 and the far-from-Goldilocks rule T 13 produce the most variation with υ and the lowest overall values at late times.
Let us next consider the complexity outcomes for alternative initial and boundary conditions (Fig. 9).We initialize each QCA in a random product state where each qubit is set to |1 with probability p and otherwise set to |0 .The QCA evolution is performed; and the reduced density matrices of single sites, two-site pairs, and the central bipartition are recorded at each time step.The process is repeated 500 times with different realizations of the initial condition.At each time step, the reduced density matrices are averaged across the trials.From the ensemble-averaged reduced density matrices, we compute the long-time-average clustering [row (a)], disparity fluctuations [row (b)], and entropy fluctuations [row (c)].We repeat these simulations for boundaries fixed to |0 's (left column), boundaries fixed to |1 's (center column), and periodic boundary conditions (right column).
Boundary conditions play a relatively insignificant role in determining complexity outcomes, as compared to initial conditions.Also, the clustering and and entropy fluctuations are near their minimum, for all rules, with initial conditions near p = 50%.Disparity fluctuations, in comparison, appear less sensitive to changes in the initial state.The Goldilocks rule T 6 generally maintains the highest clustering, disparity fluctuations, and entropy fluctuations for all p.The exception is near p = 50%, where the near-Goldilocks rule T 1 exhibits a greater clustering.The far-from-Goldilocks rule T 13 appears the least complex for all initial conditions.The near-Goldilocks rules T 1 and T 14 lie somewhere between T 6 and T 13 .
Finally, Figure 10 shows the one-point dynamics as a function of time for key sites.This visualization provides an alternative to Fig. 2 that allows one to compare more directly the site-local observables for various QCA rules.
FIG. 1.(a) The classical-elementary-cellular-automaton rule C201 flips a classical bit (black=1, white=0) if and only if its nearest neighbors are 0s [4].This L=19-site chain is initialized with a 1 centered in 0s.A periodic pattern, called a blinker, propagates upward in discrete time steps.(b) We extend C201, for QCA, to rule T1, defined in the main text.The classical bit becomes a quantum bit in an L=19-site chain.The average spin σz j has richer discrete-time dynamics, oscillating between the classical extremes of 0 (black) and 1 (white) and filling the lattice.A truly quantum analog of the classical blinker-an entangled breather-appears in Figs.2(c)-(d).(c) The quantum lattice evolves into a high-entropy entangled state.sj denotes the site-j von Neumann entropy.(d) An analog quantum computer has implemented rule T1, of which the recently named PXP model is a particular case.Optical tweezers were used to trap a Rydberg-atom chain [23].(e) A quantum circuit can realize digital T1 dynamics.One QCA time step requires two layers of controlled-controlled-Hadamard gates.The first layer evolves even-index qubits; and the second layer, odd-index.The dashed line represents the boundary qubits, which remain fixed in |0 's.
)-(b) and (c)-(d)].Their measures display persist patterns.For instance, F 4 generates an entangled breather analyzed in Sec.III D. The near-Goldilocks rule T 1 [Figs.2(e)-(f)] resists equilibration second-best: Although its patterns fade, they persist across time.The fading is stronger under the near-Goldilocks rule T 14 [Figs.2(g)-(h)].Rule F 26 [Fig.2(i) and (j)] encourages equilibration the most: σz j quickly equilibrates to zero (rose color), and s j quickly equilibrates to yellow (a high entropy value) across the lattice.F 26 is a far-from-Goldilocks rule: It evolves site j if the site has one, three, or four neighbors in |1 .
FIG. 2. QCA produce diverse dynamical behaviors, visualized here with familiar one-point observables.For (a)-(b) and (e)-(h), the digital TR rules, an L=19-site chain was initialized to |. . .000 1 000 . . ., then evolved in discrete time under various TR rules.The Goldilocks rule T6 generates the spin dynamics shown in (a) and the entropy dynamics shown in (b).The dynamics exhibit patterns that resist equilibration for long times, as shown in each panel's upper block.For (c)-(d) and (i)-(j), we extend QCA to 5-site rules, initialize the L=19-site chain to |. . .000 101 000 . . ., and evolve in continuous time.(c)-(d) The Goldilocks rule F4 generates a nonclassical variation on the blinker in Fig. 1(a), an entangled breather.(i)-(j) In contrast, far-from-Goldilocks rule F26 promotes rapid equilibration to a high-entropy state.

FIG. 3 .
FIG.3.Mutual information networks distinguish Goldilocks rules as generating the most complexity.(a) An L=19-site spin chain was initialized, then evolved for 500 time units.In each of three cases, the initialization was to |00 . . .0 1 00 . . .0 , and a TR QCA dictated the evolution.In two other cases, the initialization was to |00 . . .0 101 00 . . .0 , and an FR QCA dictated the evolution.The open circles represent the qubits.Nodes j and k are connected by a line whose thickness is proportional to the mutual information M jk .The nodes' layout is found by modeling links as springs; well-connected nodes lie close together.conditions are chosen to match those of Fig.2.The networks' labels provide a legend for the panels below.The black lines and dots show the results of initializing the qubits to a random state |R , evolving the state under one of the five QCA, then averaging the final state over the QCA.(b) Clustering coefficient C vs. time t.(c) Disparity fluctuations ∆Y vs. t.(d)-(e) Clustering coefficient C and disparity fluctuations ∆Y, averaged over late times t ∈ [5, 000, 10, 000] (denoted by the overline), vs. system size L. (f) Probability density over node strengths gj collected at t ∈ [5, 000, 10, 000].Statistics were accumulated with logarithmically spaced bins.Only bins with > 75 counts were maintained.The inset shows the same data on a linear horizontal scale and identical vertical scale.
FIG. 4.Entropy, and fluctuations in the entropy, of the qubit chain's central bond.Large entropy fluctuations are characteristic of macro-scale complex systems.Goldilocks rules display the largest such entropy entropy fluctuations of all the rules considered here.The black lines and dots, labeled |R , were generated as described in the caption of Fig.3.The digital TR QCA and analog FR QCA were initialized as described there.(a) The central cut's second-order Rényi bond entropy, s bond , as a function of time for systems of size L=19.(b) The bond entropy averaged over late times t ∈ [5, 000, 10, 000], s bond , as a function of system size.(c) Entropy fluctuations ∆s bond , computed as a rolling standard deviation with a window size of L=19 time units.(d) Time-averaged fluctuations in the entropy, ∆s bond , as a function of system size.

FIG. 5 .
FIG.5.The entangled breather produced by Goldilocks rule F4 has a characteristic structure and is robust to perturbations.Here we analyze the breather on an L=19-site lattice and denote the center site by L/2.(a) The probability of site j being in |1 , P(1)j , shows temporal oscillations occur with the greatest amplitude near the lattice's center.(b) Late-time average of fluctuations in the probability of sites being in |1 , ∆ P (1) j , quantifies the spatial localization of temporal oscillations.The coloring is as in panel (a).(c) Upon perturbation by the F26 Hamiltonian, scaled by ε, the characteristic fluctuations near the lattice's center decay exponentially.The black lines show exponential fits that yield a lifetime τ .The decaying curves are shifted by 0.05 from one another to aid visualization.The dashed horizontal lines mark the shifted x-axes.(d) The entangled breather lifetime decays as as a function of the Hamiltonian-perturbation strength according to a power law.The coloring is as in panel (c).
(k)].Related work appeared in the literature recently: A Rydberg-chain realization of QCA was proposed in [59] [95].Wintermantel et al. focus on 3-site neighborhoods.We present a more-general version of unitary Rydbergatom QCA, which encompasses 3-site and 5-site neighborhoods.

−1 1 FIG. 7 .
FIG. 7. (a)-(i) Comparison of QCA evolution and Rydberg-atom dynamics.We study σz j for the near-Goldilocks three-site rule T1 in (a)-(c), the three-site Goldilocks rule T6 in (d)-(f) and the five-site Goldilocks rule F4 in (g)-(i).Panels (c), (f), and (i) show the results for key sites j, distinguished by line style.Colored lines reflect the exact QCA results; black lines represent their Rydberg-atom simulations.(j)-(k) Lattice geometries for the three-site TR rules (j) and five-site FR rules (k).Blue circles represent the atoms; red lines, the interactions with strength Va; and gray, dashed lines, the next-nearest-neighbor interactions (under TR rules) and next-next-nearest neighbor interactions (under FR rules), which generate higher-order effects.A local detuning is applied at the boundaries to simulate the QCA boundary conditions.Parameters: L=17.In the Rydberg simulations, Ω = 2 MHz, and Va = 36 MHz.Time is expressed in units of 2π µs.

FIG. 8 .
FIG. 8.We generalize the activation unitary from a Hadamard gate to a Hadamard gate times an azimuthal phase gate of angle υ.The long-time trends that distinguish rule T6 remain.(a) Clustering, (b) disparity fluctuations, and (c) bond entropy fluctuations at averaged over times t ∈ 500, 1, 000.As in the main text, these simulations are performed with L = 19 qubits, an initial condition of a single |1 centered in |0 's, and boundary qubits fixed to |0 's.

FIG. 9 .FIG. 10 .
FIG.9.Goldilocks rule T6 generates the most complexity even under various boundary conditions and when initial states are generalized to ensembles.The long-time averages of (a) clustering, (b) disparity fluctuations, and (c) bond-entropy fluctuations are shown for various TR rules and excitation probabilities p.The results are shown for boundaries fixed to |0 's (left columns), boundaries fixed to |1 's (center column), and periodic boundaries (right column).For these simulations, the system size is reduced to L = 15, and the long-time average is taken for t ∈ (250, 500], to mitigate computational expense.