Machine learning discovery of new phases in programmable quantum simulator snapshots

Machine learning has recently emerged as a promising approach for studying complex phenomena characterized by rich datasets. In particular, data-centric approaches lend to the possibility of automatically discovering structures in experimental datasets that manual inspection may miss. Here, we introduce an interpretable unsupervised-supervised hybrid machine learning approach, the hybrid-correlation convolutional neural network (Hybrid-CCNN), and apply it to experimental data generated using a programmable quantum simulator based on Rydberg atom arrays. Specifically, we apply Hybrid-CCNN to analyze new quantum phases on square lattices with programmable interactions. The initial unsupervised dimensionality reduction and clustering stage first reveals five distinct quantum phase regions. In a second supervised stage, we refine these phase boundaries and characterize each phase by training fully interpretable CCNNs and extracting the relevant correlations for each phase. The characteristic spatial weightings and snippets of correlations specifically recognized in each phase capture quantum fluctuations in the striated phase and identify two previously undetected phases, the rhombic and boundary-ordered phases. These observations demonstrate that a combination of programmable quantum simulators with machine learning can be used as a powerful tool for detailed exploration of correlated quantum states of matter.


I. INTRODUCTION
Recent advances in realization of programmable quantum simulators (PQSs) have opened up a new era in the exploration of strongly correlated quantum matter [1][2][3][4], which calls for new approaches for analyzing large volumes of data generated by such quantum devices.Using optical techniques, it is possible to arrange a large number of qubits in arbitrary lattice geometries [5] and to control the Hamiltonian evolution of the system [2] dynamically in real time.Remarkably, these simulators can probe states within an extremely large Hilbert space.For example, in a 13 × 13 system, the quantum states live in a 2 169 -dimensional space while each measurement probabilistically projects to just a single dimension.Tomographically [6][7][8] inferring the entire manybody wavefunction from such measurements themselves is a formidable task.While certain types of many-body states can be easily identified by evaluating simple local observables, many exotic quantum phases that can be explored on programmable simulators cannot be characterized using conventional approaches.
In this work, we introduce a hybrid unsupervisedsupervised machine learning approach to analyze the * eun-ah.kim@cornell.edudata generated using programmable quantum simulators.Specifically, we apply this method to a PQS based on Rydberg atoms arrayed on a square lattice [5].We show how this framework directly reveals the correct order parameters to diagnose multiple density-wave-ordered phases, thus providing a route to the construction of (a priori unkown) order parameters for potentially more complicated symmetry-breaking phases [9].Additionally, our analysis allows us to uncover several new features including: (i ) the pattern of quantum correlations in the "striated phase"; (ii ) a previously undetected "rhombic" phase in consistency with the theoretical predictions of Ref. 10; and (iii ) a boundary-ordering quantum phase transition in which the edge of the system develops long-range order while the bulk remains trivial.
Our programmable Rydberg quantum simulator [see Fig. 1(a)] consists of neutral atoms trapped in defect-free arrays of optical tweezers with programmable geometries [5].Coherent laser excitations to atomic Rydberg states realize an Ising-like spin model [2,4] described by the Hamiltonian where atoms in the ground (Rydberg) state are denoted by |g (|r ), and n i ≡ |r i r i |.The transverse field Ω corresponds to the Rabi frequency of the laser field, the First, an unsupervised technique is used to generate a rough first-pass phase diagram.Here, we choose to measure average Fourier amplitudes |n(k)| 2 at each (∆, R b ), perform a dimensionality reduction using principal component analysis, and finally cluster using a Gaussian mixture model.The resulting phase diagram informs the starting "seeds" in the parameter space, from which snapshots are sampled in a second supervised stage.We then learn to distinguish these snapshots using interpretable classifiers, from which we can extract refined phase boundaries and key identifying features.
longitudinal field ∆ corresponds to the laser detuning, and 6 is the long-range van der Waals interactions between Rydberg excitations at x i and x j .
Density-matrix renormalization group (DMRG) calculations on the square lattice [10] have predicted a number of quantum phases of the Hamiltonian (1) arising from the interplay between coherent laser driving and the long-range van der Waals interactions [see Fig. 1(be)].These phases can be understood based on the Rydberg blockade phenomenon [11]: the strong interactions V ij can prohibit (or "blockade") the simultaneous excitation of neighboring atoms to the Rydberg state.The spatial extent of this blockade (or equivalently, the interaction strength) is captured by the blockade radius, defined as R b ≡ (V 0 /Ω) 1/6 .The full phase diagram is thus parametrized by the ratio of the longitudinal to the transverse field, ∆/Ω, and R b /a, where a is the lattice spacing.For ∆/Ω > 0, the system energetically favors maximizing the number of atoms in the Rydberg state.However, this is subject to the blockade constraint, so for R b /a 1, only one out of every pair of nearest neighbors can be excited; on a square lattice, this leads to the checkerboard phase with antiferromagnetic ordering of atoms in ground and Rydberg states.Higher values of R b /a result in various new density-wave-ordered phases.Some of these correspond to classical hard-sphere packing of Rydberg excitations [10], while others are stabilized by quantum coherence between the ground and Rydberg states [5,10].
Recent experiments [5] have demonstrated three of these predicted states [Fig.1(b-d)], namely, the checkerboard, striated, and star phases.In the experiments, different values of R b /a are accessed by tuning the lattice spacing a at fixed V 0 .By linearly ramping ∆/Ω from negative to positive values, one can quasiadiabatically prepare different ordered states, which are probed by measuring order parameters derived from the Fourier transforms of the excitation densities.However, some of the more complex ordered phases are challenging to detect and analyze directly, especially in light of the finite system size and experimental imperfections.In particular, the finite system sizes introduce significant edge effects that impede the identification of phases based on manual inspection of the Fourier transforms.
To address the challenge associated with large volumes of data produced by PQS projective measurements, we introduce an unsupervised-supervised hybrid machine learning approach: Hybrid-CCNN.In Ref. 12, some of us introduced the CCNN architecture, which modifies the standard convolutional network to allow direct interpretation of the key-feature correlations of a learned phase.While this interpretability endowed the CCNN with the ability to reveal new theoretical insights, a direct supervised learning approach is ultimately limited by the subjective bias entering through training data selection.Increasingly, the community has made efforts to develop "phase recognition" [13] algorithms to curtail such bias either by using fully unsupervised learning [14][15][16][17][18][19][20][21][22][23][24] or by adopting an element of unsupervised learning [25][26][27].However, discovery of new phases or fluctuation effects from experimental data through such efforts are rare [19] to date.Here, we introduce "Hybrid-CCNN" that prefaces CCNN's supervised learning with an unsupervised discovery stage and apply the two-stage process shown in Fig. 1(f) on voluminous Rydberg PQS data to arrive at three discoveries.
An innovative feature of the Hybrid-CCNN is the use of different data preprocessing for the two stages to investigate different facets of the high-dimensional data: overall density modulations and snapshot-to-snapshot fluctuations.At each point in the (∆, R b ) phase space, the available snapshot data consists of 250 binary maps For the unsupervised learning stage, we focus on density modulations starting from the mean density-normalized snapshots ñi (x) = n i (x) − n, where n is the average density across all sites and snapshots.Then, Fourier transforming to |ñ i (k)| 2 , subtracting the background contribution ∼ k |ñ i (k)| 2 , and averaging over all snapshots leads us to a density-shift-invariant structure factor (refer to Appendix A for details).Beginning our analysis in Fourier space allows us to benchmark the Hybrid-CCNN findings against results obtained using a traditional measure of spatially modulated order for finite-size simulations (see Appendix B).
On the other hand, much information about the structure of many-body quantum fluctuations captured in each snapshot will be lost upon this averaging.More-over, the Fourier basis is noisy in the presence of nonperiodic spatial variations and incommensurate domains.Hence, in the final supervised phase, we choose to train our CCNNs to learn phase characteristics directly from the full density fluctuation maps δn i (x) ≡ n i (x) − n(x) produced from each snapshot as independent training data.Here, n(x) = i n i (x)/N is the snapshot-averaged Rydberg excitation density at each site, where N = 250 is the number of snapshots at each (∆, R b ).As detailed later, this normalization allows the CCNN to easily focus on fully connected components of low-order correlation functions and prevents the network from learning trivial overall excitation densities.

II. UNSUPERVISED PHASE DISCOVERY
For the unsupervised learning stage, the aim is to form a collection of observables and use a clustering algorithm to map out different regions in the (∆, R b ) space where the observable set changes dramatically.The challenge is that clustering tends to be most robust in lowdimensional feature spaces.On the other hand, even after averaging over snapshots, the full Fourier-space dimension is too large for standard clustering algorithms.Specifically, working with a 16 × 16 grid in k-space [28] implies a 256-dimensional feature vector associated with each point in (∆, R b ) space.To make this manageable for clustering, we further reduce the dimensionality of the feature space using principal component analysis (PCA) [29].
PCA identifies vectors called principal components (PCs) in this 256-dimensional feature space along which the data varies the most dramatically across the full phase diagram.Specifically, each of the principal components is a linear combination of multiple points in k-space that vary together across the phase space, visualized in Fig. 2(a-c).PC1 and PC2 show considerable overlap with the theoretical order parameters for checkerboard and star phases, respectively.While PC3 and PC4 do not immediately offer interpretation as an order parameter in Fourier space, possibly capturing peak broadening, PC5 resembles the Fourier-space order parameter associated with the rhombic phase predicted in DMRG simulations on a cylinder [10].This is tantalizing since the previous manual analysis of the data only detected three phases: checkerboard, star, and striated [5].
The final step of the unsupervised learning stage is to cluster the phase-space points (∆, R b ) in the reduced feature space spanned by the first few dominant PCA components.As our clustering algorithm, we use a Gaussian Mixture Model [29] (GMM) for its robustness and invariance to the scale of each feature.The two choices we must make are the number of principal components to keep for the clustering, and the number of clusters to fit.We choose the first by increasing the number of retained principal components one-by-one until the clusters stabilize, finding 10 to be sufficient in the process.We then determine the optimal number of clusters to be six as the Bayesian information criterion [29] plateaus past this number (see Appendix A).
The clustering result shown in Fig. 1(f) partially resembles the manually obtained phase diagram based on evaluation of three target order parameters for the checkerboard, striated, and star phases [5] (see Appendix C), except that the unsupervised learning indicates that there are five phases distinct from the disordered phase [shown in grey in Fig. 1(f)].A visualization of the clusters projected to the subspace spanned by PC1, PC2, and PC5 shows that the clusters are tightly defined [Fig.2(d)].Although the clusters reside in submanifolds unrestricted to a single principal component axis, the red, green, and purple clusters are located far along the PC1, PC2, and PC5 axes, respectively.Hence, the unsupervised learning stage of Hybrid-CCNN identifies the checkerboard phase (red region and PC1) and the star phase (green region and PC2) with no prior knowledge.Moreover, it discovers evidence of a new phase associated with a large PC5 component (purple).At this stage, the identities of the orange and blue phases also remain unclear from the unsupervised learning alone.

III. SUPERVISED PHASE CHARACTERIZATION A. Architecture and training
To better characterize each of the phases from the full dataset of snapshot-to-snapshot fluctuations, δn i (x; ∆, R b ), we now turn to a supervised learning stage.The unsupervised learning results indicate five phases distinct from the disordered phase and inform the starting choice of the training points.However, since the CCNN has access to the full information about quantum fluctuations, we anticipate changes to the precise phase boundaries.Hence, we optimize the choice of training points by iteratively exploring to minimize training loss, overlap between different phases, and understandibility [30].To learn the distinct identity of each phase, we train multiple neural networks, with each given the task of identifying snapshots of a single phase against the rest through a binary classification [31].For the neural network architecture, we adapt the CCNN introduced in Ref. 12 to learn the spatial structure of cor-relations specific to each phase [see Fig. 2(e)].CCNNs recognize m-site correlations [12] from site-normalized density fluctuation snapshots δn i (x) through m th order polynomials of convolutional maps produced with learned filters f α (a), C α (x) allows the CCNN to discover higher-order and spatially inhomogeneous correlations that are often overlooked, but which can be crucial for identification of a given phase [12].For the present dataset, we find truncating to m = 3 to be sufficient.Specifically, we have For interpretability, we fix the filters to be positive definite.In addition, we learn spatial structures beyond the length scale of the filters using a follow-up spatial weighting w(x) [32] applied to the correlation maps C The sign and magnitude of the learned weights β (m) α can then be investigated to reveal the ferromagnetic (+) or antiferromagnetic (−) nature and significance of the identified correlations.This whole process is pictorially summarized in Fig. 2(e).
For training, the snapshots from the target training points are labeled with y = 1 and those from the remaining training points are labeled with y = 0.During training, the filters f α (a), logistic weights β (n) α , bias , and spatial weighting w(x) are all simultaneously learned by stochastic gradient descent to minimize the cross-entropy loss, which drives the predicted ŷs towards their correct labels y (Appendix A).In order to arrive at a minimally complex set of features describing each phase, we perform ablation studies [33] (described in Appendix A) to progressively determine the minimal number of components necessary for the CCNN's successful training.The resulting nuggets of each phase allows us to connect the full snapshot dataset to theoretical insight.

B. Interpreting the learned phases
We first focus on the red and green phase-space regions in Fig. 2(f).We find that for these phases, a fixed uniform spatial weighting (w(x) = 1) and truncation to second-order correlations are sufficient for characterization.In Appendix C, we show by the convolution theorem that a linear combination of uniform c (2) α features amounts to a sum of the structure factor |δn(k)| 2 weighted by the Fourier transform of the filters f α (a).Inspecting the learned effective weightings in Fourier space (shown in Appendix C), we remarkably find that these CCNNs identify the correct order parameters traditionally used to characterize the respective density-wave orderings.Indeed, comparing these two regions to the earlier result [5] based on manual evaluation of order parameters in Fourier space, the red and green regions of Fig. 2(f) clearly map to checkerboard [Fig.1(b)] and star [Fig.1(d)] phases, respectively.For these phases, the key advantages of CCNN-based phase recognition are an unbiased discovery of the simplest order parameter suitable for the complexity of fluctuations present in the data (see Appendix B), as well as noticeably sharper phase boundaries.We emphasize the nontriviality of the former observation as the Hybrid-CCNN is able to reveal the correct order parameters without any prior input about the physics or structure of the density-wave-ordered ground states.This highlights the utility of our method for applications to potentially more complicated symmetry-  breaking phases, for which the correct order parameters may not be immediately obvious.

Fluctuations in the striated phase
Next, we turn to the blue region in Fig. 2(f).From the CCNN trained to identify this phase, we produce the confidence map shown in Fig. 3(a) by "classifying" each point in parameter space by averaging the correlator features c (m) α across all available snapshots and then predicting using Eq. 5. We find that resulting map overlaps significantly with the region previously demarcated as the striated phase in Ref. 5 using an approximate Fourier order parameter [Fig.3(b)].However, the CCNN confidence map produces a remarkably sharper region of support for the phase.In the parameter regime of interest, the striated phase is unique in that it is an intrinsically quantum phase stabilized by the quantum fluctuations driven by the transverse field ∼ Ω; in fact, the Rydberg interaction term alone energetically favors star order, so the striated phase does not exist in the classical limit of Ω/∆ → 0 [10].Motivated by this, previous work [5] characterized this phase using a mean-field ansatz consisting of a product of single-particle states in quantum superpositions of |g and |r .However, such an ansatz describes a product state.In contrast, the two and three-site connected correlations that the CCNN learned to recognize using two filters f 1 and f 2 [Fig.3(c)], combined with the long-ranged structure learned by the weight map w(x) [Fig.3(d)], offer a first glimpse at potentially nontrivial quantum many-body correlations in this region.
First, the weight map w(x) learned to identify a specific sublattice in the bulk, and correspondingly activates the filters only when they are centered on this sublattice.As a result, the CCNN measures correlations within repeating 3 × 3 blocks that span the system [see Fig 3(e)].The learned positive sign of the logistic weight β (2) 1 associated with the filter f 1 implies that the corner sites will tend to be jointly excited, which matches the expected pattern of excitations in the striated phase.Moreover, the weight map identifies correlations within this effectively larger unit cell via the filter f 2 .Specifically, the learned negative logistic weight β (2) 2 , combined with the pattern in filter f 2 , indicates that the CCNN identified pairwise anticorrelations between three sites in the unit cell as a relevant feature of the striated phase.Note that such pairwise anticorrelation of three sites is not possible for any perfectly polarized product state, in which the atom at each site is exclusively in either |g or |r .Rather, Fig. 3(c The learned spatial weighting and filters for the model trained to identify the edge-ordered phase, with the spatial extent of the snapshot indicated by the dashed line.The outermost pixels correspond to where the filter is centered on zero-padding but "clip" the edge pixels of the snapshot.For display purposes, the filter weights are normalized such that the maximum is 1 within each filter.(c) Measured performance discrepancy between second-order models with a fixed uniform w(x) = 1 and a freely learned spatial weighting.The central lines and bands show the mean and standard deviation across five randomly initialized models of each type, respectively.The persistence of nonzero second-order connected correlations (which should ideally vanish deep in an ordered phase) calls for a description beyond a simple productstate ansatz.Quantum correlations and entanglement can arise from two primary sources.Firstly, they might be present in the ground state itself, particularly in the vicinity of a second-order quantum phase transition [34].Second, they might be generated in the dynamical state preparation process due to the quantum Kibble-Zurek mechanism [35], where nonadiabatic processes can coherently generate superpositions including excited states that generically result in entanglement.Our Hybrid-CCNN approach cannot discern between these scenarios as it is agnostic to the actual origin of the correlations.Nonetheless, to gain further insight into the potential entanglement structure, we inspect the bipartite entanglement entropy and correlations within a 9 × 9 system using the density-matrix renormalization group (DMRG) in Appendix D. We find that the calculated von Neumann entanglement entropy S peaks sharply along transition lines, before plateauing to a small but nonzero value within the phase.This is accompanied by anticorrelations between the excited sublattices as found by the CCNN, though the state preparation process appears to significantly extend the support of these correlations as compared to the DMRG ground state (see Appendix D).

The boundary-ordered phase
Having analyzed the properties of the striated phase, next, we turn to the first of the two mysterious phases, depicted in dark orange in Fig. 4(a).The first clue regarding the identity of this phase comes from the learned weight map w(x), which focuses strongly on the edges of the snapshots.As shown in Fig. 4(b), the CCNN learned to measure the differences in correlations between the bulk and the boundary by having large w(x) > 0 along the edge and predominantly w(x) < 0 in the interior.The learned filters focus the CCNN's attention on specific short-range two-point correlations that differ significantly between the edge and bulk.Figure 4(c) demonstrates the dramatic performance gain enabled by the learned edge-centered weight map.Inspection of the experimental snapshots in this orange phase [Fig.4(d)] indeed confirms a regular occurrence of local Z 2 patterns of (• • •) along the edges of the snapshots.In contrast, the bulk of the snapshots appear disordered, further evidenced by explicit measurements of correlation functions along the edge and in the bulk in Fig. 4(e).
While the local Z 2 pattern is commensurate with the neighboring checkerboard and striated phases, the reduced energetic cost of having Rydberg excitations along the boundary (relative to the bulk), due to the presence of fewer neighbors, actually allows the edge to order before the bulk.Hence, we identify this mystery phase as a boundary-ordered phase characterized by the edge ordering in the absence of long-range order in the bulk.The  to be nearly zero, so we omit two-point correlators stemming from the filter f2.Our CCNN is symmetrized, (see Appendix B) and so measures all correlators symmetry-equivalent under rotations and flips to those shown.(j) Identification of these two-and three-point motifs in the idealized rhombic ordering with boundary defects (light blue).(k) Identification of local occurrences of these motifs in experimental snapshots sampled from the training set.
subsequent onset of bulk order, in the presence of preexisting edge order, defines an "extraordinary" boundary universality class [36].We highlight that the existence of this boundary-ordered phase is a new discovery of the present work since this phase cannot be obtained on geometries with fully periodic or cylindrical boundary conditions, as was used for earlier DMRG calculations [10].Critical to the identification of this phase is the real-space nature of the CCNN analysis as the edge ordering introduces a large number of artifacts into |n(k)| 2 , which can challenge traditional Fourier-based analysis.Interestingly, a complementary work [37] independently detected this edge ordering in quantum Monte Carlo simulations of the system with open boundary conditions, and confirmed the first-order nature of several transitions.

The rhombic phase
Finally, we examine the other mystery phase identified by the Hybrid-CCNN: the purple swath in Fig. 5(a).To identify the defining characteristic of this phase, we restrict the CCNN to learn positive correlation functions by enforcing β (n) α ≥ 0 during training, increase the filter size to 4 × 4, and fix uniform w(x) = 1 (see Appendices C and D).The CCNN learns the two filters f 1 and f 2 shown in Fig. 5(b,c) and uses a combination of second-and third-order correlations c α to recognize this phase as shown in Fig. 5(d-i).Remarkably, these learned motifs strongly point towards the rhombic phase from among the candidate ordering patterns in Fig. 1(b-e).
The rhombic phase is an intricately patterned densitywave-ordered phase characterized by Fourier peaks at ±(π, π/4), ±(2π/5, π) (and their C 4 -rotated copies) [see Fig. 1(e)], which was originally predicted by Ref. 10.However, given its large unit cell comprising 40 sites, the robustness of this phase in the actual experimental system of Ebadi et al. [5] is a priori unclear due to both the long-ranged tails of the van der Waals interactions and the incompatibility of the ideal ordering pattern with the dimensions of the lattices used.Our results illustrate that, interestingly, we can still find characteristic remnants of this phase.In particular, the three-point motif of Fig. 5(d) provides a unique signature of the rhombic phase as a fragment of a full rhombic crystal while Fig. 5(g) occurs as an edge defect when the rhombic pattern is embedded in the finite incommensurate system as shown in Fig. 5(j).Additionally, the shorter-range threepoint motifs of Fig. 5(h) and (i) occur most frequently in the rhombic phase (see Appendix E).The virtue of these motifs is that they signify the tendency of fluctuations towards rhombic ordering even when extended ordered portions cannot form inside a finite system (due to the large and incommensurate 6 × 5 unit cell).Indeed, these motifs are ubiquitous in the experimental snapshots sampled from this phase region as we showcase in Fig. 5(k).
Hence, we can unambiguously identify the second mystery phase to be the finite-size manifestation of the rhombic phase.

IV. DISCUSSION
In summary, we developed a supervised-unsupervised hybrid machine learning approach, the Hybrid-CCNN, to reveal collective quantum phenomena in voluminous quantum snapshot datasets and applied our approach to square-lattice Rydberg PQS data.The initial unsupervised stage used Fourier intensities |n(k)| 2 and clustering in a low-dimensional feature space obtained using PCA to reveal a rough initial phase diagram.This first pass reveals the number of phases to expect and informs the initial location of training points for the supervised stage.The phase diagram is then refined in the second supervised stage by training nonuniform CCNNs to define sharper phase boundaries and uncover the identities of each phase.The identities revealed in this way not only confirmed the previously known phases [5,10] but also resulted in new insight into potential quantum entanglement structures in the striated phase and the discovery of two previously undetected phases: the edge-ordered phase and the rhombic phase.To the authors' knowledge, this is the first time new insights and discoveries were gained from quantum simulator data using machine learning tools trained entirely on experimental data.
With the rapid progress towards probing more exotic quantum many-body phenomena using PQS [38][39][40], the need for new data-centric approaches to extracting insight from large volumes of quantum snapshot data will only grow.Here, we demonstrated that the Hybrid-CCNN can not only meet this need but also enable the discovery and identification of new quantum states.The Hybrid-CCNN's ability to extract spatial structures of a quantum state at multiple length scales is particularly valuable given the limited spatial extent and incommensurate domains of phases produced by finite systems under quasiadiabatic state preparation [35].In particular, this ability enabled the discovery of the edge-ordered phase and rhombic phases in our present work.Given the challenges to establishing and detecting new ordered phases due to the inevitable breakdown of adiabaticity when sweeping across a quantum phase transition [35], it will be interesting to compare the motifs learned by Hybrid-CCNN with rigorous finite-size scaling studies of corresponding quantum Monte Carlo simulations.Also tantalizing is the Hybrid-CCNN's ability to detect nonclassical correlations (that defy description in terms of product states) in the snapshot data of pure states; for instance, this hinted at nontrivial quantum entanglement among diagonal bonds in the striated phase.Comparing CCNN-detected entanglement against traditional measures is an interesting direction for future exploration.
We anticipate such probes of entanglement to critically contribute to the identification and understanding of new exotic states that are beginning to become experimentally accessible [40].The authors acknowledge helpful discussions with Soonwon Choi, Marcin Kalinowski, and Roger Melko.We would also like to thank H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein for the use of experimental data presented in this work.The DMRG calculations in this paper were performed using the ITensor package [41] and were run on the FASRC Odyssey cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

Appendix A: Unsupervised learning details
To perform our initial rough-estimate unsupervised phase discovery, we first collect features at each experimental parameter point (∆, R b ).For easy interpretability, we choose to work with simple features that represent all second-order fluctuations, but are blind to the overall Rydberg excitation density present in the snapshots.We find that blinding our entire machine learning pipeline to the overall density results in phase boundaries that are more closely aligned to meaningful changes in spatial orderings rather than simple increases in Rydberg excitation densities.More complex phase transitions would require higher-order correlations or nonlinear unsupervised learning techniques, but we find the following approach to be sufficient for our data.To allow for direct comparison to previous experimental work [5], we would like to work with the average squared Fourier amplitudes of the snapshots at each (∆, R b ).In order to make the process blind to overall density changes, we must take two steps.
First, in this section, we work with snapshots which are overall density-normalized as δn i (x) ≡ n i (x) − n , where the expectation value is computed over all sites of all snapshots available at the same (∆, R b ).Note that this is different than the per-site density normalization used in the supervised follow-up.In the supervised phase, we found that per-site density normalization was necessary to build interpretable order parameters for phases with subtle correlation structures such as the striated and rhombic phases, which could be masked by average density modulations induced by the edge ordering.Per-site density normalization subtracts out this average density modulation, allowing connected correlators to be easily measured.However, while difficult to interpret, these very same average density modulations appear to be key to the success of our unsupervised phase, which is restricted to measuring raw Fourier amplitudes.
We then measure the average squared Fourier amplitudes of these normalized snapshots: Our first density normalization does not make δn(k) invariant under an overall shift in the density of the n(x)s, even away from k = 0, due to the input data being sampled from a bosonic system with the values of n(x) restricted to 0 or 1.In this case, the marginal distribution of each individual site's density is a Bernoulli distribution whose variance is linked to its mean as δn(x) 2 = n(x) (1 − n(x) ).In this way, information about the density can "bleed" through to the learning process at all k-points.We can see this explicitly by expanding out the square in Eq. (A1) and rewriting it as To make our Fourier-space features invariant under density shifts we need to remove the first term.Using Plancherel's theorem, we can achieve this this by normalizing once more in k-space as with L being the number of k-points sampled along each direction of the discrete Fourier transform.These resulting features p(k; ∆, R b ) are then provided to the principal component analysis, the results of which are summarized in Fig. A6. Figure A6 shows the resulting top 12 principal components, which together capture > 99.9% of the variance of the dataset.Notably, since all of the input features lie within the k p(k) = 0 subspace, so do the resulting principal component vectors.
Below each PCA vector in Fig. A6, we show its average inner product with each p(k; ∆, R b ) across parameter space.These maps reveal which areas of parameter space have average Fourier-space intensity patterns that match the principal component vectors well.Importantly, many of these maps vary smoothly in parameter space and as such, correspond to meaningful changes in Fourier structure.Meanwhile, at first glance PCA6, and partially PCA10, seem visually noisy in parameter space.However, this is due to these PCA components breaking a rotational symmetry in Fourier space.In phase regions where the relevant Fourier peaks are present, these PCA components are either strongly postive or strongly negative depending on which way the symmetry is broken, while in other regions they remain close to zero.These components could be improved by further postprocessing; however, for simplicity, we do not do so here.
To robustly perform the GMM clustering, we initialize the clusters using the k-Means algorithm [29], repeat the initialization and clustering with different random seeds until 500 sequential clusterings do not improve the final log-likelihood, and keep the best-found clustering.Interestingly, random initializations (rather than k-Means) often produce clusterings with higher final log-likelihoods but poor structures in parameter space, with one cluster often a seemingly-random collection of points across parameter space.This ambiguity points to the persisting need to have a physicist "in the loop" verifying machine learning (ML) results at each stage, which is made easier when using interpretable techniques.
To determine the appropriate number of clusters and PCA components to retain, we perform the entire clustering process while varying the number of clusters and the number of retained PCA components.In Fig. A7(a), we find that starting from two clusters, each additional cluster meaningfully captures the next-strongest phase in the dataset, with topology similar to the clustering of the main text with N clust = 6.As expected, each additional cluster improves the final log-likelihood [Fig.A7(b)]; however, a distinct kink and reduction in slope is observed after N clust = 6.This is reflected additionally by the Bayesian information criterion (BIC) [29], a standard heuristic metric for determining the correct number of clusters, plateauing past this point.As we should pick the simplest model which minimizes the BIC, this points to N clust = 6 as being optimal.Clusterings with N clust > 7 are increasingly noisy and can be seen to simply be chunking off transition regions between phases, rather than dramatically changing the phase diagram's topology.
In Fig. A8, we show clustering results with a fixed N clust = 6 but varying number of PCA components kept.Across all settings, we predominately see the topology presented in the main text (similar to N PCA = 9), with exceptions being N PCA = 5, 6, 7 which substitute the rhombic phase with a cluster around the broadenedcheckerboard region characterized by PCA3.The relative stability and understandability of these clusterings give us confidence that they indicate real regions of varying orderings suitable for a follow-up analysis.As truncating the number of PCA components kept is primarily a costsaving and robustness-improving technique, we increase the number of principal components kept until the clustering is seen to stabilize, which can be seen in Fig. A8 to be roughly 10.We implement our Convolutional Correlator Neural Networks using the Pytorch [42] library, and our code is made available at github.com/KimGroup/QGasML.Each CCNN is trained to identify a single phase: all snapshots sampled within that phase are labeled 0, while snapshots sampled from all other phases are labeled 1.The list of parameter-space points sampled to form each phase's training set is given in Table I.From each point, we randomly take 90% of the snapshots as the training dataset shown to the network, and take the remaining 10% as a validation dataset which is not shown during training, and is only used to verify that the network has not overfit.To handle the uneven distribution of snapshots available from each phase, when training a CCNN to identify a phase P each presented snapshot has a 50% probability to be sampled from phase P, and a 10% probability to be sampled from each of the 5 other phases.This ensures that there is an equal representation of "within-phase-P" snapshots and "out-of-phase-P" snapshots, as well as equal representation among all classes within the "outof-phase" distribution.
The training points in Table I are obtained by starting with points suggested by the unsupervised clustering, and modifying iteratively based on the results of training until all phases overlap minimally and have clear distinctions.For the striated phase, whose support was very narrow in the unsupervised learning results, we explored a wider range in phase space to identify the training points yielding physically meaningful features.
During training, the free parameters of our CCNNs, {f α (a), β (m) α , , w(x)}, are all simultaneously learned so as to minimize the training loss averaged across the dataset, chosen to be the standard cross entropy loss used for classification tasks with an additional L1 regularization on the filter weights: and N is the total number of snapshots.We optimize this loss using the ADAM [43] optimization algorithm, with a minibatch size of 128 snapshots, an initial learning rate of 0.01, and a cosine-annealed learning rate schedule as implemented by Pytorch's CosineAnnealingLR [42].
As in Ref. 12, we also place a BatchNorm [44] layer (without the optional affine transformation) after the correlator features c (m) α , which introduces no extra free parameters but aids in rapid and stable convergence of the training.Similar to Ref. 12, we additionally spatially symmetrize C (m) α (x) and w(x) to improve generalization and interpretability.Specifically, on each forward pass, we symmetrize the maps by summing over all symmetry transformations as with the sum running over all group elements g ∈ D 4 acting on the convolutional maps, and the same transformation additionally applied to w(x).This symmetrization procedure aids in generalization, as the model's predictions are made invariant under symmetry transformations of the input snapshots, and reduces the effective total number of parameters to learn.In particular, this allows the model to avoid having to learn symmetryequivalent versions of the filters f α (a), and makes the spatial weighting w(x) easier to visually interpret.To make interpretation simpler, we also restrict f α (a) to take only positive values by applying an absolute value function on every forward pass.
To ensure that all pixels of the learned filters hit all pixels of the snapshot in the convolution, we zero-pad each δn i (x) with a sufficient number of zeros.In particular, if the convolutional filters f α (a) are of spatial extent F × F , we pad the input snapshots with F − 1 zeros on all edges.
After training is completed, a CCNN has learned a collection of convolutional filters f α (a), as well as a set of coefficients β (m) α connecting the m th order feature c (m) α derived from filter f α to the output, and an overall bias .The output of the model is then where c (m) α are constructed from the learned filters f α and the input by Eqs. ( 2) and (3) of the main text.

Ablation testing
When building a specific CCNN, there is a lot of flexibility in the architectural choices.These architectural hyperparameters include the order m to truncate at, whether to include a learnable spatial weighting w(x), TABLE II.Final validation accuracies for models identifying snapshots of each phase using various architectures.Measurements are made with 10-fold cross-validation [29], with 5 random seeds run at each train/validation split.Reported errors for each model type are the standard error across all runs.All models are trained with 3 convolutional filters.Figure 2(f) of the main text is produced by the final, most expressive, row of models.Table entries in italics are simpler models studied in Figs.4,5 of the main text.Models sampled from the "2nd Order, w(x) = 1" row produce the order parameters shown in Fig. A10 of Appendix C.
and how many filters (and of what size) to use.For the current work, we take the approach of building the simplest architecture (second-order, uniform w(x)) first, and adding on architectural complexity piece-by-piece.If a large gain in accuracy is achieved by a single architectural addition, then we keep that piece and attribute a quality of the phase to requiring the additional expressibility.This approach is commonly referred to as ablation testing in the ML literature [33], though for large neural networks, ML practitioners commonly remove modules piece-by-piece rather than adding them as we do for our shallow networks.
To improve interpretability, an additional hyperparameter that we have at our disposal is the coefficient γ on the filter L1 loss in Eq.A1.Intuitively, larger γ results in simpler filters with more pixels deactivated but decreased classification performance.For all models except for the uniform second-order models (which we can easily interpret in Fourier space, see Appendix C), we increase γ until the filters identifying all phases are sparse enough to easily interpret while performance is maintained sufficiently high.
Along these lines, in Table II, we summarize validation accuracy measurements of several variations of our CCNNs trained to identify each phase.From these measurements, in conjunction with observing the quality of the resulting phase diagram, we can determine what is required to form a good order parameter for each of the identified phases.For example, the checkerboard, star, and rhombic phases show roughly uniform or even decreasing (due to overfitting) validation performance as the CCNN's truncation order is increased, indicating that second-order features are sufficient to distinguish these from other phases.However, we find the second-order features identifying the rhombic phase to be somewhat uninteresting as they primarily just measure the tendency for longer-range density correlations; see Appendix C. The star phase shows some improvement when incorporating spatial inhomogeneity, but we find that the phase diagram changes little while the second order model is simplest to interpret (see Appendix C).
The most striking changes are observed for the two re-maining phases, both of which benefit heavily from learning a nonuniform spatial weighting w(x).For the striated phase, we observe a dramatic 11.6% jump in classification accuracy between uniform and nonuniform second-order models.This difference reflects that due to finite-size effects, many striated-like correlations persist throughout the system within the star phase.Appendix C presents an alternate interpretation in Fourier space, where the difficulty originates from the relevant Fourier peaks being overly diffuse.However, many of these striated-like correlations in the star phase occur on the even sublattice, while the striated ordering entirely occupies the odd sublattice.If the network has the ability to focus on a specific sublattice matching the ideal striated ordering (making all other sites contribute with a negative weight to the output), it can better distinguish the striated from the star phase.Alternatively, we find that increasing the order of the model while keeping w(x) = 1 also results in a good classifier for the striated phase but one which is more difficult to understand.
In contrast, for the edge-ordered phase, we find that simply increasing the order of the model while keeping uniform w(x) = 1 makes negative changes to generalization performance.Meanwhile, keeping the model at second order, we see a 4.9% jump in accuracy when allowing for a spatially varying w(x).While this may not seem dramatic, we find that all models with uniform w(x) produce order parameters which persist deep into the disordered region.This indicates that successfully identifying this phase requires measuring spatially-inhomogeneous correlation functions, which reflects that this phase is itself defined by this inhomogeneity.
Appendix C: Extracting Fourier-space order parameters from uniform second-order CCNNs Second-order CCNN features (under a uniform spatial weighting w(x) = 1) have an alternate interpretation as performing weighted sums of the Fourier spectrum of the input.Specifically, each feature c (2) α can be written as where Eq. ( A1) is just the expanded definition of C (2) (x), Eq. (A2) is the Fourier-equivalent form, and discrete Fourier transforms are defined by where L f , L n are the lengths of each dimension of the convolutional filter f and the input snapshot δn, respectively.The wavevectors in the discrete Fourier transform are defined as k i ∈ {0, 2π/N, 4π/N, . . ., (N − 1)2π/N } with N ≥ L defining the resolution of the Fourier transform.Choosing N > L encodes no new information in the result, but produces higher-resolution Fourier spectra for plotting.
Note that for all terms in Eq. (A1) to be well-defined, we take an "infinite zero-padding" convention, where the domain of δn is expanded to all x by defining δn(x) = 0 Under this convention, we take the x sum to be over all space.In practice, our CC-NNs equivalently pad δn with enough zeros to recover all nonzero C (2) α (x).Equation (A2) can be obtained from Eq. (A1) by straightforward application of the convolution theorem and Plancherel's theorem.
A simple interpretation of this result is that a uniformweighted c (2) α measures a weighted sum of δn(k), with the weights given by the Fourier transform of f α , normalized to be zero-mean in k-space.Since second-order CCNNs linearly combine multiple filters f α with coefficients β α to produce the input to the final logistic function, we can understand the full action of the network by the weighted sum of the input in k-space with effective weights f α fα (k), with normalized Fourier intensities of each filter defined as fα (k Due to the symmetrization procedure outlined in Appendix B, this map must be then symmetrized over all symmetries of D 4 .This "order parameter map" defines our effective second-order CCNN order parameter for a given phase as with σ(x) = (1 + exp(−x)) −1 the logistic sigmoid function, and the learned bias.
In Fig. A9(a), we provide a visual diagram demonstrating this process to produce an interpretable f (k) sym .In Fig. A9(b-d), we show exemplary resulting Fourier weighting maps learned for each of the checkerboard, striated, and star phases, and demonstrate the intuition to interpret these maps.Applying the weighting f (k) sym to the average Fourier intensity of the target phase should produce a large positive number, while applying it to the intensities of all other phases should produce smaller or negative numbers (after subtracting the learned , only the target phase should remain positive).
We find that the checkerboard, striated, and star Fourier weightings look strikingly similar to smearedout versions of the hand-crafted order parameters discussed in Refs. 5 and 10, which can be interpreted in our framework as f (k) that are nonzero only at finitely many k-points.In particular, our learned checkerboard order parameter is strongly positive at (π, π), the striated at (±π, 0), (0, ±π), and the star order parameter spreads positive weight towards (π/2, 0), (0, π/2) while putting negative weight at (π, π).This observation is reassuring as it demonstrates that these models identify each phase by an ordering similar to the ideal density waves.
In Fig. A10(a-c), we show measurements of Fourierspace order parameters manually crafted to capture each of the checkerboard, star, and striated phases.We can see that due to the blurring of the relevant Fourier peaks resulting from the finite-size system, the order parameters for the star and striated phases [Fig.A10(b,c)] heavily overlap.For comparison, in Fig. A10(d-h), we show confidence maps from uniform second-order CCNNs with 4 × 4 convolutional filters trained to recognize each phase examined in the main text.We can observe that only the checkerboard (red) and star (green) phases are wellresolved by these reduced models, in accordance with the intuition that these phases are classical crystals easily identified in Fourier space.Meanwhile, just as with manually crafted Fourier-space order parameters, the uniform second-order CCNN's predictions for the striated phase significantly overlap with those of the star phase.These combined observations point to a fundamental limit to phase resolution using Fourier-space order parameters in finite, open-boundary, noisy systems.The rhombic (purple) phase is somewhat well-resolved, but it too suffers from overly broad peaks in Fourier space.Meanwhile, we see that uniform second-order models fail dramatically α , then symmetrized to form the final map.(b-d) The order parameter maps f (k) sym learned to identify the checkboard, striated, and star phases.We apply each order parameter map as weightings to the idealized Fourier intensities for each of the checkerboard, striated, and star phases.If learning is successful, applying this weighting and then summing in k-space should produce large positive numbers for the target phase, and small or negative numbers for every other phase.
for the edge-ordered (orange) phase, which by its nature fundamentally requires the ability to learn spatially inhomogeneous functions, as measured in Appendix B.
To uncover what is being learned by the rhombic second-order CCNN, we repeat the above analysis to produce the Fourier-space order parameter shown in Fig. A11(a).We see that the learned order parameter attempts to identify Fourier intensity along the diagonals connecting the (±2π/5, 0) and (0, ±π/2) peaks (and symmetry equivalents), as the intensity uniquely appears here only in this phase.Due to broadening resulting from experimental noise and the finite-size system, the CCNN does not attempt to directly measure the rhombic (±2π/5, 0) peaks as these blur too strongly into the star phase's (±π/2, 0) peaks.Nevertheless, these peaks can be visually resolved when averaging over a large number of snapshots as in Fig. A11(c).
Manual inspection of the learned filters and β α coefficients reveals that the rhombic CCNN obtains this order parameter by placing negative β (2) α on filters which contain short-range patterns.This makes intuitive sense, as density fluctuations between nearby sites are anticorrelated at high R b .However, this does not give us clear insight into the actual Rydberg crystal being realized, other than that it is constructed from longer-range displacements.This motivated our choice in the main text of investigating third-order models which enforce β (m) α > 0 to ensure that the rhombic CCNN must rely on positive correlations to identify the phase.Due to three-point correlators having a nontrivial sign structure (see Appendix E), the CCNN still measures some short-range patterns, but also learns several key longer-range patterns uniquely characterizing the phase.As often exploited by Rydberg-atom quantum simulators, the Rydberg blockade effect can generate entanglement between interacting atoms [45][46][47].Here, we perform simple calculations to examine the interplay between entanglement and correlations in the vicinity of the transition to the striated phase and draw connections to the correlations in the experimental data uncovered by our CCNN analysis.To this end, we perform densitymatrix renormalization group (DMRG) calculations using a "snakelike" matrix product state ansatz for a 9 × 9 system with open boundaries.The Rydberg Hamiltonian realized by the experiment examined in this work is given by Eq. ( 1), and for notational simplicity, we will refer to the first, second, and third terms therein as Ω, ∆, and V , respectively.Entanglement can be generated in the ground state of Ĥ(≡ H/ ) due to energetic competition between these terms.For any two sites, the interaction term V prefers small overlap with basis states containing |rr , while ∆ desires large overlap with all of |gr , |rg , and |rr .Crucially, Ω favors weight to be present with opposite phase between basis states with a single site flipped as |g ↔ |r .As a result, for a two-site system, the ground state as R b /x 12 → ∞ places weight across all of the |gg , |rg , |gr basis states, but not |rr , resulting in an entangled, anticorrelated state.
To examine this behavior more closely, we now turn to the results of the DMRG computations.In Fig. A12, we show the bipartite entanglement entropy between two halves of the system.Within the disordered phase (∆ 1), the entanglement entropy of the ground state increases monotonically as one approaches the quantum critical points, and, at large R b , the density on any site is anticorrelated with that of its next-nearest neighbors in the corners [Fig.A12(d)].As we transition deep into the classically ordered phases, both the entanglement and the connected correlations vanish due to the density on each site approaching either 0 or 1 [Fig.A12(b,c)].However, in a narrow region at R b /a ≈ 1.4 where quantum fluctuations stabilize a significant density on the (1, 1) sublattice [see Fig. A12(d)], both entanglement and diagonal anticorrelations survive.We emphasize that this entanglement is dependent on the state on each sublattice remaining in a quantum superposition of |g and |r , as is uniquely characteristic of the striated phase.
Indeed, as shown in Fig. A13, in the experiment, the connected part of many short-range correlations remains finite and negative upon transitioning into the striated phase.In particular, all nearest-neighbor correlators remain anticorrelated, along with next-nearest-neighbor correlators between the two excited sublattices.However, there are many confounding effects which make it difficult to pinpoint with certainty the origin of these correlations.First, due to decoherence and experimental noise, the experiment, in principle, produces a mixed rather than a pure state.Given a mixed state ρ, nonzero connected correlations can emerge either by entanglement, or by "classical" correlations between the pure states comprising ρ [48].Efficient means for unambiguously revealing entanglement in experimental settings without full tomography is an active field of research [49][50][51][52][53][54][55].
Moreover, quasiadiabatic sweeps across phase boundaries produce final states which are not perfect ground states but are dependent on the original state starting from which the phase boundary was crossed.As the system is theoretically transitioning from a region of high entanglement, where nearby neighbors have anticorrelated densities, it is likely that some of the magnitude of correlations captured by the CCNN is not directly representative of the true ground state, instead having been "frozen in" from before the transition [35].Nevertheless, the nature (and, in particular, the signs) of these correlations still reveals qualitative structures of each identified phase.I), performed on each individual nearest-neighbor (NN), next-nearest-neighbor (NNN), or next-to-next-nearest-neighbor (NNNN) bond, respectively.The colors of each site are for visual aid, showing the sites expected to be in the mostly excited (pink), mostly ground (purple), and entirely ground (black) states in the ideal striated limit.(d, e) Averaged correlations across all bonds of different symmetry classes, tracked as a function of ∆/Ω for a cut at R b = 1.56.

Appendix E: Sign structures of third-order correlators
By inspecting the learned β α and the patterns in the learned associated filters f α , we can determine which three-point correlations are being measured for a target phase, and whether they should be positive or negative within the phase.However, the sign on a three-point contribution can be somewhat confusing to interpret, as there are multiple ways to obtain positive/negative three-point correlations.The CCNN itself does not inherently point out how to interpret these correlatorssimply that they are strongly positive/negative within the phase.Manual follow-up and investigation is always necessary to understand the underlying physics.This section attempts to clarify the subtleties of these measured three-point correlators, and presents explicit measurements thereof from the data to confirm that the CCNN's identification was meaningful.
Within this section, for notational brevity, spatial indices are written as (vector) subscripts.Given densitynormalized snapshots δn i ≡ n i − n i , uniform thirdorder CCNN features measure weighted sums of connected three-point correlation functions, averaged over all spatial translations: ijk=(0,0) f i f j f k δn x+i δn x+j δn x+k , (A1) where the inner sum runs over all configurations of displacements (i, j, k) within the spatial extent of the filter f of length L f .Due to the symmetrization process outlined in Appendix B 1, this must also be averaged over all rotations and flips of the three-point pattern.
In terms of individual density configurations, each of these correlators acquire positive contributions either when all of (δn i , δn j , δn k ) are positive (+++), or when only one is (+−−), (−+−), (−−+).The CCNN does not directly inform us as to how the correlator became positive-manual follow-up is necessary to uncover this information.For example, in Fig. A14, we show the statistics of these different sign contributions for the dominant three-point correlators learned by the CCNN within the rhombic phase.For conciseness, the heights of the bars corresponding to mixed-sign contributions are averaged over all configurations which produce the same sign.In Fig. A15, we show the value of a handful of these key correlators across the (∆/Ω, R b /a) parameter space in the data.
From this, we can see that the short-range patterns learned by the CCNN are actually producing positive signals due to a large number of (+−−) density configurations, as shown in Fig. A14(a,b).This points to the longer-range packing of Rydberg excitations within this phase-if one the sites in the indicated triples is excited, it is more likely that the other two are in the ground state.Meanwhile, the patterns of Fig. A14(c,d) are positive dominantly due to (+++) configurations, indicating that these motifs signal actual common configurations of joint Rydberg excitations.Figure A14 shows that the star-like configurations also have a large (+++) signal in this phase, as expected from the idealized pattern, but the other sign contributions cause this correlator to turn negative and not be picked up by the CCNN.Together, we can infer from these correlations that we are probing a rhombic-like phase which is failing to develop long-range order due to the incommensurate geometry of the system.
In Fig. A15, we show the extent in parameter space of the key identified connected three-point correlators, and observe that the region colored by the CCNN as the rhombic phase does indeed correspond to the region where the long-range three-point motifs uniquely produce a positive signal.Further theoretical and numerical analysis is needed to confirm that this phase exists in the thermodynamic limit and that the rough region of parameter space identified by our CCNN corresponds to the true region hosting this phase, as well as to determine if these higher-order signals remain good indicators of the phase in the thermodynamic limit.

FIG. 1 .
FIG. 1.(a) Defect-free square lattices of neutral atoms undergo coherent quantum evolution for different values of blockade extent R b /a and linear detuning sweeps' endpoints ∆/Ω, followed by projective readout in which atoms excited to the Rydberg state are detected as loss (red circles).(b-e) Idealized real-space patterns corresponding to phases predicted to be present at various regions of parameter space.Dark pink and white sites indicate |r and |g states, respectively, while the light pink sites in the striated phase are in a quantum superposition of |r and |g .(f ) A diagram outlining the Hybrid-CCNN approach.First, an unsupervised technique is used to generate a rough first-pass phase diagram.Here, we choose to measure average Fourier amplitudes |n(k)| 2 at each (∆, R b ), perform a dimensionality reduction using principal component analysis, and finally cluster using a Gaussian mixture model.The resulting phase diagram informs the starting "seeds" in the parameter space, from which snapshots are sampled in a second supervised stage.We then learn to distinguish these snapshots using interpretable classifiers, from which we can extract refined phase boundaries and key identifying features.

FIG. 2 .
FIG. 2. (a-c) Principal components 1, 2, 5 in Fourier space.(d) Results of the GMM clustering performed in the reduced 10-dimensional PCA space and projected for visualization into the space spanned by (a-c).(e) The CCNN architecture for the supervised learning stage, constructed here up to third order, C (3) α , with three learned filters fα and a learned spatial weighting w(x).First, a density-normalized snapshot δn(x) is convolved with the filters fα to produce a convolutional map C (1) α (x).The δn(x) are zero-padded to allow a convolution of the filters over the entire snapshot.Then, a series of polynomials are applied [Eqs.(2),(3)] to produce maps C (m) α (x), which measure m th -order correlators near each x.These maps are summed with a learned spatial weighting w(x) to produce features c (m) α , which are used by a final logistic layer for classification.(f ) Resulting phase diagram produced by supervised learning, obtained by cropping the classification confidence maps at level-set contours ŷ = 0.75 and overlaying them.
α (a j )δn(x + a j ).The use of local features C (m)

FIG. 3 .
FIG. 3. (a) CCNN-learned region of support for the striated phase, with highlighted boxes indicating the training points.(b) Previously used approximate order parameter detecting the striated phase.Red markers indicate phase boundaries obtained from DMRG simulations on a 9 × 9 array [5].(c) The filters learned by a third-order nonuniform CCNN to identify the striated phase in (a) and the signs on the β (2) α coefficients connecting the corresponding c

( 2 )
α to the output.For ease of display, the filter weights are normalized such that the maximum is 1 within each filter.(d) The spatial weighting w(x) learned by a third-order CCNN identifying the striated phase.A single-pixel outer layer, corresponding to where the filter lands on the zero-padded region, is omitted for clarity.(e) A diagram showing example patches of the idealized striated phase whose correlations are measured by the CCNN of (c,d).

FIG. 4 .
FIG. 4. (a) CCNN-learned region of support for the edge phase, with highlighted boxes indicating the training points.(b)The learned spatial weighting and filters for the model trained to identify the edge-ordered phase, with the spatial extent of the snapshot indicated by the dashed line.The outermost pixels correspond to where the filter is centered on zero-padding but "clip" the edge pixels of the snapshot.For display purposes, the filter weights are normalized such that the maximum is 1 within each filter.(c) Measured performance discrepancy between second-order models with a fixed uniform w(x) = 1 and a freely learned spatial weighting.The central lines and bands show the mean and standard deviation across five randomly initialized models of each type, respectively.(d) Experimental snapshots from the training set, showing • • • motifs primarily only along the single-site border, with the interior highly disordered.(e) Measurement of the third-nearest-neighbor δni,jδni+2,j connected correlation function within the edge and bulk (all sites but the outermost two-site strips) along a cut at R b = 1.46, averaged across translations and other symmetries.
FIG. 4. (a) CCNN-learned region of support for the edge phase, with highlighted boxes indicating the training points.(b)The learned spatial weighting and filters for the model trained to identify the edge-ordered phase, with the spatial extent of the snapshot indicated by the dashed line.The outermost pixels correspond to where the filter is centered on zero-padding but "clip" the edge pixels of the snapshot.For display purposes, the filter weights are normalized such that the maximum is 1 within each filter.(c) Measured performance discrepancy between second-order models with a fixed uniform w(x) = 1 and a freely learned spatial weighting.The central lines and bands show the mean and standard deviation across five randomly initialized models of each type, respectively.(d) Experimental snapshots from the training set, showing • • • motifs primarily only along the single-site border, with the interior highly disordered.(e) Measurement of the third-nearest-neighbor δni,jδni+2,j connected correlation function within the edge and bulk (all sites but the outermost two-site strips) along a cut at R b = 1.46, averaged across translations and other symmetries.

FIG. 5 .
FIG. 5. (a) The region of support for the rhombic phase as learned by the full CCNN model of Fig. 1(n), with 5 × 5 filters and highlighted boxes indicating the training points.(b,c) The two learned 4 × 4 convolutional filters for a simplified model (w(x) = 1, β (n) α ≥ 0) trained to identify the rhombic phase.(d-i) High-weight two-and three-point connected correlators measured by c (2) α , c (3) α resulting from the filters in (a,b).We find β . acknowledges funding support by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC002034.K.W. M.G. and E-A.K. acknowledge support by the National Science Foundation through grant No. OAC-1934714.This research was supported in part by a New Frontier Grant from Cornell University's College of Arts and Sciences.R.S. and S.S. acknowledge support by the U.S. Department of Energy, Grant DE-SC0019030.S.E., T.T.W., H.P., and M.D.L. acknowledge financial support from the Center for Ultracold Atoms, the National Science Foundation, the U.S. Department of Energy (DE-SC0021013 & LBNL QSA Center), the Army Research Office, ARO MURI, an ESQ discovery grant, and the DARPA ONISQ program.
A1) where N (∆, R b ) is the number of available snapshots at the parameter value (∆, R b ).
FIG.A6.The top twelve principal components, shown as weightings in k-space, and the average projection onto each component across the experimental parameter space.
FIG. A13.(a, b, c) Measurements of connected correlations from the striated experimental training dataset (see TableI), performed on each individual nearest-neighbor (NN), next-nearest-neighbor (NNN), or next-to-next-nearest-neighbor (NNNN) bond, respectively.The colors of each site are for visual aid, showing the sites expected to be in the mostly excited (pink), mostly ground (purple), and entirely ground (black) states in the ideal striated limit.(d, e) Averaged correlations across all bonds of different symmetry classes, tracked as a function of ∆/Ω for a cut at R b = 1.56.
FIG. A14.(a-e)The three-point connected correlators identified from the third-order CCNN trained to identify the rhombic phase (a-d) in the main text, along with another notable correlations characterizing the star and rhombic phases (e).For each, we show the contributions to the correlator from three-site density fluctuation patterns of different sign configurations, measured from the snapshots in the rhombic training dataset (see TableI).Contributions are summed across all 8 reflection/rotationtransformed correlators, and averaged across all translations which leave all three sites within the snapshot.Red bars show magnitudes of positive contributions, while blue bars correspond to negative contributions.For simplicity, the bar graph shows the (++−) and (+−−) contributions averaged over all three possible sign combinations.Above each bar graph, we show the total value of the correlator, obtained by summing all contributions with the appropriate sign and multiplying the averaged mixed-sign contributions by 3.

TABLE I .
) where i runs over all snapshots in the training dataset, Parameter-space points used to train models identifying each phase.Training sets are formed by all points in the Cartesian product of the ∆/Ω and R b /a columns.