Network science Ising states of matter

Network science provides very powerful tools for extracting information from interacting data. Although recently the unsupervised detection of phases of matter using machine learning has raised significant interest, the full prediction power of network science has not yet been systematically explored in this context. Here we fill this gap by providing an in-depth statistical, combinatorial, geometrical and topological characterization of 2D Ising snapshot networks (IsingNets) extracted from Monte Carlo simulations of the $2$D Ising model at different temperatures, going across the phase transition. Our analysis reveals the complex organization properties of IsingNets in both the ferromagnetic and paramagnetic phases and demonstrates the significant deviations of the IsingNets with respect to randomized null models. In particular percolation properties of the IsingNets reflect the existence of the symmetry between configurations with opposite magnetization below the critical temperature and the very compact nature of the two emerging giant clusters revealed by our persistent homology analysis of the IsingNets. Moreover, the IsingNets display a very broad degree distribution and significant degree-degree correlations and weight-degree correlations demonstrating that they encode relevant information present in the configuration space of the $2$D Ising model. The geometrical organization of the critical IsingNets is reflected in their spectral properties deviating from the one of the null model. This work reveals the important insights that network science can bring to the characterization of phases of matter. The set of tools described hereby can be applied as well to numerical and experimental data.


INTRODUCTION
Networks [1][2][3][4][5] encode the information present in a large variety of natural and artificial interacting systems by representing them as graphs, i.e. a set of nodes (representing the element of the system) connected by links or edges (representing typically the interactions).In particular, the underlying architecture of complex systems is encoded in networks that strongly deviate from random graphs whose information content can be mined by exploiting their statistical, combinatorial as well as the geometrical and topological nature [6].Networks are hence a simple yet very powerful framework that has been able in the last twenty years to transform our understanding of complex systems.These complex networks obey relevant organization principles while retaining a stochastic nature.
Recently great attention has been addressed to formulating unsupervised machine learning algorithms to detect different phases of matter [7][8][9][10][11][12].In this research line, important progress in characterizing phase transitions has recently been made by using interpretable machine learning tools such as graphical models and restricted Boltzmann machines [13][14][15][16][17].The core of such approaches is that learning methods -in particular, unsupervised -shall be capable of revealing hidden struc-tures in data sets, that correspond to physical information (such as, e.g., the existence of order parameters).However the very advanced, complementary tools of network science to extract information from complex data of interacting systems have not yet been systematically employed for this task.
Here we want to show how network science can enrich and enhance our unsupervised characterization of phases of matter.Indeed we will show evidence that network science provides a very transparent set of methods to investigate the characteristics of different phases across critical phase transitions and allows the determination of unsupervised indicators of their critical points.
Historically networks have been used in condensed matter for describing physical interactions existing among the elements of a system, as well as structured in configuration and Hilbert space.In principle, they can be used as well to represent abstract data structures coming from numerical simulations or directly from experiments, giving access to a whole new toolbox to define and interpret many-body correlations (of arbitrary order, if expressed in terms of local observables).
Reflecting this two-fold possible application of network science, recently several works have explored the use of networks to model or represent interactions in condensed matter systems.For instance, networks can be consid-ered to define Hamiltonians whose interaction terms are determined by a complex network rather than by a lattice.In particular networks can be used to define different critical phenomena including the Ising model [18][19][20] and the inverse Ising model [21], as well as quantum critical phenomena including the transverse Ising model [22][23][24], the Bose-Hubbard model [25], the Jaynes-Cummings-Hubbard model [26], in addition to others classical collective phenomena and inverse problems [27][28][29][30].Alternatively, networks can define quantum environments [31,32] or even multilayer couplings between interdependent superconductor networks [33].Networks also can be used to represent correlations in complex and financial networks [34,35] as well as in quantum systems.In particular recently network structures have been shown to encode the quantum long-range mutual information (and other measures of quantum correlations) existing among the nodes of quantum lattice models in one [36][37][38] and two dimensions [39], providing in some cases indicators for quantum critical points.
As we will see in this work, weighted networks [40] are amenable to be analysed and treated with Topological Data Analysis (TDA) [41][42][43][44] which provides a very efficient way to probe topological, large scale and global network properties.TDA, although until now only applied to point clouds, is raising significant interest to do unsupervised inference of phase transitions [45][46][47][48][49][50][51], and has wide applications, including the characterization of universal dynamics in quantum gases [52] and of confinement in lattice field theory [53][54][55].
Finally and most relevantly for our work, networks have been proposed to capture the underlying structure of quantum spin systems as revealed by wave function snapshots that can be probed experimentally as well as sampled from Monte Carlo simulations [56].However, despite these very pioneering works [36,56], little attention has been so far addressed to study phases of matter using network science.
Here we launch a large-scale systematic study of the phase of matter based on network science.We leverage a multiplicity of tools developed in network science and we reveal the combinatorial, statistical, geometrical and topological network representation of different phases of matter.We provide an in-depth characterization of the networks generated from single snapshots of spin system configurations.
In Ref. [56] it was shown that wave function networks constructed starting from quantum wave function snapshots are strongly deviating from random graphs and for a wide range of values of the threshold distances they give rise to networks with very broad degree distribution.An open question is whether the complex properties of these networks are inherently quantum effects or they can be observed in classical systems as well.
To characterize phase transitions and critical behaviors, the Ising model is arguably the most studied model in statistical physics.Apart from the Ising model defined on lattices, the model has also been generalized on networked structures such as small-world networks [57,58], random scale-free networks [18-20, 59, 60] and spatially embedded scale-free networks [61].The phase diagram of the model is shown to be highly sensible on the value of the branching ratio of the network.These results have been also extended to the Transverse Ising models [22,23].Moreover some of the machine learning tools to study the Ising model such as restricted Boltzmann machine and graphical models [13][14][15][16][17] strongly leverage on their underlying network structure.However, to the best of our knowledge, the tools from network science have not been used to the analysis the simulation snapshots of the Ising model.
In this work we consider 2D Ising snapshot networks (IsingNets) following a construction proposed in Ref. [56] applied to classical 2D Ising model snapshots and we characterize their structure using advanced statistical, combinatorial, geometrical and topological tools of network theory.In order to provide an in-depth analysis of the IsingNet, across different phases we focus our attention on IsingNets obtained starting from Monte Carlo simulations of the 2D Ising model performed across the phase transition.IsingNets are obtained from a sample of state configurations of the 2D Ising model which constitute the set of nodes of the IsingNets.Each pair of nodes of an IsingNet is associated with a distance, here taken to be the Euclidean distance between the configuration snapshots.IsingNets are constructed starting from the fully connected distance matrix between the nodes by connecting only the nodes whose distance is smaller than a threshold value of the distance.
Our in-depth network analysis of the IsingNets will allow us to go well beyond the characterization of these networks based solely on the degree distribution.Possibly in the future, this in-depth analysis can be conducted also on networks built from quantum wave function snapshots in order to assess which are the properties inherently quantum in the latter networks.
Our analysis is conducted following two main directions whose goal is different but complementary.First, we will perform an analysis of the IsingNets that is agnostic about the choice of the distance threshold.In particular, we will study network properties as a function of the distance threshold.These include percolation properties, persistence homology, network embedding and statistical characterization of the distance matrices.Secondly, we will consider specific choices of the distance threshold and we will characterize the statistical, combinatorial and geometrical/spectral properties of the IsingNets, showing the important roles of degree-degree and weight-degree correlations in these systems.
Anticipating our main results we have found that IsingNets reflect the symmetry of the configuration space of the 2D Ising model in a prominent way.In particular, the percolation properties of the IsingNets strongly deviate from the percolation properties of networks in which the same distances among the nodes are distributed randomly.In fact, below the critical temperature of the 2D Ising model, the IsingNets are characterized by two giant components whereas their randomized counterpart displays only a single giant component.When the threshold distance defining the IsingNets is raised significantly these two giant clusters merge, but interestingly they keep a rather compact structure as revealed by our persistent homology results highlighting that the Betti number of their clique complex are strongly suppressed with respect to their random counterpart.The Weisserstein distance between the persistence diagrams of the IsingNets and their randomized couterpart is here shown to be an unsupervised indicator of criticality as it displays a maximum in correspondence of the critical temperature.Our statistical and combinatorial analysis of the IsingNets strongly demonstrates the complex organization of these networks which display strong heterogeneity in both their topological (degree, clustering coefficient, degree correlations) and their weighted network properties.In particular nodes of higher-degree are characterized by having neighbours connected by stronger affinity weights (smaller distances).Finally, the IsingNets possess a significant geometrical organization at criticality as it is revealed by their interesting spectral properties.
Among the most important benefits of our approach is that the proposed analysis is fully interpretable.Moreover, we emphasize the applicability of the approach.Indeed while our analysis is here performed on data coming from Monte Carlo simulations the approach can be readily applied also to experimental data.However, a possible limitation of the approach is due to the significant computational cost of our proposed unsupervised TDA analysis.

THE 2D ISING MODEL MONTE CARLO SIMULATIONS
We consider a square 2D lattice of dimension L × L where on each site n is located the spin S n ∈ {−1, 1}.The nearest neighbour spins are interacting through the Hamiltonian The 2D Ising model is characterized by Z 2 spontaneous symmetry breaking and undergoes a second-order phase transition at T c = 2/ ln(1 + √ 2) ≈ 2.269 [62].Starting from Markov Chain Monte Carlo simulations of this model, for each temperature single snapshots ⃗ x i = {S 1 , S 2 , . . ., S L 2 } of the spin system are sampled at equilibrium [63].More specifically, we use the Wolff cluster algorithm [64,65], starting from the configuration with either all up spins or all down spins, chosen at random.Next, 30000 to 50000 'cluster flips' are performed for the system to equilibrate.After this, we collect snapshots every 1000 to 1500 cluster flips to ensure that the collected state configurations are as uncorrelated as possible [63].In total, we gather 10000 snapshots during a Monte Carlo run.
For each temperature, five independent Monte Carlo simulations are performed as prescribed.By combining the sampled configuration snapshots of these runs, we thus obtain a data set with N r = 50000 independent thermal configurations {⃗ x i } Nr i=1 .The starting point to construct the IsingNets is a set of N configuration snapshots i ∈ {1, 2, . . ., N } randomly selected from the data set described above, and the fully connected distance matrix d of elements d ij between these states.Here the distance d ij between two generic snapshots ⃗ x i and ⃗ x j is taken to be their Euclidean distance.As discussed in previous works, such manifolds are typically living in very high dimensional subspaces [12,63], so that simple dimensional reductions are not applicable, and a full-fledged network analysis is needed.

NETWORK CHARACTERIZATION ACROSS THE DISTANCE FILTRATION Weight filtration
As anticipated in the introduction, in this first Section, our analysis focuses on the properties of IsingNets observed as a function of the distance filtration.We consider IsingNets which are graphs G = (V, E) formed by a set of N nodes V and a set of links E with (i, j) ∈ E only if the nodes (state configurations) i and j have distance d ij < r.Here r determines the distance filtration and indicates a tunable parameter that ranges from the minimum of the distances r min = min i,j d ij between the N nodes to their maximum distance r max = max i,j d ij .We are thus here completely agnostic about the best choice of r as we study the properties of IsingNets across all possible choices of the distance threshold r.
In particular as a function of r we will explore the percolation properties of the IsingNets which define an agglomeration from N disconnected nodes to a single connected component as r is raised from r min to r max and we will compare this process with the corresponding null model obtained from the same process applied to a randomized distance d rand matrix.The randomized distance matrix d rand is constructed by randomly reshuffling the upper triangular elements of d and subsequently symmetrizing the matrix.Therefore the null model networks display the same distribution of "distances" as the true IsingNet while being completely randomized.Note however that one of the main differences between the distance matrix d and the randomized distance matrix d rand is that the entries of d rand are not proper distances as they do not obey the triangular inequality.
In this section, we will use a combination of tools coming from network science to analyse the IsingNets described above.In particular, anticipating the detailed description of the methods used to perform this analysis in the following paragraphs, we will use persistent homology [42][43][44] to further characterize topologically the mentioned aggregation process.In this way, we will show that persistent homology is able to detect the position of the critical temperature of the 2D Ising model under study.This analysis will be accompanied by the visualization of the network using Minimum Spanning Trees [35,66], the results of the network embedding conducted using the UMAP (Uniform manifold approximation and projection for dimension reduction) algorithm [67], and the statistical characterization of the distance matrix as a function of the temperature conducted using the closeness centrality [68] distribution.

Percolation process
We start our investigation exploring the percolation process [2,3,[69][70][71] monitoring the connected component of the network as a function of the filtration parameter r.To contrast the behavior of the percolation process of IsingNets with a null-hypothesis percolation process, we consider the process defined on the actual IsingNet distance matrix d with the same process defined starting from a matrix d rand obtained by randomly permuting distances among pair of nodes.
The percolation process of IsingNets reveals a major difference with the percolation process on the randomized distance matrix: mainly the IsingNets obtained for the 2D Ising model below the phase transition, i.e.T < T c display for a very significant range of values of the filtration parameter r, two giant components while the randomized process only displays one giant component.This phenomenon is evident from the plot in Figures 1 and  2 showing the relative size of the largest component R and the second largest component R 2 , which are both giant, i.e. extensive for a wide range of r values.Indeed below the critical temperature T c , the IsingNets display two transitions as the value of the filtration parameter is raised (see Figures 1 and 2).The first transition is characterized by the emergence of the two equal size giant components corresponding to the symmetry of the configuration snapshots of the 2D Ising model for T < T c and the second one is characterized by the merging of these two giant components for very large values of r, characterized by the disappearance of a significant second largest connected component (orange line in Figure 1 (a) and Figure 2 (a).This phenomenology is dramatically different from the percolation obtained in the randomized null model where the giant component is unique for every value of r (see Figures 1 and 2).For temperatures above the critical one (see Figure 3), instead, only one giant component is observed corresponding to the paramagnetic state of the 2D Ising model.In order to further characterize the percolation process we also monitor as a function of r the average size of finite components ⟨s⟩, the number of components n s and the inverse participation ratio Y whose inverse determines the number of typical clusters.Specifically the inverse participation ratio Y is defined as   On an Erdös-Renyi random graph the average size of finite component ⟨s⟩ plays the role of the percolation susceptibility [69], diverging in correspondence with the emergence of the (single) giant component.In the IsingNets (see Figures 1-3) the average size of finite com-ponent ⟨s⟩ has a very suppressed maximum with respect to the same quantity measured on the considered null model, indicating that the agglomeration of the two giant clusters proceeds by subsequent agglomeration of very small components and isolated nodes rather than by the agglomeration of finite clusters of diverging average size as in a random graph.This is also confirmed by the behavior of the number of clusters n s as a function of the filtration parameter r which for the IsingNets decays less steeply than in the randomized null model.Finally, the inverse participation ratio Y for low temperatures reveals a significant plateau at Y = 1/2 indicating the existence of two giant components of approximately equal size (see Figure 1).Above the critical temperature, for T > T c as the filtration parameter is raised only one giant component emerges and the difference with respect to the randomized null model is reduced (see Fig. 3).

Persistent homology
Topology is the study of shapes and their invariant properties under continuous deformations (see for an introduction [6,41]).Major examples of topological invariants are the Betti numbers.The Betti number β 0 indicates the number of connected components, the Betti number β 1 indicates the number of one-dimensional holes, the number of β 2 indicates the number of twodimensional holes, etc.For instance, a point has Betti numbers β 0 = 1 and β n = 0 for any other value of n, a circle has non-zero Betti numbers β 0 = β 1 = 1 and a sphere has non-zero Betti numbers β 0 = β 2 = 1.An important result of algebraic topology [6,[41][42][43][44] is that the n-dimensional Betti number β n is the rank of the n-dimensional homology group of the considered topological space.
In the discrete setting, the Betti numbers are defined in general for simplicial complexes.Simplicial complexes are a type of higher-order network formed by a set of simplices such as nodes, links, triangles, tetrahedra, etc.They have the additional property of being closed under the inclusion of faces of each simplex.This last property implies that if a triangle belongs to the simplicial complexes also all its links and nodes belong to the simplicial complex.
Topological data analysis [40][41][42][43][44]72] and in particular persistent homology allows to characterize the topological properties of data as a function of a filtration parameter and is becoming increasingly important in network and data science with applications ranging from the study of gene-expression to the investigation of brain networks.As mentioned in the introduction TDA is recently becoming a very popular computational tool to study phase transitions as well [45][46][47][48][49][50][51][52][53][54][55].However, in this context, most of the TDA so far are performed on point clouds rather than on networks.In our setting, when each pair of nodes (i, j) is assigned a distance d ij , persistent homology characterizes the topology of data by forming a simplicial complex representation of the data and characterizing its homology, i.e. the connected components (H 0 homology classes); the independent cycles -one-dimensional holes-(H 1 homology classes); the twodimensional holes (H 2 homology classes) etc.The simplicial complexes [6,40,42,43] that we use to perform the topological data analysis are the so-called Vietoris-Rips complexes of the network generated by filling all the simplices having all links at distance d ij < r.Thus as a function of r the set of simplicial complexes forms a filtration.
As the filtration parameter r is raised, first each node belongs to a different connected component, then connected components merge progressively as described also by the previous paragraph (see Figure 4 for a schematic description of the filtration).However, as r increases there is the potential also for one-dimensional holes (or network cycles) to emerge with each independent cycle represented by a different H 1 homological class.Eventually, as r increases these cycles will become filled and thus disappear.Moreover, also higher-dimensional holes represented by higher dimensional homological classes H n might arise and eventually coalesce as r is further increased giving rise to barcodes representing the topology of the data.
By monitoring the homology classes as a function of r, the results of persistent homology are typically summarized by barcodes where each homology class is represented as a bar that extends through the corresponding range of values of r for which the homology class is observed (see Figure 4).The barcodes are then represented by persistent diagrams where for each homology class the value of r where the homology class is disappearing (death) is plotted versus the value of r corresponding to the first appearance of the homology class (birth), see for instance Figure 5. Significant homology classes are the ones represented by points far from the diagonal which last for a large interval of values of the filtration parameter r.
The investigation of the persistent diagrams for the IsingNets can further characterize these discrete structures by revealing important properties about how the clusters merge as a function of r.In particular, the persistent diagrams of the IsingNets allow us to show that clusters remain compact with a suppression of the Betti numbers with respect to the corresponding persistent diagram of the randomized null models.This result is evident from Figure 5 and Figure 6 where we plot the persistent diagram (for homology H n with n ∈ {0, 1, 2, 3} and the Betti numbers β n with n ∈ {0, 1, 2, 3} as a function of the filtration parameter r for IsingNets and their randomized null models as a function of the temperature T .Indeed the persistent diagrams of IsingNets corresponding to homology H n with n ∈ {0, 1, 2, 3} (see Figure 5) where β [I] (r) and β [N ] (r) indicate the Betti number of the IsingNet and the corresponding null model with filtration parameter r.These distance measures provide good indicators of the critical points as they display a maximum approaching the critical temperature of the 2D Ising model as the size L of the 2D lattice increases (see Figure 7).Note that here, due to the computational cost of calculating the persistent diagrams corresponding to higher-order homological classes, we focus here only on homology classes H 0 and H 1 .This is rather clear evidence that the IsingNets have a topology that encodes fundamental properties of the underlying spin system.

Minimum Spanning Tree Visualization
Our unsupervised analysis of the IsingNets includes their visualization which reveals their highly heterogeneous structure below the critical temperature T c and at criticality.A very efficient way to visualize the IsingNets is by plotting their Minimum Spanning Trees (MST) [35,66].The MST is the subtree of the network whose sum of the distances between the connected nodes is minimal, and its topology reveals important properties of the fully connected weighted distance matrix of the IsingNets.In particular, as shown in Figure 8 for T deep in the ferromagnetic phase the topology of the MST of the IsingNets is dominated by very relevant hub nodes, and the network displays a clear partition between the two clusters detected by the K-means algorithm with K = 2 indicated in the figure by two different colors of the nodes.As the temperature is raised to the critical region, T ≃ T c the hubs of MST become less dominant.Above the phase transition, the MST becomes clearly more random with a suppression of the hubs in the MST.

Network embedding
Our analysis of the IsingNets is here enriched by considering the UMAP 2-dimensional embedding of the fully connected network endowed with the distance matrix d.UMAP is a widely used embedding algorithm exploiting dimension reduction (for more detail see for instance Ref. [73]).The embedding is here performed as a function of the temperature T and the nodes are colored according to their clustering in two groups performed using K-means (see Figure 9).
The UMAP embedding provides a clear visualization of the two main clusters of the nodes of the IsingNets present for T < T c and corresponding to the symmetry among configuration snapshots and of their merging as the temperature T is raised above the critical temperature.

Metric-based centrality measures
To conclude our characterization of the IsingNets without imposing a fixed value of the filtration parameter, we present here the statistical properties of some important geometrical aspects of the IsingNet captured by the closeness centrality of the nodes.
The closeness centrality Cl i of a node i measures how close is the node to the other nodes of the network, and is defined as the inverse of the average distance of node i to the other nodes of the network, i.e.
We investigate the statistical properties (the mean ⟨Cl⟩, the standard deviation σ(Cl), the skewness Sk(Cl) and the kurtosis Ku(Cl)) of the distribution of the closeness centrality in IsingNets as a function of the temperature T .
In Figure 10 we show that the average closeness centrality decreases as a function of the temperature, demonstrating that on average the distances between the nodes are higher at higher temperatures.The higher-order statistics of the closeness centrality distribution are even more revealing of the IsingNets organization and the skewness and kurtosis of the closeness distribution provide a good unsupervised indicator of the critical point (see Figure 10).Indeed the standard deviation of the closeness centrality displays a maximum for temperatures close to the critical temperature; the skewness of the closeness centrality is negative below the critical temperature and positive above the critical temperature and the kurtosis becomes negative above the critical temperature, while at very high temperature is strongly affected by noise.

NETWORK CHARACTERIZATION OF ISINGNETS AT A GIVEN VALUE OF THE THRESHOLD DISTANCE r
IsingNets at given threshold distance r In this section, we study the properties of IsingNets where we fix a given choice of filtration parameter r.In particular we will consider the IsingNets, whose N × N adjacency matrix A has elements with θ(x) = 1 for x > 0 and θ(x) = 0 otherwise.In the following, we will indicate with i ∼ j two neighbour nodes for which A ij = 1.
We adopt the same choice of the parameter r proposed in Ref. [56] where it was shown that for a wide range of choices of r the IsingNets are scale-free presenting often power-law exponents less than two which are known to occur in a variety of context [74][75][76][77].In particular, here we take r equal to the average distance of the 5 th nearest neighbour.
Similarly, the randomized network forming our null model in this section is performed by thresholding the randomized distance matrix d rand with the same threshold used for the corresponding IsingNet.
We consider the statistical and combinatorial properties of these networks where we assign to each link (i, j) a weight w ij equal to the inverse of the distance d ij , provided this distance is smaller than r, i.e.
We provide an in-depth network analysis of these networks investigating their degree and strength distribution, and going beyond these statistical properties providing evidence of the presence of degree correlations, of a nontrivial k-core structure and of weight degree correlations.Moreover, the investigation of the spectral properties of the IsingNets will demonstrate non-trivial signatures of criticality.This analysis provides clear evidence that not only the degree distribution of IsingNets is strongly different from the one of an Erdös-Renyi random graph, but the network also obeys important higher-order correlations reflecting the correlations existing in the spin system configuration snapshots.

Degree and strength distribution
One of the most simple yet fundamental property of a network is its degree distribution P (k) which characterize globally the network starting from the knowledge of the node's degrees where the degree k i of the node i indicates the sum of the links incident to it, i.e.
Thus while the degree is a local property of the nodes the degree distribution P (k) is able to characterize the global properties of the networks.In particular scale-free [78][79][80] and in general degree distribution with second k 2 (and eventually first ⟨k⟩) moment diverging with the network size have been shown to have a significant role in determining the response of the network to random damage and the critical behavior of the dynamics defined on top of these networks, such as epidemic spreading and the Ising model [1,79].
For weighted networks it is also possible to define weighted degree also called strength s i of the generic node i [81] given by the sum of the weights of its incident links, i.e.
From the sequence of the node's strengths, it is possible to extract also the strength distribution P (s/s 0 ) where s 0 is the minimal strength of the links.This distribution has been shown to display broad distributions in a number of real weighted networks, such as collaboration networks and airport networks [82].
In Figure 11 we plot the degree P (k) and strength P (s/s 0 ) distribution for temperatures below and above the phase transition demonstrating that these distributions are broad.We note that below the critical temperature T c , the networks are not only broad but also dense, i.e. having an average degree growing with the network size (for models of these networks see [74][75][76][77]).
To quantify the scale-free nature of these distributions we plot the ratio between the second moment, the average degree ⟨k⟩ and the average strength ⟨s⟩ as a function of the temperature, showing that that ⟨k⟩, and ⟨k 2 ⟩/⟨k⟩ are good indicators of the critical point displaying a maximum for T = T c without noticeable finite size effects for the sizes investigated in this work (see Figure 12).We note that origin of such power laws is not necessarily related to the presence of a critical point, as the data sets in Figure 11 are representative of regimes where the correlation length is not much larger than L (and, for T = 2.5, is smaller).We attribute these features of P (k) to the strong correlations present in the system (so: not critical behavior, but rather, just correlation length much larger than the lattice spacing).

Degree correlations
In order to go beyond single node statistics and to explore how far IsingNets are from random networks with a given degree distribution, in this section we characterize the degree correlations [2,3] of the IsingNets.Degreedegree correlations measure to what extent the degree of two end nodes of the same link are correlated.Degree correlations are typically classified as either assortative or disassortative.Assortative degree correlations imply that highly connected nodes are more likely to be connected to highly connected nodes while low degree nodes are more likely to be connected to low degree nodes than in a maximally random network with the same degree distribution.Conversely, disassortative networks are net-works in which highly connected nodes are more likely to be connected to low degree nodes than in the maximally random network with the same degree distribution.Examples of assortative networks are social networks while examples of disassortative networks include the Internet and the protein interaction networks.The degree-degree correlations can be quantified by considering the average degree of the neighbour of a node k nn (i) defined as When k nn (i) tends to be higher for nodes of higher degree k i the network is classified as assortative.Instead when k nn (i) is typically lower for nodes of higher degree k i then the network is classified as disassortative.The IsingNets are clearly displaying a disassortative behavior for T < T c that deviates strongly from the behavior of the null model in which the distance matrix d is reshuffled (see Figure 13).On the contrary, for T > T c the trend of k nn versus k does not appear to be fully monotonic, while nodes of larger degrees remain more likely to connect to low degree nodes.
Interestingly the degree correlations are also affecting the average clustering coefficient C(k) [3,83,84] of nodes of degree k which display a decay as a function of k typical of networks with disassortative degree correlations (see Figure 13).
Degree-degree correlations can also be detected by the Pearson correlation coefficient r [2] between degrees of linked nodes which is plotted as a function of the temperature in Fig. 14(a) showing a clear negative (disassortative) correlations for low temperatures which strongly deviates from the null model.In Fig. 14(b) we also report the average clustering coefficient C as a function of the temperature showing that IsingNets display a much larger average clustering coefficient than the null model counterpart and that this average clustering coefficient is higher deep in the ferromagnetic phase (lower temperatures).

K-core structure
Networks can be decomposed in nested K-cores characterizing their core-periphery structure [85][86][87].A Kcore is a subgraph of the network formed by a set of M (K) nodes each having at least K connections with the other nodes of the set.Power-law networks with exponent γ ∈ (2, 3] display a significant K-core structure with the maximum K diverging with the network size and a power-law decay of M (K) as a function of K. On the contrary Erdös and Renyi networks have a finite number of K-cores also when the average degree diverges.Here we show that IsingNets have a very rich K-core structure having statistical properties that change below and above the critical temperature (see Figure 15).Indeed above the critical temperature, we observe a behavior similar to the expected behavior for sparse scale-free networks with power-law exponent between two and three presenting a broad (seemingly power-law straight line on a log-log plot) decay of M (K) versus K.However, below the critical temperature where the average degree diverges, the K-cores include more nodes while the decay of M (K) versus K is better approximated by an exponential (straight line in a log-linear plot) rather than by a power-law.

Weight-topology correlations
Interestingly in weighted networks not only the network topology can reveal relevant degree correlations showing that the networks deviate from maximally random networks, but also the weights can be distributed in a non-random way.In particular, there are two main network analyses that are able to detect weight-topology correlations.The first analysis [81] involves studying the normalized strength s/k versus the degree k for each node of the network.If the weights are distributed randomly and independently on the degree of the two end nodes there should not be any significant dependence of s/k with k.Conversely if s/k increases with k it implies that nodes with higher degrees are incident in average to links with larger weights.The second analysis [88] investigates the weight-topology correlations aiming at revealing the weight heterogeneity among links that connect to the same node.This heterogeneity if any can be quantified by calculating the inverse participation ratio Y for the weights of the links ending to node i, defined as If the weights w ij of the links (i, j) incident to node i are homogeneous, Y i ∼ 1/k i .If the weights are highly heterogeneous, Y −1 i indicates the effective number of links with significant weight.Interestingly when we measure Y i for the IsingNets we observe that this second type of weight heterogeneity is missing in the data and that Y i ∼ 1/k i indicating that for each node i the weights of the links incident to it have all weights of comparable order or magnitude (see figure 16).

Spectral properties of Ising networks
The IsingNets do not only have very interesting combinatorial and statistical properties encoded in their highly correlated structure but display also relevant geometrical properties reflected in their spectrum.In particular, the critical IsingNets display non-trivial spectral properties characterized by a power-law scaling close to criticality and a highly degenerate spectral gap.The spectral properties of networks are usually probed by considering the spectrum of the graph Laplacian describing diffusion processes.The graph Laplacian ∆ is defined as ∆ = D − A where D is the diagonal matrix whose diagonal elements are the degrees of the nodes (i.e.D ii = k i ), and A is the adjacency matrix of the network.The graph Laplacian ∆ is semi-definite positive and the spectrum always includes a zero eigenvalue with degeneracy given by the number of connected components of the network, i.e. given by the Betti number β 0 .The smallest non-zero eigenvalue of the graph Laplacian of a network is also called the Fiedler eigenvalue and is typically indicated as λ 2 (as it is the second smallest eigenvalue in a connected network).In the literature, often one distinguishes between network models displaying a finite Fiedler eigenvalue λ 2 → λ ⋆ 2 > 0 in the limit N → ∞ and network models in which λ 2 → 0 as N → ∞.In the first case, we say that the networks display a spectral gap whereas in the latter case, we say that the "spectral gap closes".Examples of networks with finite spectral dimension are random graphs above the percolation threshold and examples of networks in which the spectral gap closes are finite dimensional lattices.
In several networks in which the spectral gap closes, it is possible to observe the spectral dimension d S [6,89].The spectral dimension d S is the dimension perceived by diffusion processes on the networks encoded in the graph Laplacian.On a lattice, the spectral dimension coincides with the Euclidean dimension d of the lattice while on general network topology, the spectral dimension can be distinct from the Hausdorff dimension of the network.Interestingly also small world networks with infinite Hausdorff dimension can have a finite spectral dimension d S ≥ 2 [6,90].The spectral dimension determines the scaling of the cumulative density of the eigenvalues ρ c (λ) of the graph Laplacian for λ ≪ 1 in networks where the spectral gap closes.In particular we have that networks with a spectral dimension d S have a cumulative distribution ρ c (λ) that obeys for λ ≪ 1 where C is a constant.In Figure 17 we show that the IsingNets display nontrivial spectral properties that have very peculiar characteristics strongly deviating from their corresponding null model.Particularly noticeable are the spectral properties of IsingNets close to the critical point where one observes the coexistence of a highly degenerate finite Fiedler eigenvalue λ 2 with a power-law scaling of the cumulative distribution for λ > λ 2 with a exponent given by d ≃ 0.78 ± 0.04 for L = 40.Above the critical temperature, the degeneracy of the Fiedler eigenvalue is reduced and one observes a nontrivial spectrum reminiscent of the scale-dependent spectral dimension discussed in Refs.[91] within the critical region that at higher temperatures converges with the spectrum of the null model.Below the critical dimension, the spectral gap remains highly degenerate while the rest of the spectrum remains broadly distributed.The spectrum of the graph Laplacian is also key to characterizing the network von Neumann entropy S V N [92][93][94] defined as The von Neumann entropy strongly departs from the von Neumann entropy of the null model for low temperatures displaying a local maximum for T = T c (see Figure 18).An interesting open question that will be addressed in the following works is the relation between these spectral properties of the IsingNets graph Laplacians and the intrinsic dimension and the entropy measures that have been recently proposed starting from the unsupervised PCA analysis of spin systems [12,63].We have addressed the characterization of IsingNets following two different approaches.In the first approach, we have studied the structure of IsingNets as the filtration parameter r is increased, which enforces an effective percolation process in which nodes are subsequently aggregated by considering connected pairs of nodes at increasing distances.This percolation process reveals the presence of two giant components in IsingNets in the ferromagnetic state, each one corresponding to configurations with different magnetization.The same filtration scheme has also been used here to study the topology of the data by constructing the clique complexes of the IsingNets and calculating their persistent diagram.Interestingly, the persistent diagrams reveal that real IsingNets are formed by compact clusters as the Betti numbers of their clique complex are strongly suppressed with respect to the clique complex of the corresponding randomized null models.This network analysis conducted across the filtration is also enriched by effective visualization of the network embedding conducted using MST and the UMAP embedding and by the statistical characterization of the distribution of closeness centralities.
Secondly, our investigation of IsingNets has been conducted by considering only links at a distance less than a threshold value r taken to be the average distance of the 5th nearest neighbours nodes according to the (fully con-nected) distance matrix.These networks display a broad degree and strength distribution but their complexity extends well beyond the degree and strength distribution because IsingNets have strong degree-degree correlations and weights-degree correlations, and a rich core structure that changes significantly across the phase transition reflecting the highly nontrivial structure of the spin system that they describe.This work opens new perspectives for the unsupervised characterization of the study of phases of matter using the tools of network science.This work can be extended in different directions.The analysis performed here for the 2D Ising model can be extended to the study of other classical critical phenomena with the goal of characterizing the possible presence or the lack of universalities among the networks constructed from spin system configuration snapshots.Similarly, the same toolbox can be utilized to attack out-of-equilibrium critical behavior.Widening the class of models where our analysis can lead to physical insight is fundamental to establish ground for possible analytical treatments.In particular, before this is done, a key point to be understood is the relation between correlation length and sampling in models in different universality classes.The reason is that a naive connection (number of samples is proportional to the correlation length resolved, modulo dimensional factors) is likely not correct, based on earlier manifold learning characterization of path integrals [95] (which did not observe such relation, at least in simple quantities such as the intrinsic dimension).
Another natural extension is to consider path integrals of quantum systems [95].While the data structures of such objects might feature anisotropies due to the different roles of space and imaginary time correlations, they shall be equally amenable to the analysis discussed here.Moreover, single-sliced path integrals can also be represented as networks, as discussed in Ref. [56]: this last route provides a very promising venue for future investigation is the application of the proposed network science tools to study directly experimental data of many-body wave function snapshots.
) where sp indicates the size of the p-th largest component.The abrupt increase of Y at large r further indicates the merging of two giant components (Figure 1 (g) and Figure 2 (g)).

FIG. 1 .
FIG. 1.The percolation properties of the IsingNet (left panels) generated from 2D Ising model Monte Carlo simulations on spin systems of linear size L = 40 at temperature T = 2.12 < Tc are shown as a function of filtration parameter r.Nodes are connected if their distance is less than r.Five quantities are measured: the fraction of nodes in the largest connected component (the first row, blue line) and the fraction of nodes in the second largest connected component (the first row, orange line), the average size of components that are smaller than the second largest component ⟨s⟩ (the second row), the number of components ns (the third row) and the inverse participation ratio Y (the fourth row).The results are compared with these quantities obtained from corresponding percolation properties obtained from a randomly permuted distance matrix (right panels).The number of nodes of the IsingNets is N = 6000.

r
FIG.4.A schematic illustration of the filtration process.The barcodes are used to show the appearance and disappearance of topological features corresponding to different homology classes as the filtration parameter r is increased.The filtration process ends at r = rmax when all N nodes are fully connected and forms a N -simplex.Simplices of different dimensions are indicated by different colors.

2 ∞ 1 / 2 ( 3 )
FIG. 5.The persistent diagram corresponding to homology classes in H0, H1, H2, and H3 of the IsingNet clique complexes are plotted as a function of the filtration parameter r.Panels (a), (b), and (c) show the persistent diagrams of IsingNets obtained from the spin system of linear size L = 40 at T = 2.12 (a), T = 2.25 (b), and T = 2.50 (c).Panels (d), (e), and (f) show the persistent diagram of corresponding randomized null models obtained at T = 2.12 (d), T = 2.25 (e), and T = 2.50 (f).The networks are formed by N = 100 nodes.

FIG. 6 .
FIG. 6.The Betti numbers β0, β1, β2, and β3 of the IsingNet clique complexes are plotted as a function of the filtration parameter r.Panels (a), (b), and (c) show the persistent diagrams of IsingNets obtained from the spin system of linear size L = 40 at T = 2.12 (a), T = 2.25 (b), and T = 2.50 (c).Panels (d), (e), and (f) show the persistent diagram of corresponding randomized null models obtained at T = 2.12 (d), T = 2.25 (e), and T = 2.50 (f).The networks are formed by N = 100 nodes.

FIG. 7 .
FIG. 7. The Betti distance d β (panel (a) and (b)) and the Wasserstein distance dW (panel (c) and (d)) between the persistent diagrams of the IsingNets and its randomized null models are plotted as a function of the temperature T for homology H0 and H1.The plot shows the finite size scaling of the distances on systems of size L = 32, L = 40, and L = 48.The dashed line indicates the critical temperature Tc.The IsingNets on which the persistent diagrams have been calculated have N = 300 nodes.

FIG. 10 .
FIG. 10.The mean ⟨Cl⟩, standard deviation σ(Cl), skewness Sk(Cl), and kurtosis Ku(Cl) of closeness distribution at different temperatures T and sizes L are plotted on the IsingNets where all the distance between each pair of nodes is retained.The closeness is calculated via Eq. 5.The corresponding random networks are formed by randomly permuting the distances between node pairs.The average closeness centrality of the IsingNets coincides by construction with the average of the null model (not shown) however the higher-order statistics strongly depart from the null model behavior.The vertical green dashed lines indicate the critical temperature Tc the horizontal black dashed lines indicate Sk(Cl) = 0 in panel (c) and Ku(Cl) = 0 in panel (d).

FIG. 11 .FIG. 12 .
FIG. 11.Degree distribution P (k) and strength distribution P (s/s0) of IsingNets obtained from the spin system of linear size L = 40 formed by N = 10 4 nodes at temperature T = 2.12 and T = 2.50.Panel (a) shows the degree distribution of IsingNet at T = 2.12.Panel (b) shows the strength distribution of the same network as panel (a).Panel (c) is the same as panel (a) but obtained at T = 2.50.Panel (d) is the same as panel (b) but obtained at T = 2.50.

FIG. 13 .
FIG. 13.Average nearest neighbour degree knn and clustering coefficient C(k) on IsingNets obtained from Monte Carlo simulations of the spin system of linear size L = 40 and corresponding null models with N = 2000 nodes.The threshold distance of connecting two nodes is the 5th nearest neighbour average distance.The random network is obtained by randomly permuting distances between node pairs and nodes are connected with the same threshold.The average neighbour degree knn is shown as a function of degree k on IsingNets and corresponding null models.Panel (a) and (b) show knn obtained at T = 2.12 (a) and T = 2.35 (b).Panel (c) and (d) show knn of the corresponding null model at T = 2.12 (c) and T = 2.35 (d).Panel (e) and (f) show C(k) obtained at T = 2.12 (e) and T = 2.35 (f).Panel (g) and (h) show C(k) of the corresponding null model at T = 2.12 (g) and T = 2.35 (h).

FIG. 14 .FIG. 15 .
FIG. 14. Pearson correlation coefficient r (a) and average clustering coefficient C (b) on IsingNets and corresponding null models formed at different temperatures.The Pearson correlation coefficient r is calculated on networks with 10 4 nodes and the average clustering coefficient C is calculated for IsingNet with N = 2000 nodes coming from Monte Carlo simulation with linear size L = 40.The threshold distance of connecting two nodes is the 5th nearest neighbour average distance.The random network is obtained by randomly permuting distances between node pairs and nodes are connected with the same threshold.

FIG. 16 .
FIG. 16.The ratio of strength and degree and the inverse participation ratio on IsingNets obtained from Monte Carlo simulations of the spin system of linear size L = 40 are shown versus degree on networks formed at different temperatures.The networks are formed by N = 5000 samples.The threshold distance of connecting two nodes is the 5th nearest neighbour average distance.The random networks are obtained by randomly permuting distances between node pairs and nodes are connected with the same threshold.The first row shows the strength degree ratio s/k versus degree k on IsingNets and the second row shows which on corresponding random networks.The third row shows the inverse participation ratio Y (k) versus degree k.The networks are formed by simulation obtained at temperature T = 2.12 (left column), T = 2.35 (right column).

FIG. 17 .
FIG. 17.The cumulative distribution ρc(λ) of the eigenvalues λ of graph Laplacian ∆ of IsingNets and corresponding null models at different temperatures T .The distribution is shown at T = 2.12 (a), T = 2.27 (b), T = 2.50 (c), and T = 3.50 (d).Data are shown for Isingnets generated from Monte Carlo simulations of spin system of linear size L = 40.In panel (b), the dashed lines indicate a power-law growth shown in Eq. 13 with exponent d = 0.78 ± 0.04.

FIG. 18 .
FIG. 18.The von Neumann entropy SV N of different system size L is shown versus the temperature T .The IsingNets are formed by N = 10 4 nodes.The dashed line indicates the critical temperature Tc.