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 , 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ő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, aiming 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.

1. Introduction

The study of diluted ferromagnets [16] dates back to several years ago, following two main paths sometimes overlapping, that is, 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 [911]. In particular, the latter extends the classical graph theory to the analysis of networks characterized by nontrivial 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 approximating its behavior, with the advantage of being also analytically tractable (i.e., it can be exactly solved for any size of system). As a result, in some conditions, solving the Ising model might require performing 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 readswhere the summation is extended to all the nearest neighbors in the lattice (realized with periodic boundary conditions, so actually becoming, in topological terms, a toroid). As a result, the value of Hamiltonian 1 depends on the set , that is, the configuration of spins in the lattice. Accordingly, the two ground states of the system correspond to the spin configurations and . Therefore, considering a lattice with sites and starting with a random configuration , defined as , the Metropolis algorithm leads the system towards a state of equilibrium which, for a temperature , corresponds to one of the two ground states. This algorithm is based on two simple steps: (1)Randomly select a site , and compute the local associated with its spin flip.(2)If , accept the flip;else, accept the flip with probability .

They are repeated until the equilibrium state is reached. We recall that and , appearing in the probability shown in step of the Metropolis algorithm, refer to the Boltzmann constant and to the system temperature, respectively. In addition, the term “local” , used in step , indicates that the difference in energy is computed considering only the site 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 are quenched, that is, the thermalization is fast enough to allow considering the interactions as constant. In the opposite case, that is, with nonconstant interactions, we have different scenarios. For instance, a spin system can become glassy by introducing antiferromagnetic interactions (i.e., ) or can undergo a dilution process by removing interactions (i.e., setting ). 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 recalling that previous investigations (e.g., [2326]) 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 quenched, while for studying thermalization processes after a dilution, the quenched variables are the interactions . 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 envisioning potential applications in the area of complex networks. The reminder of the paper is organized as follows: Section 2 introduces the proposed strategy. Section 3 shows results of numerical simulations. Eventually, Section 4 provides a description of the main findings.

2. Modeling Dilution on Ferromagnets

Let us consider ferromagnets of dimension , modeled via the Curie-Weiss (CW hereinafter) model. The latter is composed of sites, with a position and a spin . Here, the interactions are not limited to the nearest neighbors (like in the Ising model) but are extended to all the system; that is, every site interacts with all the others. Accordingly, the Hamiltonian of the CW model readswith , that is, combination of spins . Then, like in the Ising model, the two ground states correspond to the spin configurations and , that is, those that minimize the value of (see (2)). It is worth highlighting that the CW model can be represented as a complete (i.e., fully connected) graph, where each site corresponds to a node and each interaction corresponds to an edge. In addition, the number of interactions is equal to . The mapping from the physical object (i.e., the ferromagnet) to the mathematical entity (i.e., the graph) allows mapping a dilution to a pruning process. However, before presenting the dilution strategy developed with the framework of graph theory, we discuss the application of a more classical method, that is, the previously mentioned Metropolis algorithm.

2.1. Dilution by the Metropolis Algorithm

In principle, the Metropolis algorithm, and similar methods, may be used for modeling the dilution of ferromagnets. Notably, since this algorithm modifies spins from to , and vice versa, according to the energy difference resulting from the spin flipping, an opportune variant—say Metropolis-like—might be used for flipping the interaction variables . In the case of spin flipping the possible values that can take are whereas, in the case of interactions , the latter may take three different values: (i.e., ferromagnetic), (i.e., antiferromagnetic), (i.e., removal). Thus, a Metropolis-like algorithm devised for flipping interactions may, in principle, generate a spin glass [2729] (flipping from to ) and perform a dilution (flipping from to ). In addition, both processes (i.e., from to and to ) can be combined, modeling the emergence of a diluted spin glass. Hence, focusing on dilution, from now on, we consider only the case . 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 between two sites, and compute the local associated with its flip to .(2)If , accept the flip;else, accept the flip with probability .

As in the thermalization processes, the M-L strategy depends on the Hamiltonian of the system. Furthermore, one might consider also flipping of from to , that is, modeling a kind of (edge) repopulation. However, since the addition of interactions between inverse spins would increase the Hamiltonian, the actual realization of flipping would be quite rare.

2.2. 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 proving the effectiveness of network theory in a further application. As a result, the proposed model has a double valence; that is, the process of ferromagnet dilution can be analyzed by the tools developed in network theory and allows envisioning 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 description of ferromagnets with the formal language of graph theory. In general, a graph is an entity composed of two sets: (i.e., nodes) and (i.e., edges). As above reported, the maximum number of edges (i.e., ) depends on . 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 they refer to (a ferromagnet or a real network as a social network [30], a biological network [31], an 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 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 . Notably, decreasing (increasing) entails removing (add) edges in the graph. The algorithm for generating an E-R graph is very simple:(1)Define the number of of nodes and the probability .(2)Draw each edge with probability .

In particular, in step (2), all possible edges are considered. Therefore, an E-R graph generated with contains exactly edges (i.e., it is complete) and constitutes the graphical counterpart of the CW model previously described. The tuning of the parameter allows representing ferromagnets with different amounts of interactions. Thus graphs generated with (i.e., having ) represent diluted ferromagnets. This last observation constitutes the base of our model, that is, an E-R-like (ER-L hereinafter) model devised for dilution processes. For the sake of clarity, now we provide a pictorial representation for highlighting their main differences between two mentioned strategies, that is, M-L and ER-L—see Figure 1. So, a quick glance to the pictorial representation allows observing what follows: (i) the ER-L strategy starts with nonconnected nodes and 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 obtaining more configurations than the M-L strategy, the latter being “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.

2.2.1. ER-L Strategy

We are now ready to present the ER-L strategy in detail. Firstly, the ER-L uses a parameter , representing a kind of control in the dilution process. Notably, entails the process is not controlled, while 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 step function and equal toIn doing so, when has a null or a negative value and, at the same time, the normalization condition (i.e., ) 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 it 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 a 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 both representing dilution processes in physical systems, as one can do also with a more classical approach (e.g., that described before), and considering 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 are considered quenched. So, with the aim of analyzing the behavior of diluted ferromagnets, the variation of must be faster than that of ; that is, the latter is quenched. Now, having defined the ER-L model, we discuss how it can be used. First, one can realize a diluted ferromagnet via ER-L, studying then its properties. Second, one can analyze the behavior of the ferromagnet during the dilution process. While a single realization, of a diluted ferromagnet, is very similar to the realization of a graph via the E-R model, that is, drawing edges according to the probability, it is worth explaining how to implement a continuous dilution process by the ER-L model. To this end, let us consider a dilution as the motion of a graph over a phase space, along an axis . The phase space is larger for low values of and becomes narrower as increases, until it contains only one state when (i.e., fully connected). The lower , the larger the number of possible realizations of with the same amount of edges (i.e., the 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 , to a state without losing information about the edges existing at . For instance, if and , the ER-L must account for the removal of a density of edges equal to , preserving the remaining structure of the graph. Therefore, the simple generation of a first graph with and then with , is not allowed because the two resulting graphs are not correlated. Thus, in the considered example, the continuous dilution process entails moving in the phase space of the graph by removing each edge with a probability much smaller than , in order to consider also the effect of . To generalize, given and , with , if an edge (i.e., connecting sites and ) belonging to the graph in the state has to be confirmed in the state , one cannot simply use because, after the process, the edge would be present with probability , that is, obviously smaller than . For this reason, we need to compute the factor such that . In this way, at , each edge remains in the graph with probability . In a similar fashion, we implement the inverse process, that is, repopulating the graph with missing edges, from the state to state , now having . In particular, a new edge (always defined ) must be added to with a probability , coming from the following relation: , so that . In summary, 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 and define the sampling rate for the dilution, that is, .(2)While .(3)_ _ Remove each edge in with probability .(4)_ _  .

The parameter represents the final edge probability , that is, the probability one should use for generating via the E-R model a graph similar to that resulting from the dilution process. corresponds to the starting value of for generating the initial graph, and corresponds to the value of at step . The inverse process, that is, the graph repopulation, can be summarized as follows:(1)Generate a graph and define the sampling rate for the repopulation, that is, .(2)While .(3)_ _ Add each new potential edge in with probability .(4)_ _  .

As , represents the final value of one should use for generating a similar graph (i.e., with same statistical properties), achieved after repopulation, 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; that is, it refers only to missing edges. Eventually, we analyze also thermalization processes (considering, after each dilution, the variables as quenched). To this end, the system magnetization defined asoffers 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 that 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.

3. Results

The proposed model is studied by means of numerical simulations, considering ferromagnets composed of 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 repopulation 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 , normalized over the actual number of edges , which readswith denoting a specific spin configuration. It is important to emphasize that (6) is normalized in order to consider only those connections that survive during the dilution process. As expected, the Hamiltonian (see (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 —we recall that we are considering the absolute value of the Hamiltonian, so that its actual value is . For , the maximum of can be reached spanning within well-defined ranges. Notably, the latter enlarges by increasing . For instance, when the optimal is obtained with , while when the optimal is obtained with . It is worth clarifying that we cannot study the dilution in function of , since its value depends on the involved spins, consequently being potentially different from edge to edge. This preliminary investigation constitutes an early indicator of the existence of a critical edge probability , that is, the highest value of to have . Now, we start analyzing continuous dilution and repopulation processes. Before to show results of simulations, let us highlight that modeling a continuous dilution entails starting the process from a fully connected graph. Therefore, considering Figure 3, here the direction of a dilution corresponds to that of the M-L strategy while, obviously, a repopulation follows the inverse path. The pictorial representation of Figure 3 aims to give an overview about the continuous dilution process (via ER-L) in two conditions: noncontrolled () and fully controlled (i.e., ). Provided that the starting graph and the landing one are equal in both cases, a quick glance to the pictorial allows appreciating the influence of . In particular, the intermediate graphs, between the starting and the ending ones, are related to those achieved via the ER-L method setting . Notably, in the case , the graph is obtained setting , while in the case the plot illustrates a graph achieved with and one with . Remarkably, for , the resulting graph appears perfectly divided between the two communities (i.e., spins of separated from spins of ). Instead, for values of slightly higher than , as represented in Figure 3, the two communities are connected by few edges. This observation is very important, because of being strongly related to thermalization processes (i.e., when the variables are taken as quenched after the dilution step). Numerical simulations, shown in Figure 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; that is, dilution and repopulation cover two perfectly overlapping paths in the plot of Figure 4. Only in the case with , we found an observable difference between the two paths, which can still be considered negligible.

Thermalization Processes on Diluted Ferromagnets. Now, we study thermalization processes on ferromagnets diluted with different . To this end, ferromagnets can be diluted both by implementing single realizations of the ER-L strategy (as we did here) and by performing the continuous dilution, that is, considering the resulting graph obtained at each step. Moreover, we recall that thermalization is analyzed by studying the average magnetization (i.e., (5)) of the system—see Figure 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 . In particular, for we found , while for the value is smaller than . When , 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 , where the disordered phase is reached for the emergence of two well separated, and ordered, communities having opposite spin (i.e., one with and one with )—see the inset of Figure 5. In addition, we found that with and the average (absolute) value of the magnetization and the variance are equal, proving its role as the critical . Then, an aspect shown in Figure 5, that is, the first order phase transition, is worth further clarifying. Notably, when the dilution is strongly controlled (i.e., ) the first edges to be removed are those linking nodes with opposite spin, so that once half of the edges are removed only those connecting nodes with same spins survive, leading to 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 considering the value of related spins, so that the transition occurs for lower values of . To conclude, we recall that these simulations have been performed on ferromagnets containing an equal amount of positive and negative spins.

3.1. 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, that is, 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, that is, a predictive analysis, based on the properties of nodes, for evaluating the risk of its disappearance. Notably, in a number of models studied in social dynamics [3841], often properties and behaviors are mapped onto binary spins [42]. So, in principle, one could use the Hamiltonian defined in (6) for measuring the stability of a community [43], that is, the higher its the higher the probability that the community survives over time. In particular, the value of reflects the degree of similarity between the nodes connected in the same community. Even if only the analysis of real datasets would allow confirming 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 sharing 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; that is, individuals are more likely to interact with their likes. In addition, recalling that the value of (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 referring 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 model [45] as reference. In doing so, we are able to represent situations where there are more than opinions (e.g., [46, 47]), states, or behaviors. The first investigation considers, at the beginning, the combination named ; then some spins flip increasing the density of nodes with spin (this process leads to the same results also considering the inverse case, that is, starting with and then flipping spins to ). Results are illustrated in Figure 6. As expected, the minimum of is reached when the number of spins is equal to that of spins. Finally, we analyzed the Hamiltonian of a community using the model. Figure 7 reports the related results, for different , considering both different states and different states. Here, the pairs of spins are evaluated according to the cosine similarity , with and representing the value of the involved spins. We observe two main differences from the classical binary spins. In particular, using the model, the decrease of the Hamiltonian is smoother in the model than in the Ising model, where it appears less monotonous. Furthermore, the maximum value of is always smaller than . Accordingly, we note that communities are more stable (or robust to spin flipping) when there are more than states characterizing the related nodes. On the other hand, many possible states do not allow reaching a perfect stability (i.e., ), exposing a community to a higher risk of disappearing.

4. 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 Hamiltonian-independent, that is, while a Metropolis-like strategy can dilute a ferromagnet according to energy-based rules, the ER-L strategy depends on a probabilistic (nonphysical) 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 , 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 along the axis . In doing so, undergoes a kind of phase transition (see also [4851]), 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 to 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 highlighting that ER-L allows diluting a graph also after its Hamiltonian has been minimized, while a Metropolis-like strategy, being “Hamiltonian-dependent,” would not be able to do so. 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; that is, after a dilution the interactions become quenched, and the spins can flip towards an equilibrium state (see also [5254]). 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 Figure 5. Once the outcomes of the ER-L strategy are analyzed, 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 model, that is, admitting spins with more than 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, etc.). In addition, we found that communities whose nodes have more possible states (e.g., opinions) never reach a perfect stability (i.e., ) 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; for example, once two different individuals interact, one might imitate the other, relaxing the system. Moreover, considering the two strategies here described, that is, ER-L and M-L, we deem that it is 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 referring 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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


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