Finding hidden order in spin models with persistent homology

Persistent homology (PH) is a relatively new field in applied mathematics that studies the components and shapes of discrete data. In this work, we demonstrate that PH can be used as a universal framework to identify phases in spin models, including hidden order such as spin nematic ordering and spin liquids. By converting a small number of spin configurations to barcodes we obtain a descriptive picture of configuration space. Using dimensionality reduction to reduce the barcode space to color space leads to a visualization of the phase diagram.


I. INTRODUCTION
In condensed matter physics we often deal with materials where collective degrees of freedom such as magnetization, density modulations, coherent condensate formation, etc. can undergo order to disorder transitions as a function of external parameters. The distinctly different orderings of these degrees of freedom constitute different phases of the system, characterized by observables known as order parameters.
In the Landau theory of phase transitions, the transition between phases involves the breaking of a symmetry of the system. A classical example is the breaking of rotational symmetry in ferromagnetic materials. In the paramagnetic (high temperature) phase the spins are fully disordered and the system has the symmetry of the crystalline structure. At the onset of ferromagnetic order (low temperature) the symmetry is reduced. Order parameters reveal key symmetries of the phases and serve as indicators for the phase transitions where the condensed matter system undergoes a qualitative change as a result of change in control parameters such as temperature, pressure, and electromagnetic fields. Conventionally, these order parameters are constructed by hand.
In the presence of frustration between atomic spins, magnetic materials feature rich phase diagrams with different forms of collinear and noncollinear antiferromagnetic orderings, as well as more unconventional phases that lack long range magnetic orderings. At the same time we know that both the easily identifiable phases, such as the ferromagnet, and more complicated and not easily visualized orders, like spin nematic, do describe strong correlations between spins. Of relevance for the present work are classical and quantum spin liquids 1-5 , the spin-ice phase [6][7][8] , and nematic ordering that involves breaking of spin rotational symmetry while in the absence of magnetic ordering 7,[9][10][11][12] . The spin nematic ordering is an example of a hidden order parameter which can be difficult to observe directly in an experiment, but that still reveals itself in the heat capacity signature and other observables at the phase transition into nematic state 4,13 .
The commonly used indicators of phase transitions are heat capacity, susceptibility of the relevant order parameter(s), and derived quantities such as the reduced fourth-order Binder cumulant 14 . While a peak in the heat capacity as a function of temperature can reveal the presence of a phase transition, it can not alone be used to characterize the phases. Understanding of the order parameter with concomitant susceptibility is a required step.
A new paradigm to investigate phase transitions based on neural networks recently emerged. It has been shown that neural networks can classify spin configurations sampled from a Monte Carlo simulation and thus can function as order parameter identifiers 15,16 . The neural network is trained in a supervised fashion, i.e. it is given labeled microstates, where the labels are determined by a known conventional order parameter. Intriguingly, recent work has shown that it is also possible to identify phases without prior information of the phase diagram, including phases without conventional order parameters [17][18][19] .
In this work, we use a qualitatively different approach to automatically construct order parameters, based on a persistent homology (PH), a recently developed field in applied mathematics. In this approach the objective is to identify components and shapes given discrete data and a distance metric. These features are captured in a persistence diagram or barcode, two different graphical representations of the same information. Roughly speaking, these diagrams show at what scale different features are present in the data, and over what range of scale they persist.
PH is used in many fields including image processing 20 , complex networks 21 and natural language processing 22 . However, it has not seen much use in condensed matter yet. In cosmology, it has been applied to study the cosmic web and its structure 23 . In materials science, it is used to study nanoporous materials 24 and other larger structures such as silica glass and polymers 25,26 . A recent work proposes persistent homology observables to study dynamics of a quantum many-body system, using the two-dimensional Bose gas as an example 27 . There are two works that discuss the application of persistent homology for the detection phase transitions. Firstly, Donato et al. 28 study the mean-field XY model and a classical Φ 4 model. The phase transitions are detected by computing the persistent homology of the configuration space with molecular dynamics simulations. Secondly, the recent work by Tran et al. 29 demonstrates the detection of phase transitions in the classical XY model and in quantum models.
We demonstrate how PH provides an unsupervised method for constructing phase diagrams, capturing both conventional and unconvential phases with a single framework. Each microstate is converted into a persistence diagram (or barcode), after which two different microstates can be compared with a a similarity metric on their corresponding persistence diagram. Defining a similarity metric (or kernel) for persistence diagrams is still an open problem and we examine the capability of a sliced Wasserstein distance for spin models. This provides a unbiased way to construct a phase diagram.
The paper is organized as follows. We introduce persistent homology as an order parameter in Sec. II, describe our method to construct phase diagrams in Sec. III and the XXZ model on pyrochlore lattice in Sec. IV. The results from the application of persistent homology to the XXZ model are presented in Sec. V followed by a discussion in Sec. VI and conclusion in Sec. VII.

II. PERSISTENT HOMOLOGY AS AN ORDER PARAMETER
Persistent homology (PH) is a method that identifies the qualitative features of a finite metric space (also called a point-cloud dataset) -for an introduction see 30,31 . Given a finite metric space and a distance parameter r, the point-cloud is converted into a graph where edges are added between vertices with distance smaller than r. A simplicial complex K 1 is built from the graph, and the homology of K 1 can be computed efficiently by linear algebra. Performing this construction for increasing r yields a sequence of nested subcomplexes, where the sequence K is called a filtered simplicial complex. The computed homology group H k (K i ) identifies k-dimensional holes in the simplicial complex. For example, k = 0 relates to connected components (clusters) and k = 1 to 1-dimensional holes and so on. The rank of H k counts the number of k-dimensional holes, also referred to as the Betti number β k . The persistence (or lifetime) of k-dimensional holes over the distance parameter r is visualized as a barcode, persistence diagram (PD), or lifetime diagram, all of which represent the same information.
Common types of complexes are the Vietoris-Rips (VR),Čech and Alpha complex 30,31 . The complexes differ in number of simplices, affecting the computational complexity. We examine the capability of Alpha complex (also α-complex) that is often used when the input data is 2D or 3D, in our case spins on a 3D lattice. The α-complex is equivalent to theČech complex when studying persistent homology, but it contains fewer simplices.
The α-complex is constructed from a Delaunay triangulation on the input point-cloud dataset 32 . For a 3dimensional point cloud, the highest dimensional simplex is the 3-simplex (i.e. tetrahedron). Due to the input space being 3-dimensional, there are three homology groups of interest: H 0 , H 1 and H 2 .
With this brief introduction we now move to outline how we apply PH to spin models. The spin configurations (microstates) from classical Monte Carlo simulations form a metric space in multiple ways. Firstly, a single spin configuration can be considered a point in state space (e.g. discrete in the case of an Ising model, continuous in case of a Heisenberg model). Given a distance metric between two states, the corresponding persistence diagram (PD) is describing the topology of the phase diagram, where the homology groups describe the shapes of the phases in parameter space. Secondly, a single site in the spin system can be considered a point, leading to a single PD per microstate. A phase diagram can be constructed by comparing PDs.
In this work, we focus on the latter concept, i.e. the homology groups describe a microstate (spin configuration) directly. Averaging over microstates is possible by simply merging barcodes and discarding the the information from which microstate each bar (or point in PD) originated.
The concept introduced above requires a distance metric between two sites, where we can use the information of the lattice and the spin. The symmetry of the barcode is dictated by the choice of distance metric and determines its capability of capturing distinct phases.
As a simple example, consider the ferromagnetic square lattice Ising model (J = 1) hosting a lowtemperature symmetry-breaking ferromagnetic phase and a phase transition at T c ≈ 2.269 J to a paramagnetic phase 33 . Choose a distance parameter r where aligned spins are closer than opposite spins, with the distance being an Euclidean distance between spin tops. This definition of a distance metric leads to domains forming connected components in homology group H 0 . Hence the Betti number β 0 (r) counts the number of domains present in the spin configuration, similar to the conventional order parameter of magnetization density M . In general, any metric capturing the shapes present in the spin texture is able to identify the two phases.
The phase diagram of a 2D Ising model discussed above is a convenient example, but the phase diagram is too simple to benefit from an approach with PH. In this work we study a more complex system, namely the XXZ model on a pyrochlore lattice, known to host a large number of competing phases 4 .

III. PHASE DIAGRAM CONSTRUCTION
Traditionally, phase diagrams are constructed by observing a change in the order parameter that typically ranges from zero (phase A) to finite (phase B). In our proposed scheme, shown schematically in Fig. 1, the persistence diagram (PD) plays the role of the order parameter, and therefore we need to quantify the change in the PD.
Common distance metrics for PDs are the bottleneck distance and the Wasserstein distance 34 . Both essentially measure the similarity by attempting to match the points (i.e. persistent features) in the two diagrams. These metrics are too computationally expensive for our purposes. Instead, we use the sliced Wasserstein (SW) distance which is an approximation of the Wasserstein distance 35 . The distance is measured between two PDs of same kth homology group, i.e. SW(PD k i , PD k j ). For the total distance, we consider holes (H 1 ) and voids (H 2 ): The distance matrix D is established by calculating the SW distance between all pairs of system parameters (e.g. all combinations of interaction J and temperature T ). Finally, the phase diagram is visualized by dimensionality reduction on D to a 3-dimensional color space (red, green and blue). There are obvious drawbacks to this simple approach such as the non-linearity of color perception, but we find it sufficient. When distances are small and more color space is required, it is possible to limit the visualization to a subregion of the phase diagram, as shown in Sec. V. Dimensionality reduction algorithms aim to construct a low-dimensional image where distances of the original high-dimensional space are preserved as best as possible. Principal component analysis (PCA) is a common method based on matrix factorization, when given a set of input data points. In our case we only have a distance matrix D. Therefore PCA cannot be used directly, but multidimensional scaling (metric MDS) is the equivalent basic technique in this case 36 . In summary, metric MDS algorithm is used to reduce the distance matrix D to three RGB color channels that color-code the different phases in the phase diagram. Similar phases with small SW distances will appear close in color space and are therefore color-coded with similar colors.

IV. XXZ MODEL ON THE PYROCHLORE LATTICE
The pyrochlore lattice is a cornerstone in research on frustrated magnetism, and is having a lead role in experimental and theoretical explorations of spin-ice physics 1,4,[6][7][8]  pyrochlore materials, a variety of spin Hamiltonians with long or short ranged interactions on the pyrochlore lattice have been investigated. The XXZ model on a pyrochlore lattice is a system with short ranged interactions having a rich phase diagram with competing antiferromagnetic, spin-ice, spin-liquid and spin-nematic phases 4,6 . This model has also been used previously to test a machine learning model (support vector machine, SVM) that identifies phases 17 . The Hamiltonian of the XXZ model is given by with S ± i = S i,x ± iS i,y and S i = (S i,x , S i,y , S i,z ), S i = 1, and has the cubic symmetry of the pyrochlore lattice. In the following, energies and temperatures are expressed in terms of the antiferromagnetic J zz = 1 exchange interaction.
Extensive classical Monte Carlo simulations, field theoretical analysis, and spin dynamics simulations by Taillefumier et al. 4 have established a J ± − T phase diagram with six phases which here will be briefly recapitulated (numbering corresponding to 17  IV Easy-plane spin-nematic (SN⊥) with algebraic spinspin correlations, no long range spin dipole ordering. The rotational U (1) symmetry of the local z axes is broken by the onset of a higher order multipolar ordering.
V Spin ice (SI). Each tetrahedron on the lattice has the "two-in, two-out" spin configuration with the spins aligned along their local z axis (see Fig. 2 (a)). The phase has algebraic spin correlations, no long range spin ordering.
VI Pseudo-Heisenberg antiferromagnet (pHAF) with algebraic spin correlations distinct from phases III and V, no long range spin ordering.
In summary, one phase has long range antiferromagnetic ordering (I AF⊥), one phase is disordered (II PM), one phase has spin-nematic ordering (IV SN⊥), and three of the phases (III SL⊥, V SI, and VI pHAF) are classical spin liquids. Relative to the paramagnetic phase, in which the system has the symmetries of the pyrochlore lattice, breaking of symmetry occurs only on entering the easy-plane antiferromagnetic phase (I AF⊥) or on entering the easy-plane spin-nematic phase (IV SN⊥), hence only when crossing in and or out from one of these two phases does the system go through a phase transition, with characteristic peak structure in the heat capacity and order parameter susceptibilities 4 . Transitions into or from one of the three spin liquid phases constitute a gradual evolution (crossover) of spin configurations, involving no breaking of symmetry, and lacks sharp peak structures in the heat capacity. The latter quantity nevertheless contains information that can be used to define a distinct criteria for where the crossover occurs 4 . The three classical spin liquid phases can be distinguished by the different form of algebraic spin correlation reported for each of the phases 4 , correlations that can be measured in neutron scattering experiments. Rather remarkably, the easy-plane spin-nematic phase have the same algebraic spin correlation as the easy-plane spin liquid phase, yet the phases differ in the regard that in the latter the rotational symmetry in the local z axes is broken by the nematic order parameter.
To test the proposal of PH for detection of the phases, we use representative spin configurations as an input.

V. PERSISTENT HOMOLOGY ON THE XXZ MODEL
A single spin configuration is converted into a point cloud by placing a point at each spin tip, where the spin length is set to 1/4 of the lattice tetrahedron edge length d = a 2 √ 2 . Euclidean distance is used once the point cloud is established using the distance metric where r(i, j) is the Euclidean distance between two sites. This corresponds to the distance of tips of the spin arrows (see Fig. 2). Note that S i is in the global crystal frame, not in the local frame of Equation 3 (see Appendix A for coordinate frame conversion). For nearest neighbors, this means D(i, j) is in the range (d/2, 3d/2). Figure 2 shows the exact barcode for a single unit cell in the spin ice (V SI) phase, where spins are set ac-cording to the "two-in, two-out" rule. The connected components (H 0 ) reveal that there are two length scales present. The smaller length scale (α 2 ≤ 0.010) corresponds to neighboring spins forming 1-simplices (e.g. Fig. 2 (a)). The larger length scale (α 2 ≥ 0.034) leads to a single connected component, as shown in Fig. 2 (b). The 1-dimensional holes (H 1 ) appear at α 2 ≥ 0.034 as 9 bars. This includes 8 empty triangles (i.e. complex of 1-simplices rather than a single 2-simplices) and one large hole (also visible at α 2 = 0.05, Fig. 2 (b)). The empty triangles are filled at α 2 = 0.037 and become 2simplices, shown as shaded triangles in Fig. 2 (b). At α 2 = 0.062, the 1-simplices form a large tetrahedron with 2-simplices at its corners. A tetrahedron of 1-simplices has β 1 = rank(H 1 ) = 3, which corresponds to 3 bars for H 1 in Fig. 2 (c). All 1-dimensional holes are filled at α 2 ≥ 0.131. At this point, a single 2-dimensional void (H 2 ) appears for α 2 = [0.132, 0.168], shown in Fig. 2 (d).
Generally, each phase has persistent features at different characteristic length scales. We sample 32 sta- tistically independent spin configurations for a single (T /J zz , J ± /J zz ) parameter combination. The αcomplexes are calculated by GUDHI 32 . The resulting barcodes are merged into a single barcode, discarding the information of the microstate origin from each bar. The similarity between barcodes is measured by the sliced Wasserstein (SW) distance, implemented in persim 38 , see Equation 2 .
The distance matrix D is dimensionality reduced with MDS to a 3-dimensional space which is interpreted as red, green and blue channels. Figure 3 shows the phase diagram. The paramagnetic phase II PM and the phases V SI and I AF⊥ are clearly distinct due to large pairwise sliced Wasserstein distances. The small distances between III SL⊥ and IV SN⊥ require a rescaling of the colormap to be visible, in the separate left pane. Also, the VI pHAF phase has small SW distances to the paramagnetic phase and is therefore difficult to detect, but a small color gradient is present. Figure 4 shows the lifetime diagrams and confirms that the phases II PM, III SL ⊥ , IV SN ⊥ and VI SN ⊥ are similar, leading to smaller SW distances. A more detailed view of the lifetime diagrams as temperature decreases is shown in the Appendix Fig. 6. Note that the spin ice phase V SI exhibits two H 1 features (with lifetimes around 0.07 and 0.1), corresponding to the two longer H 1 bar lengths in Fig. 2. In general, each phase has its own characteristic features and differences in its fingerprint are captured by the sliced Wasserstein distance. The qualitative similarities of phases and significant differences between lifetime diagrams, I AF vs V SI vs IV SN ⊥ vs III SL ⊥ vs VI pHAF phase, leads to the proposal that these lifetime diagrams can serve as respective order parameters (i.e. lifetime PH portraits of the phases).

VI. DISCUSSION
The PH method outlined here is using only 32 spin configurations for each point in parameter space (T /J zz , J ± /J zz ) and does not require data on energies or heat capacities, making it an efficient technique. The α-complexes are faster to calculate than theČech and Vietoris-Rips complex and capture the holes and voids created by spins locally ordering. The computational bottleneck is the calculation of sliced Wasserstein distances between persistence diagrams. This could be improved by a faster distance metric for persistence diagrams, which is still an open problem in persistent homology. Finally, we note that persistent homology has very few free parameters in comparison to neural networks, which facilitates the interpretation of results.
In comparison to the work by Greitemann et al. 17 our method has a couple of key differences. Their work relies on the support vector machine (SVM) and the construction of monomials of spin components. We use persistent homology instead, involving no regression coefficients, and operate directly on the spin configurations. For the phase diagram construction, Greitemann et al. use the Fiedler vector to partition the phase diagram given a distance matrix D. We instead use the simpler technique of multidimensional scaling (metric MDS) directly on the distance matrix D.
Recent work has shown that PH can detect phase transitions in a variety of classical and quantum models 28,29 . We show the universality of the approach by identifying 6 different phases at once, including three classical spin liquids and construct a full two-dimensional phase diagram.
The detection of phases is possible by studying persistence diagrams (Fig. 4) and changes are quantified by the sliced Wasserstein distance. We reduce the persistence diagrams space to color space with matrix factorization, which proves to work best when SW distances are large and of similar magnitude. Alternatively, dimen- sionality reduction based on the neighbor graph approach (e.g. Uniform Manifold Approximation and Projection, UMAP 39 ) could highlight local structure better at the expense of representing global structure of the original high dimensional barcode space.

VII. CONCLUSION
We demonstrate that persistent homology (PH) can be used to capture different types of order in spin models. The application of our method to the XXZ model on the pyrochlore lattice with its six phases demonstrates the versatility of this approach. Both the phase with long range order and the phases with only local ordering are identified.
The barcodes reveal the characteristic length scales present in the spin model. Using sliced Wasserstein distance, we obtain a distance matrix for all pairwise system parameters. The distance matrix can be visualized using dimensionality reduction, revealing similar regions in parameter space. Alternatively, the persistence diagrams can be inspected directly to observe changes as the system temperature decreases.
In summary, persistent homology provides a new general computational framework to study both long range and local spin ordering. Extending the persistent homology framework to quantum spin models will be the topic of future work.