Dilution of Ferromagnets via a Random Graph-based Strategy

The dynamics and behavior of ferromagnets have a great relevance even beyond the domain of statistical physics. In this work, we propose a Monte Carlo method, based on random graphs, for modeling their dilution. In particular, we focus on ferromagnets with dimension $D \ge 4$, which can be approximated by the Curie-Weiss model. Since the latter has as graphic counterpart a complete graph, a dilution can be in this case viewed as a pruning process. Hence, in order to exploit this mapping, the proposed strategy uses a modified version of the Erd\H{o}s-Renyi graph model. In doing so, we are able both to simulate a continuous dilution, and to realize diluted ferromagnets in one step. The proposed strategy is studied by means of numerical simulations, aimed to analyze main properties and equilibria of the resulting diluted ferromagnets. To conclude, we also provide a brief description of further applications of our strategy in the field of complex networks.


I. INTRODUCTION
The study of diluted ferromagnets [1-6] dates back to several years ago, following two main paths sometimes overlapping, i.e. the statistical mechanics approach to lattices, and the graph theory approach to networks [7,8]. A notable result, coming from their combination, is the modern network theory [9][10][11]. In particular, the latter extends the classical graph theory to the analysis of networks characterized by non-trivial topologies and containing a big amount of nodes. So, the role of statistical mechanics is to offer methods and strategies for investigating the properties and the dynamics of these 'complex networks' [12,13]. Usually, investigations on ferromagnets are performed using the Ising model [14], mainly because the latter constitutes a simple and powerful tool for studying phase transitions and further applications, also beyond the domain of statistical mechanics (e.g. Data Science [15] and Machine Learning [16,17]). Despite its simplicity, the Ising model becomes, itself, a very hard problem (not yet solved) when studied in dimensions greater than 3. In those cases, the Curie-Weiss [18,19] model allows to approximate its behavior, with the advantage to be also analytically tractable (i.e. it can be exactly solved for any size of system). As result, in some conditions, solving the Ising model might require to perform numerical simulations using Monte Carlo methods [20]. For instance, the Metropolis algorithm [21] constitutes one of the early, and most adopted, strategies for simulating thermalization processes over a lattice. This latter algorithm is based on the optimization of the Hamiltonian function representing the energy of the system. Notably, the Hamiltonian of the Ising model reads where the summation is extended to all the nearest neighbors (i, j) in the lattice (realized with periodic boundary conditions, so actually becoming, in topological terms, a toroid).
As result, the value of the Hamiltonian 1 depends on the set s, i.e. the configuration of spins σ in the lattice. Accordingly, the two ground states of the system correspond to the spin configurationsŝ + = [+1, +1, . . . , +1] andŝ − = [−1, −1, . . . , −1]. Therefore, considering a lattice with N sites, and starting with a random configuration s x ∈ S, defined as s x = [σ x 1 , ..., σ x N ], the Metropolis algorithm leads the system towards a state of equilibrium which, for a temperature T = 0, corresponds to one of the two ground states. This algorithm is based on two simple steps 1. Randomly select a site i, and compute the local ∆E associated to its spin flip 2. IF (∆E ≤ 0): accept the flip; ELSE: accept the flip with probability e −∆E kT repeated until the equilibrium state is reached. We remind that k and T , appearing in the probability shown in the step (2) of the Metropolis algorithm, refer to the Boltzmann constant and to the system temperature, respectively. In addition, the term 'local' ∆E, used in the step (1), indicates that the difference in energy is computed considering only the site i and its nearest neighbors. Thus, in principle, some flips may increase the global energy of the whole system. In general, the process simulated by the Metropolis algorithm takes into account the fact that the ferromagnetic interactions J are quenched, i.e. the thermalization is fast enough to allow to consider the interactions as constant. In the opposite case, i.e. with non-constant interactions, we have different scenarios. For instance, a spin system can become glassy by introducing anti-ferromagnetic interactions (i.e. J = −1), or can undergo a dilution process by removing interactions (i.e. setting J = 0). In this work, we focus on dilution of ferromagnets introducing a strategy, based on the Erdős-Renyi model [22], for modeling this process. It is worth to recall that previous investigations (e.g. [23-26]) highlighted the critical behavior of diluted ferromagnets, including for example the ergodicity breaking and the vanishing of a giant component. So, beyond providing a novel method for dilution, we give also a description of some statistical properties of the resulting system, of the dynamical processes living on it, and on potential applications. To this end, the analyses are performed in two different conditions: for introducing the dilution strategy and studying some properties of the ferromagnets, the spin variables (i.e. σ) are considered as quenched, while for studying thermalization processes after a dilution, the quenched variables are the interactions J. Finally, the proposed strategy and the related analyses are performed by means of numerical simulations. Beyond describing the behavior of our model, we emphasize that the achieved results allow also to envision potential applications in the area of complex networks. The reminder of the paper is organized as follows: Section II introduces the proposed strategy. Section III shows results of numerical simulations. Eventually, Section IV provides a description of the main findings.

II. MODELING DILUTION ON FERROMAGNETS
Let us consider ferromagnets of dimension D ≥ 4, modeled via the Curie-Weiss (CW hereinafter) model. The latter is composed of N sites, with a position i and a spin σ ± 1.
Here, the interactions are not limited to the nearest neighbors (like in the Ising model), but are extended to all the system, i.e. every site interacts with all the others. Accordingly, the Hamiltonian of the CW model reads anti-ferromagnetic), 0 (i.e. removal). Thus, a Metropolis-like algorithm devised for flipping interactions may, in principle, generate a spin glass [27][28][29] (flipping J from +1 to −1), and perform a dilution (flipping J from +1 to 0). In addition, both processes (i.e. from +1 to −1, and to 0) can be combined, modeling the emergence of a diluted spin glass. Hence, focusing on dilution, from now on, we consider only the case J = +1 → J = 0. In doing so, starting with a random distribution of spins, a Metropolis-like algorithm (M-L hereinafter) can be defined as follows: 1. Randomly select an interaction J between two sites, and compute the local ∆E asso- As in the thermalization processes, the M-L strategy depends on the Hamiltonian of the system. Furthermore, one might consider also flipping of J from 0 to +1, i.e. modeling a kind of (edge) re-population. However, since the addition of interactions between inverse spins would increase the Hamiltonian, the actual realization of flipping 0 → +1 would be quite rare.

B. Dilution via a Random Graphs-based Strategy
As mentioned above, modern network theory and its methods are spreading in many other scientific fields. Then, it is interesting to see whether and how network theory can be useful for facing the problem of diluting ferromagnets. Notably, our work, beyond to introduce a further method for this task, allows also to prove the effectiveness of network theory in a further application. As result, the proposed model has a double valence, i.e. the process of ferromagnet dilution can be analyzed by the tools developed in network theory, and allows to envision new applications. For instance, as shown later, the subfield of community analysis can benefit from the proposed strategy. Given this premise, we can now proceed with a brief the description of ferromagnets with the formal language of graph theory. In general, a graph G is an entity composed of two sets: N (i.e. nodes) and L (i.e. edges). As above reported, the maximum number of edges (i.e. L M ) depends on N . In addition, the edges can be provided with some properties, as a direction, a weight, and so on, in order to represent specific characteristics of the object their refer to (e.g. a ferromagnet, or a real network as a social network [30], a biological network [31], a immune network [32,33], a financial network [34], and many others). In the proposed model, edges have no particular properties (i.e. they are indirect and unweighted), and the graph is implemented via the E-R model. The latter is realized by defining a number of nodes N and a parameter β, which represents the probability of each edge to exist. Thus, the expected number of edges in an E-R graph is equal to E(L) = L M · β. Notably, decreasing (increasing) β entails to remove (add) edges in the graph. The algorithm for generating an E-R graph is very simple: 1. Define the number of N of nodes and the probability β lows: i ) the ER-L strategy starts with non-connected nodes then populates the graphs with new edges, while the M-L strategy starts with a complete graph and then removes the edges; ii ) the ER-L strategy allows to obtain more configurations than the M-L strategy, being the latter 'Hamiltonian-dependent'. In particular, once the Hamiltonian has been optimized, further actions (i.e. edge removal) have very low probability. On the other hand, the ER-L strategy, being (partially) 'Hamiltonian-independent', allows the realization of ferromagnets with higher degree of dilution. For this reason, M-L is closer to a physical realization of a dilution than ER-L.

ER-L Strategy
We are now ready to present the ER-L strategy in detail. Firstly, the ER-L uses a parameter γ ∈ [0.0, 1.0], representing a kind of control in the dilution process. Notably, γ = 0.0 entails the process is not controlled, while γ = 1.0 entails a fully controlled process.
It is worth noting that, while the dilution of a ferromagnet does not require any control, being driven towards the optimization of a Hamiltonian, using a probabilistic model (i.e. the ER-L), whose dynamics depends only in part on the local energy, the so-called control parameter γ becomes fundamental for approaching the behavior of a physical dilution. Therefore, γ compensates for the partial energy independence of the proposed strategy. Accordingly, the edge probability β is 'corrected' as follows with F s step function and ω equal to In doing so, β = 0 when ω has a null or a negative value and, at the same time, the normalization condition (i.e. ω ≤ 1) is respected for any value of β, and of σ. Thus, varying the parameter β , we can study the Hamiltonian of the resulting diluted ferromagnet and its behavior. In few words, the parameter γ makes the proposed method closer to a physical dilution, since combines β with the contribution of the two spins involved in the interaction.
For instance, from a physical point of view, an interaction between two opposite spins must be removed with a probability higher than an interaction between two equal spins. At the same time, the resulting parameter ω can take values smaller than zero, hence it cannot be directly adopted as the probability to remove an edge. As result, we introduced a 'corrected' probability β , that takes as input any possible value of ω and has a range limited between zero and one. The degree of freedom offered by the parameter γ allows to represent dilution processes both in physical systems, as one can do also with a more classical approach (e.g. that before described), and to consider other systems, as social networks, where further properties and mechanisms can be involved in the process. In particular, in the case of social networks, dynamical processes like dilution might consider both the node similarity (e.g. the spin) and a probabilistic process mapped to the β parameter. Before illustrating the results of numerical simulations, it is important to elucidate a further aspect of our investigation. As previously reported, when studying the equilibrium configuration of a spin system, the interaction variables J are considered quenched. So, with the aim to analyze the behavior of diluted ferromagnets, the variation of J must be faster than that of σ, i.e. the latter is quenched. Now, having defined the ER-L model, we discuss how it can be used. realizations of G with the same amount of edges (i.e. higher its entropy, see also [35] for further details). So in a continuous dilution, beyond considering the effect of γ, one can be able to move from a state, say G(β t1 ), to a state G(β t2 ) without losing information about the edges existing at t1. For instance, if β t1 = 0.9 and β t2 = 0.8, the ER-L must account for the removal of a density of edges equal to 0.1, preserving the remaining structure of the graph.
Therefore, the simple generation of a first graph with β = 0.9, and then with β = 0.8, is not allowed because the two resulting graphs are not correlated. Thus, in the considered example, the continuous dilution process entails to move in the phase space of the graph by removing each edge with a probability much smaller than β t2 , in order to consider also the effect of β t1 . To generalize, given β t1 and β t2 , with β t1 > β t2 , if an edge e ij (i.e. connecting sites i and j) belonging to the graph in the state G t1 has to be confirmed in the state G t2 , one cannot simply use β t2 because, after the process, the edge e ij would be present with probability P (e ij ) = β t1 · β t2 , that is obviously smaller than β t2 . For this reason, we need to compute the factor such that, P (e ij ) = β t1 · = β t2 . In this way, at t2, each edge remains in the graph with probability = β t2 β t1 . In a similar fashion, we implement the inverse process, i.e. repopulating the graph with missing edges, from the state G tn to G tn+1 , now having β tn < β tn+1 . In particular, a new edge (always defined e i,j ) must be added to G with a probability P (e i,j ) = , coming from the following relation: (1 − β tn ) · (1 − ) = (1 − β tn+1 ), so that = 1 − βtn β tn+1 . Summarizing, while a diluted ferromagnet can be realized with a single instance of the ER-L model, a continuous dilution can be implemented as follows: 1. Generate a graph G(N, β 0 ) and define the sampling rate for the dilution, i.e. ∆β; 2. While β t > θ:

3.
Remove each edge in G with probability = βt The parameter θ represents the final edge probability β, i.e. the probability one should use for generating via the E-R model a graph similar to that resulting from the dilution process. β 0 corresponds to the starting value of β for generating the initial graph, and β t corresponds to the value of β at step t. The inverse process, i.e. the graph re-population, can be summarized as follows: 1. Generate a graph G(N, β 0 ) and define the sampling rate for the re-population, i.e. ∆β; 2. While β t < ζ:

3.
Add each new potential edge in G with probability = 1 − βt βt+∆β 4. β t = β t + ∆β As θ, ζ represents the final value of β one should use for generating a similar graph (i.e. with same statistical properties), achieved after re-population, using the E-R model. Moreover, we clarify that 'new potential edge' refers to the edges that can be added to the graph for making it again fully-connected, i.e. it refers only to missing edges. Eventually, we analyze also thermalization processes (considering, after each dilution, the variables J as quenched ).
To this end, the system magnetization defined as offers a macroscopic view on the process. Notably, we recall that the magnetization is an order-parameter and allows both to observe the emergence of a phase transition, and to evaluate its nature (e.g. first order). In addition, it is worth emphasizing that quenched spins, randomly initialized with a uniform distribution, entail the magnetization is on average always null (i.e. the system remains in a disordered phase). Further analyses devised for studying the behavior of our model are introduced in the following section.

III. RESULTS
The proposed model is studied by means of numerical simulations, considering ferromagnets composed of N = 1000 sites. In particular, we aim to obtain diluted ferromagnets with single realizations of the ER-L strategy, and to use the latter for modeling continuous dilution and re-population processes. In addition, we analyze thermalization processes on the resulting diluted ferromagnets and, eventually, we present a potential application in the field of complex networks, i.e. in the evaluation of community stability [36,37].

Dilution via the ER-L Model
We start considering different realizations of ferromagnets via the ER-L model, on varying β and γ. Figure 2 shows the (absolute value of) Hamiltonian H, normalized over the actual number of edges L a , which reads with s denoting a specific spin configuration. It is important to emphasize that eq. 6 is normalized in order to consider only those connections that survive during the dilution process. As expected, the Hamiltonian (eq. 6) is equal to zero when there is no control in the dilution process, since interactions are removed without considering the spin of related nodes. On the contrary, increasing γ, we observe that the Hamiltonian increases up to 1 -we remind that we are considering the absolute value of the Hamiltonian, so that its actual value is −1. For γ > 0.0, the maximum of |H| can be reached spanning β within well defined ranges. Notably, the latter enlarges by increasing γ. For instance, when γ = 0.5 the optimal H is obtained with 0 ≤ β ≤ 0.3, while when γ = 1.0 the optimal H is obtained with setting β = 0.5, while in the case γ = 1.0 the plot illustrates a graph achieved with β > 0.5 and one with β = 0.5. Remarkably, for β ≤ 0.5, the resulting graph appears perfectly divided between the two communities (i.e. spins +1 separated from spins −1). Instead, for values of β slightly higher than 0.5, as represented in Fig. 3, the two communities are connected by few edges. This observation is very important, because strongly related to thermalization processes (i.e. when the variables J are taken as quenched after the dilution step). Numerical simulations, shown in fig. 4, demonstrate that the ER-L strategy is able to dilute and to repopulate a graph, no matter the value of γ. In addition, implementing the two processes as a cycle, we did not find any form of hysteresis, i.e. dilution and repopulation cover two perfectly overlapping paths in the plot of fig. 4. Only in the case with γ = 0.0, we found an observable difference between the two paths, which can still be considered negligible.
Now, we study thermalization processes on ferromagnets diluted with different γ. To this end, ferromagnets can be diluted both implementing single realizations of the ER-L strategy (as we did here), and by performing the continuous dilution, i.e. considering the resulting graph obtained at each step. Moreover, we remind that thermalization is analyzed by studying the average magnetization (i.e. eq. 5) of the system -see Fig. 5. In all cases, it seems that the order-disorder phase transition occurring in the ferromagnet is of first order, no matter the value of γ. At the same time, the latter strongly affects the critical β c . In particular, for γ = 0.0 we found β c ∼ 10 −3 , while for γ = 1.0 the value is smaller than 0.5 + 10 −4 . When γ = 0, the transition is caused by the relevant reduction of edges, so that without interactions the thermalization cannot take place. On the contrary, increasing γ, the order-disorder phase transition is caused initially by a combined effect of edge reduction and community separation, until γ = 1.0, where the disordered phase is reached for the emergence of two well separated, and ordered, communities having opposite spin (i.e. one with σ = +1 and one with σ = −1) -see the inset of fig. 5. In addition, we found that with β = 0.5 + 10 −5 and γ = 1.0, the average (absolute) value of the magnetization and the variance are equal, proving its role as the critical β c . Then, it is worth to further clarify an aspect shown in fig. 5, i.e. the first order phase transition. Notably, when the dilution is strongly controlled (i.e. γ → 1) the first edges to be removed are those linking nodes with opposite spin. So, that once the half of the edges is removed, only those connecting nodes with same spins survive, leading towards a total magnetization equal to zero (i.e. summing the magnetization observed in the two separated communities, which in turn reach opposite states of full order). Instead, for poorly controlled dilution, edges are removed without to consider the value of related spins, so that the transition occurs for lower values of β. To conclude, we remind that these simulations have been performed on ferromagnets containing an equal amount of positive and negative spins.

A. Community Stability
The proposed strategy aims to perform dilution processes on ferromagnets using, as reference, a well-established random graph model (i.e. the E-R model). The latter is widely used in the modern theory of networks for studying dynamical processes, and structural properties of complex networks. Now, we want to evaluate if a modification of the E-R model that we introduced, i.e. the ER-L strategy, can be useful for extracting information from a complex network. In particular, we envision a potential application in the task of measuring the stability of a community, i.e. if according to the properties of its nodes, it risks to disappear after a while. Notably, in a number of models studied in social dynamics [38][39][40][41], often properties and behaviors are mapped onto binary spins [42]. So, in principle, one could use the Hamiltonian defined in eq. 6 for measuring the stability of a community [43], i.e. the higher its |H| the higher the probability that the community survives over time.
In particular, the value of |H| reflects the degree of similarity between the nodes connected in the same community. Even if only the analysis of real datasets would allow to confirm the validity of this hypothesis, and then also the usefulness in the area of complex networks of the proposed strategy, our assumption is based on the simple observation that groups of individuals are more likely to cluster together when share common interests, opinions, and so on. Moreover, beyond observations of real scenarios, this mechanism is confirmed by the positive assortativity [44] that social networks show, i.e. individuals are more likely to interact with their likes. In addition, recalling that the value of eq. 6 can be related to single communities, it can be viewed as an alternative form of assortativity at community level, since the higher its value the higher the fraction of connections between similar nodes.
So, since the case with binary spins has been previously studied, even if referred to dilution processes, here we focus on two main analyses. First we study the influence of heterogeneity in a complete community, measuring how the Hamiltonian decreases while increasing the amount of nodes with different spins. Second, we study the Hamiltonian of a graph, considering the XY model [45] as reference. In doing so, we are able to represent situations where there are more than 2 opinions (e.g. [46,47]), states, or behaviors. The first investigation considers, at the beginning, the combination named s + , then some spins flip increasing the density of nodes with spin −1 (this process leads to the same results also considering the inverse case, i.e. starting with s − and then flipping spins to +1). Results are illustrated in fig. 6. As expected, the minimum of |H| is reached when the number of +1 spins is equal to that of −1 spins. Finally, we analyzed the Hamiltonian of a community using the XY model. Figure 7 reports the related results, for different γ, and considering both 4 different states, and 360 different states. Here, the pairs of spins are evaluated according to the cosine similarity cos(θ a − θ b ), with θ a and θ b representing the value of the involved spins.
We observe two main differences from the classical binary spins. In particular, using the XY model, the decrease of the Hamiltonian is smoother in the XY model than in the Ising model, where it appears less monotonous. Furthermore, the maximum value of |H| is always smaller than 1.0. Accordingly, we note that communities are more stable (or robust to spin flipping) when there are more than 2 states characterizing the related nodes. On the other hand, many possible states do not allow to reach a perfect stability (i.e. |H| = 1), exposing a community to a higher risk to disappear.

IV. DISCUSSION AND CONCLUSION
This work introduces a strategy, named ER-L, for modeling the dilution of ferromagnets using the framework of modern network theory. In particular, we adopt as reference the E-R graph model, since the latter, under opportune conditions, constitutes the graphical representation of the Curie-Weiss model. The proposed method is partially Hamiltonianindependent, i.e. while a Metropolis-like strategy can dilute a ferromagnet according to energy-based rules, the ER-L strategy depends on a probabilistic (non-physical) parameter β and, only in part, on the local energy via a parameter γ, which represents a kind of control in the dilution. Notably, diluting a ferromagnet can be thought as pruning a graph G, moving the latter in a phase space composed of all its possible realizations. The amount of edges (i.e. interactions) depends on the parameter β of the ER-L model, so our strategy moves G along the axis β. In doing so, G undergoes a kind of phase transition (see also [48][49][50][51]), where different structures can be obtained, from sparse nodes to a complete graph. In addition, the parameter γ ensures that the motion along the β axis, on the phase space, corresponds to that followed by a ferromagnet during a spontaneous dilution. In particular, like during thermalization, a system tends to naturally reach an equilibrium state that minimizes its energy. In a similar fashion, spontaneous dilutions should lead the system towards a ground state. Results indicate that ER-L is able to perform this task for different values of γ, depending on the considered β. In addition, we also analyzed the closed path (i.e. the cycle) from a complete graph to single nodes, and then to a complete graph by repopulating with new edges the diluted graph. It is worth to highlight that ER-L allows to dilute a graph also after its Hamiltonian has been minimized, while a Metropolis-like strategy, being 'Hamiltonian-dependent', would not be able. Therefore, not all the structures obtained via ER-L have a physical meaning. However, during a continuous dilution, we can discriminate those that appear without a physical meaning, computing the difference of the Hamiltonian between the two structures. The related analyses have been performed considering the spin variables σ as quenched. So, we studied also the opposite case, i.e. after a dilution the interactions J become quenched, and the spins can flip towards an equilibrium state (see also [52][53][54]). In order to study this process, and to make a relation with the parameter β, we analyze the average magnetization achieved at equilibrium, which provides an indication about the phase transitions occurring in the system -see fig. 5. Once analyzed the outcomes of the ER-L strategy, we performed a further analysis for evaluating the opportunity to apply it to further tasks, in particular considering the measure of stability of communities in complex networks. First we analyzed the variation of the Hamiltonian turning an ordered system to a disordered one. Then, we studied the Hamiltonian considering as reference the XY model, i.e. admitting spins with more than 2 values. This preliminary investigation suggests that ER-L may, in principle, be useful for evaluating the risk that a community will dissolve after a while, according to the degree of heterogeneity of its individuals (e.g. in terms of opinions, interests, and so on). In addition, we found that communities whose nodes have more possible states (e.g. opinions), never reach a perfect stability (i.e. |H| < 1) but can be more robust than those with binary spins to the emergence of interaction between different individuals. Obviously, we are not taking into account all the 'social' processes that may occur in real systems, e.g. once two different individuals interact, one might imitate the other, relaxing the system. Moreover, considering the two strategies here described, i.e.
ER-L and M-L, we deem important to mention that, in principle, they might constitute also the base for developing learning algorithms [55]. Notably, almost all simulations have been carried on considering an equal distribution of positive and negative spins, however different combinations (i.e. patterns) might be used. Therefore, the optimization of the Hamiltonian during the dilution, in our view, even if referred to only one pattern, can be actually interpreted as a form of learning in a neural network [56]. On the other hand, further investigations are required for evaluating whether the proposed model may allow the graph to learn and store more than one pattern. Finally, we remark that, in order to assess the actual usefulness of ER-L for evaluating the community stability, further investigations based on real datasets are definitely mandatory.
The authors declare that there is no conflict of interest regarding the publication of this paper.

ACKNOWLEDGMENTS
MAJ would like to acknowledge support by the H2020-645141 WiMUST project, and to thank the mobility funds of the Faculty of Psychology and Educational Sciences of Ghent University. Authors wish to thank Adriano Barra for the priceless suggestions.