Kinetic Ising models with various single-spin flip dynamics on quenched and annealed random regular graphs

TL;DR

This study examines how different single-spin flip dynamics affect the behavior of kinetic Ising models on quenched and annealed random regular graphs, revealing that some dynamics are more robust to disorder than others.

Contribution

It compares the effects of various single-spin flip algorithms on Ising models in quenched and annealed networks, highlighting differences in steady states and phase transition types.

Findings

On quenched graphs, all algorithms produce identical results.

Annealed graphs cause different dynamics to lead to distinct steady states.

Some algorithms are resistant to annealed disorder, showing only small quantitative changes.

Abstract

We investigate a kinetic Ising model with several single-spin flip dynamics (including Metropolis and heat-bath) on quenched and annealed random regular graphs. As expected, on the quenched structures all proposed algorithms reproduce the same results since the conditions for the detailed balance and the Boltzmann distribution in an equilibrium are satisfied. However, on the annealed graphs situation is far less clear -- the network annealing disturbs the equilibrium moving the system away from it. Consequently, distinct dynamics lead to different steady states. We show that some algorithms are more resistant to the annealed disorder, which causes only small quantitative changes in the model behavior. On the other hand, there are dynamics for which the influence of annealing on the system is significant, and qualitative changes arise like switching the type of phase transition from continuous to discontinuous one. We try to identify features of the proposed dynamics which are responsible for the above phenomenon.

Equations

$$\Delta E=E_{\text{new}}-E_{\text{old}}=2J\sum _{⟨i,j⟩}{S}_i{S}_j.$$

$$\frac{W_{\text{old}\to \text{new}}}{W_{\text{new}\to \text{old}}}=\exp (-\beta \Delta E).$$

$$W_{\text{M}}=\min [1,\exp (-\beta \Delta E)],$$

$$W_{\text{HB}}=\frac{1}{1+\exp (\beta \Delta E)}.$$

$$m=\frac{1}{M}\sum _{j=1}^M{m}_j,$$

$$m_j=\frac{1}{N}\sum _{i=1}^N{S}_i^j,$$

$$\frac{J}{T_c}=\frac{1}{2}\ln \left(\frac{⟨k^2⟩}{⟨k^2⟩-2⟨k⟩}\right).$$

$$T_c(q)=2{\left[\ln \left(\frac{q}{q-2}\right)\right]}^{-1}.$$

$$c(t)=\frac{N_{\uparrow }(t)}{N}=\frac{m(t)+1}{2}.$$

$${\gamma }^{+}(c)=\sum _{k=0}^q\left(\genfrac{}{}{0ex}{}{q}{k}\right)c^{q-k}(1-c{)}^{k+1}\frac{1}{1+e^{\frac{-2}{T}(q-2k)}},$$

$${\gamma }^{-}(c)=\sum _{k=0}^q\left(\genfrac{}{}{0ex}{}{q}{k}\right)(1-c{)}^{q-k}{c}^{k+1}\frac{1}{1+e^{\frac{-2}{T}(q-2k)}}.$$

$$c(t+1)=c(t)+\left[{\gamma }^{+}(c)-{\gamma }^{-}(c)\right],$$

$${\gamma }^{+}(c)-{\gamma }^{-}(c)=0.$$

$$F(c,T)={\gamma }^{+}(c,T)-{\gamma }^{-}(c,T).$$

$$F(c,T)=0.$$

$$\frac{\partial F(c,T)}{\partial c}+\frac{\partial F(c,T)}{\partial T}\frac{dT}{dc}=0.$$

$${\frac{dT}{dc}|}_{(\frac{1}{2},T_c)}=0.$$

$${\frac{\partial F(c,T)}{\partial c}|}_{(\frac{1}{2},T_c)}=0.$$

$$\sum _{i=0}^q\left(\genfrac{}{}{0ex}{}{q}{i}\right)(2i-q)\tanh \frac{2i-q}{T_c}=2^q.$$

$$W_{{\text{M}}^{\prime }}=\{\begin{array}{ll}\exp (-\beta \Delta E)& (\Delta E>0)\\ W_0& (\Delta E=0)\\ 1& (\Delta E<0).\end{array}$$

$$W_{\exp }=C\exp \left(-\frac{1}{2}\beta \Delta E\right),$$

$$W_{\tanh }=\min \left[1,\exp (-\beta \Delta E)\tanh \left|a\beta \Delta E\right|\right],$$

$$V(m,T)=-\int F(m,T)dm.$$

$$V(m)=A{m}^2+B{m}^4+C{m}^6+O(m^8),$$

$$A$$

$$B$$

$$C$$

Full text

I Introduction

Recently a puzzling phenomenon related to the phase transitions has been found in the modified kinetic Ising model, the $q$-neighbor Ising model, with the Metropolis algorithm on the complete graph Jed:Chm:Szn:15 . Within this model, each spin interacts only with $q$ spins randomly chosen from its neighborhood. It has occurred that only for $q=3$ there is a continuous order-disorder phase transition, whereas for $q\geq 4$ discontinuous phase transition appears. Moreover, for $q\geq 4$ the hysteresis exhibits oscillatory behavior, i.e., expands for the even values of $q$ and shrinks for the odd values of $q$. Despite the fact that most of the results presented in Jed:Chm:Szn:15 were derived from the master equation and not only from the computer simulations, the intuitive understanding of the non-monotonicity of hysteresis is still missing.

Above results are even more confusing if we realize that the $q$-neighbor Ising model on the complete graph seems to be identical with the classical Ising model with the nearest-neighbors (nn) interactions on the annealed random $q$-regular graph ($q$-RRG) Lee:09 ; Bol:01 ; Par:Hoh:17 . The confusion comes from the fact that the Ising model with ferromagnetic nn interactions without an external field has been already examined on a number of different quenched and annealed networks, and it always manifests continuous phase transitions Pek:01 ; Agata:02 ; Dor:Gol:Men:02 ; Herr:15 ; Lee:09 ; Mal:Val:Das:14 ; Lip:Gon:Lip:15 . Discontinuous or mixed-order phase transitions have been observed only if the range of interactions were infinite Bar:Muk:14a ; Bar:Muk:14b or on the coupled networks Such:Hol:09 . Extensive studies have shown not only that phase transition in the Ising model is continuous but also characterized with the same set of bulk critical exponents on quenched and annealed networks Lee:09 .

Only very recently it has been shown that the $q$-neighbor Ising model is a limiting case of a non-equilibrium system with two heat baths: the Ising spins are in thermal contact with the heat bath $B_{S}$ with temperature $T_{s}$, whereas links are in thermal contact with another bath $B_{L}$ with temperature $T_{L}$ Par:Hoh:17 . For $T_{L}=\infty$ such a generalized model reduces to the $q$-neighbor Ising model investigated in Jed:Chm:Szn:15 . It has been also shown that for $T_{L}=\infty$ there is nonzero positive heat flux, which confirms that the $q$-neighbor Ising model is indeed out of equilibrium.

Paper Par:Hoh:17 explains the puzzle of seemingly contradictory results – continuous phase transition in the equilibrium Ising model on the annealed network Lee:09 vs. discontinuous phase transition in the $q$-neighbor Ising model Jed:Chm:Szn:15 . However, many interesting questions are still open. In particular it is still unclear if we can fully explain the switch from continuous to discontinuous phase transition on the basis of non-equilibrium regime. Is it possible that other dynamics, even for $T_{L}=\infty$ (out-of-equilibrium), would lead to the same results as the equilibrium model? It seems quite probable, because results presented in Lip:Gon:Lip:15 suggest that for the heat-bath dynamics the transition would be continuous, at least on the annealed directed graph. Summarizing, in this paper we ask the question about the role of dynamics in the non-equilibrium kinetic Ising model. Referring to this question we examine the kinetic Ising model with several dynamics (including Metropolis and heat-bath) on quenched and annealed $q$-RRGs. Precisely, we study several different algorithms (dynamics) that belong to a broad class of a single-spin-flip dynamics (i.e., Glauber dynamics, in a broad sense) Glauber ; God:Luc:05 . Within such dynamics each spin is flipped ($S_{i}\rightarrow-S_{i}$) with a rate $W(\Delta E)$ per unit time, and this rate is assumed to depend only on the energy difference implied in the flip:

[TABLE]

We investigate only such transition probabilities that on the quenched graphs fulfill detailed balance and lead to the Boltzmann distribution in a stationary state, which means that they satisfy the following general condition:

[TABLE]

This is quite obvious that on the quenched graph they should all lead to the same result, which will be checked in the next section. However, on the annealed graph the situation is far from being obvious, and we will show that the results strongly depend on the dynamics. The most popular choices for $W(\Delta E)$ are Metropolis Met:etal:53 and heat-bath dynamics Jan:08 . For the Metropolis transition probability:

[TABLE]

and for the heat-bath:

[TABLE]

However, there are many other possibilities that fulfill the required condition given by Eq. (2), and several examples will be investigated later in this paper.

II The model

In the previous paper Jed:Chm:Szn:15 we have investigated the $q$-neighbor Ising model, which was basically the Ising model with the Metropolis dynamics Met:etal:53 with seemingly small modification that was inspired by the $q$-voter model Cas:Mun:Pas:09 ; Nyc:Szn:Cis:12 ; Mob:15 ; Mel:Mob:Zia:17 . Within the $q$-neighbor Ising model we consider a graph that consists of $N$ nodes, and each node is occupied by exactly one spin, described by a dynamical binary variable $S_{i}=\pm 1,i=1,...,N$. In each elementary update a single, picked at random, spin $S_{i}$ interacts with its $q$ neighbors, also randomly chosen from its entire neighborhood. Generally, two methods of selecting neighbors are possible – with or without repetitions. For example, in the original $q$-voter model repetitions were possible, i.e., a given neighbor could be selected more than once Cas:Mun:Pas:09 . However, later also the $q$-voter model without repetitions has been investigated Nyc:Szn:Cis:12 . The version without repetitions is of course suitable only for the selected topologies in which the minimal degree of a node $k_{min}\geq q$. For the complete graph both methods are applicable for arbitrary value of $q<N$ and give roughly the same results. Here, we continue the study on the $q$-neighbor Ising model started in Jed:Chm:Szn:15 ; therefore, we use the same method as previously, i.e., without repetitions.

It is worth to notice that for the version without repetitions, that we use here, instead of defining the $q$-neighbor Ising model on the complete graph, one could think about the Ising model on the annealed random $q$-regular graph ($q$-RRG) Bol:01 . Degree distribution for such a graph is fixed and equals $P(k)=\delta_{kq}$, but the arrangement of links between nodes is random and changes in each update. Within such a reformulation, the algorithm of a  single update consists of 3 consecutive steps:

Randomly choose a spin $S_{i}$ 2. 2.

Calculate the change of the energy involved in the flip of $i$-th spin given by Eq. (1) 3. 3.

Flip the $i$-th spin with probability $W(\Delta E)$

In the previous paper, the case with $W(\Delta E)=W_{\text{M}}$ given by Eq. (3) was investigated Jed:Chm:Szn:15 . In this paper we will use several different dynamics including heat-bath given by Eq. (4). Beyond this generalization of $W(\Delta E)$, formally exactly the same algorithm was used in Jed:Chm:Szn:15 . However, it should be noted that previously we could not speak about the energy or the temperature because only $q$ out of all neighbors were considered.

III Results for the quenched $q$-RRG

We investigate the model by Monte Carlo simulations and estimate the magnetization of the system as an ensemble average over $M$ samples:

[TABLE]

where $m_{j}$ is the stationary magnetization in $j$-th sample, defined as:

[TABLE]

here $S_{i}^{j}$ denotes a value of an $i$-th spin in the $j$-th sample in the stationary state. In the case of the quenched approach for each $j$-th sample we build a new random regular graph of size $N$, in which degrees of all nodes are equal $q$, but the network does not change in time. In the case of the annealed approach we change the network after each single update defined by the algorithm in the previous section.

For the quenched approach, we expect that Metropolis and heat-bath algorithms will give exactly the same results. Moreover, results should be consistent with the general analytical formula for the critical temperature $T_{c}$ for the Ising model with coupling constant $J$ on an “equilibrium” random network with the degree distribution $P(k)$, derived by Dorogovtsev et al. Dor:Gol:Men:02 :

[TABLE]

For $q$-RRG degrees of all nodes are equal to $q$. Moreover, we have assumed that $J=1$; therefore, the critical temperature should be equal:

[TABLE]

Indeed, if we look at Fig. 1 we see that results for HB and M dynamics are the same, and the critical temperature is properly described by Eq. (8), which is also clearly visible in the left panel of Fig. 2.

IV Results for the annealed $q$-RRG

On one hand, theoretical description simplifies significantly if we replace the quenched network by the annealed one because in such a case we can easily derive the transition probabilities for the system and write down the rate and/or the master equation, as it was done in the previous paper Jed:Chm:Szn:15 . On the other hand, the behavior on the annealed graph, at least for the Ising model with Metropolis dynamics, is much more complex including: (1) a switch from a continuous to a discontinuous phase transition at $q=4$ and (2) an unexpected oscillatory behavior of the hysteresis, expanding for even values of $q$ and shrinking for odd values of $q$.

Here, following the reasoning presented in Jed:Chm:Szn:15 , we derive transition probabilities for the Ising model with HB dynamics on the annealed $q$-RRG. This allows us to calculate stationary magnetization as well as the critical temperature. In a single update, three events are possible – the number of spins ‘up’:

increases by 1 ($N_{\uparrow}\rightarrow N_{\uparrow}+1$) and simultaneously $m\rightarrow m+2/N$ with probability $\gamma^{+}$, 2. 2.

decreases by 1 ($N_{\uparrow}\rightarrow N_{\uparrow}-1$) and simultaneously $m\rightarrow m-2/N$ with probability $\gamma^{-}$, 3. 3.

or remains constant with probability $1-\gamma^{+}-\gamma^{-}$.

To simplify calculations we define a new macroscopic variable – the concentration of ‘up’ spins at time $t$:

[TABLE]

Using this variable we can write down the transition probabilities $\gamma^{+},\gamma^{-}$ for the infinite system:

[TABLE]

In above equations we have replaced $c(t)$ by $c$ for brevity of description, but one should remember that $c$ evolves in time, and this evolution is described by the rate equation Kra:Red:Ben:10 :

[TABLE]

where the unit time has been defined as $N$ elementary updates described by the algorithm in section II, as usually one Monte Carlo step is defined.

In the steady state $c(t+1)=c(t)=c$, which is equivalent to the condition:

[TABLE]

Now, we are ready to derive the formula for the critical temperature $T_{c}$. Let us define a quantity which can be thought of as a net force acting on our system in the following way Nyc:Szn:Cis:12 :

[TABLE]

In such a case, the lines of equilibrium are determined by the condition

[TABLE]

Differentiating the above equation with respect to the concentration gives

[TABLE]

Note that the critical point corresponds to the extremum of $T(c)$. Hence, the first derivative of the temperature over $c$ vanishes at $(c,T)=(1/2,T_{c})$, that is to say,

[TABLE]

Therefore, evaluating Eq. (15) at the critical point gives us the condition for its derivation

[TABLE]

Using formulas for the transition probabilities, we obtain the equation from which the critical temperature can be derived numerically

[TABLE]

For $q=4$, our solution coincides with the result presented in Lip:Gon:Lip:15 . Transition temperature as a function of $q$ is presented in the left panel of Figure 2. As seen, $T_{c}$ derived from Eq. (18) for an annealed $q$-RRG approaches the critical temperature derived by Dorogovtsev et al. Dor:Gol:Men:02 in the case of a quenched graph given by Eq. (8). Furthermore, it is seen in Fig.  2 that Metropolis algorithm, considered already in Jed:Chm:Szn:15 , gives different value of $T_{c}$ for the annealed graph.

V Different dynamics

Figure 1 illustrates differences of the out-comes from two popular dynamics: Metropolis and heat-bath. Furthermore, it has been shown that the heat-bath algorithm reproduces almost the same results for annealed and quenched graphs. The difference is only quantitative and decreases with the growing average degree $q$; see also Fig. 2 (left panel). On the other hand, Metropolis dynamics causes much more profound changes in the behavior of the system. On quenched networks only continuous phase transitions are observed whereas on annealed graphs both continuous and discontinuous transitions are present depending on the parameter $q$. Furthermore, the hysteresis exhibits oscillatory characteristics shrinking for odd values of $q$ and expanding for even ones Jed:Chm:Szn:15 .

In order to comprehend why annealing for some dynamics causes once only quantitative and other qualitative differences, we introduce and investigate several new dynamics which fulfill the detailed balance condition and reproduce the Boltzmann distribution in equilibrium for the model placed on quenched $q$-RRG. In other words, the transition probabilities satisfy Eq. (2).

The most natural choice for such a dynamics comes from the modification of the Metropolis algorithm:

[TABLE]

The above modification is identical with the original Metropolis dynamics, except of the case $\Delta E=0$. For M algorithm $W(\Delta E=0)=1$, on the other hand for the HB dynamics $W(\Delta E=0)=1/2$. Here we propose a  generalized version in which $W(\Delta E=0)=W_{0}$, where $W_{0}\in[0,1]$ is a parameter; see the right panel of Fig. 2. Such a modification reminds of the generalized zero-temperature Glauber dynamics introduced in God:Luc:05 , which occurred to be particularly interesting from theoretical point of view Lip:99 ; Rad:Vil:Mey:07 ; Sko:Szn:Top:12 ; Spi:Kra:Red:01 ; Ole:Kra:Red:11 . As we will see, the above dynamics, although artificial and not motivated by any physical processes, will also turn out to be particularly intriguing. Nevertheless, other transition rates that fulfill Eq. (2) are also possible, and we will briefly present results for two of them:

[TABLE]

where $C=\exp(-\beta q)$ is a normalizing constant and

[TABLE]

where $a$ is a parameter. Transition probabilities generated by the proposed dynamics are presented in Fig. 3 for several temperatures. Note that for growing $a$ the above transition rate approaches $W_{\text{M}^{\prime}}$ with $W_{0}=0$.

Although all dynamics satisfy detailed balance and therefore give the same results on quenched random $q$-regular graphs, the outcomes differ when we anneal the network. For $W_{\exp}$ phase transition is continuous for all values of parameters, like in the case of quenching, but diagrams are shifted towards higher temperatures. Consequently, it results in the larger values of the critical points.

More thought provoking behavior of the model is induced by the modified Metropolis algorithm $W_{\text{M}^{\prime}}$ with $W_{0}=0$. Then the system exhibits two types of phase transitions. For odd values of $q$ both Metropolis dynamics, modified and original one, are consistent and produce qualitatively the same results, in a sense that all transitions for $q>3$ are discontinuous. However, for even values of $q$ all transitions become continuous; see the left panel of Fig. 4. It should be recalled here that for $W_{0}=1$, which corresponds to original Metropolis results were almost opposite – though phase transitions were discontinuous for all $q>3$, yet just for even values of $q$ discontinuity was stronger (larger jump of the order parameter and hysteresis).

The right panel of Fig. 4 presents how the phase diagram changes under the influence of $W_{0}$ for the model with $q=4$. It is seen that there is a critical value of $W_{0}$ for which the transition type for even values of $q$ switches from discontinuous to continuous one. We can derive this critical point $W^{*}_{0}$ as well as categorize phase transitions based on the classical Landau theory Lan:37 ; Nyc:Szn:Cis:12 by introducing an effective potential $V(m)$ and studying its behavior for small values of the order parameter $m$, i.e., near zero. Having the net force $F(m,T)$, from Eqs. (13) and (9), we can construct the potential in the following way:

[TABLE]

In the vicinity of zero, we can expand it into a Taylor series, and neglect terms of the higher-order than $m^{6}$. This procedure leads to the potential in the form:

[TABLE]

where coefficients $A$, $B$, and $C$ depend on the model and can be derived directly from Eq. (22) as:

[TABLE]

When $A>0$, then the second derivative of Eq. (23) near $m=0$ is positive, and at this point $V(m)$ has a minimum; thus, a disordered phase is stable. On the other hand, when $A$ is negative, a disordered phase is unstable. Therefore, a condition $A=0$ allows us to derive the boundary of a region where a paramagnetic phase is present; note that this is equivalent to Eq.(17). The order of a transition is determined by the coefficient $B$. If it is positive, the system undergoes continuous phase transition . When it is negative, a discontinuous transition appears, and there is an area, bounded by the spinodal lines (dashed lines in Fig. 5), where two phases may coexist. Within the modified Metropolis dynamics, defined by Eq.(19), there is a ferromagnetic-paramagnetic continuous phase transition for small values of $W_{0}$ and above a certain tricritical point (TCP) the continuous transition line splits up into two spinodal lines, where the transition becomes first order (see Fig. 5). Similar behavior has been observed recently for the non-equilibrium Ising model with two heat-baths, but in that case the temperature of the heat-bath for links $T_{L}$ has driven the tricritical behavior Par:Hoh:17 .

Because TCP is a place where all three sections: the coexistence area and regions with purely ordered and disordered phases meet, it might be determined by fulfilling simultaneously two conditions $A=0$ and $B=0$. Moreover, note that TCP settles the value of $W_{0}^{*}$.

The entire phase diagram for the model with the modified Metropolis dynamics for q=4 is illustrated in Fig. 5. The analogical behavior is observed for grater even values of interacting neighbors.

In this particular case ($q=4$), the critical point lies at $W_{0}^{*}=11/15$, but in general its value depends on $q$. The dependence between the type of the phase transition and model’s parameters $q$ and $W_{0}$ for the modified Metropolis dynamics is presented in the top panel of Fig. 6. For discontinuous transitions, shading indicates the width of the hysteresis; darker areas correspond to wider metastable region $\Delta T$ where two phases may coexist. Oblique lines indicates the presence of continuous phase transitions. The boundary between these two transition types for a given $q$ corresponds to the tricritical point.

Based on the above result, we can draw a conclusion that by altering transition probabilities at $\Delta E=0$ one can control the nature of the phase transition. Moreover, lower flipping rates are connected with continuous transitions. In order to check whether indeed only the central point (i.e. $\Delta E=0$) of the dynamics has this profound impact on the model behavior we investigate the system with $W_{\tanh}$ rates, defined by Eq. (21). In this case, parameter $a$ does not affect $W(\Delta E=0)$, and for all conditions $W(\Delta E=0)=0$. Nonetheless, its higher value increases in general the flipping probabilities, which is clearly seen in Fig. 3. Note that large $a$ reproduces the outcome obtained for the model drove by the modified Metropolis algorithm with $W_{0}=0$ since then both dynamics coincide; look at Fig. 6 and compare the vertical section along $W_{0}=0$ in the top panel with the section along $a=20$ in the bottom one. However, there is a certain value of the parameter $a$ below which all transitions become continuous, like for $W_{\text{HB}}$ or $W_{\exp}$ dynamics.

VI Conclusions

Ising model, although already 100 years old still inspires many researchers and raises new questions Rad:Vil:Mey:07 ; Sko:Szn:Top:12 ; Spi:Kra:Red:01 ; Ole:Kra:Red:11 ; Roy:Bis:Sen:14 ; Roy:Sen:15 ; Men:Odo:00 . Recently we have introduced, under the name the $q$-neighbor Ising model, a seemingly small modification of the kinetic Ising model with Metropolis dynamics allowing each spin to interact only with $q$ neighbors Jed:Chm:Szn:15 . Surprisingly, this modification leads to a switch from a continuous to a discontinuous phase transition at $q=4$, and it causes an unexpected oscillating behavior of the hysteresis for $q\geq 4$ – it expands for even values of $q$ and shrinks for odd values of $q$. This was very puzzling phenomenon because extensive studies have shown not only that a phase transition in the equilibrium Ising model is continuous but also is characterized with the same set of bulk critical exponents on quenched and annealed networks Lee:09 .

However, recently it has been shown that the $q$-neighbor Ising model is a limiting case of a non-equilibrium system with two heat baths – one for the Ising spins at temperature $T_{s}$ and second for links of the graph at temperature $T_{L}$ Par:Hoh:17 . For $T_{L}=\infty$ such a generalized model reduces to the $q$-neighbor Ising model investigated in Jed:Chm:Szn:15 . This means that the $q$-neighbor Ising model is out of equilibrium on contrary to the Ising model on the annealed graph considered in Lee:09 . However, the question arises whether just being out of equilibrium causes tricritical behavior, i.e., switch from continuous to discontinuous phase transition. Maybe such an exotic behavior is observed exclusively for Metropolis algorithm? Furthermore, what is the role of dynamics for non-equilibrium systems?

Referring above questions, we have reformulated the $q$-neighbor Ising model in terms of the classical Ising model on the $q$-RRG and investigated it under several different dynamics belonging to the broad class of the generalized single-spin flip Glauber dynamics Glauber ; God:Luc:05 . On the quenched graph, which within the generalized model proposed in Par:Hoh:17 corresponds to the zero temperature of links ($T_{L}=0$), all dynamics fulfill the detailed balance condition and lead to the equilibrium Boltzmann distribution. Therefore, it is not surprising that they all give exactly the same results and reproduce analytical formulas for the critical temperature derived by Dorogovtsev et al. Dor:Gol:Men:02 .

On the annealed non-equilibrium network, which within the generalized model proposed in Par:Hoh:17 corresponds to the infinite temperature of links ($T_{L}=\infty$), each dynamics leads to completely different results. For example, the model with heat-bath algorithm on the annealed network displays always a continuous phase transition, and the critical temperature only slightly deviates from the equilibrium one. On the other hand, the generalized Metropolis algorithm gives qualitatively different behavior, which depends on the model’s parameter $W_{0}$. For $W_{0}=1$, we obtain the original Metropolis dynamics, which means that the transition switches to discontinuous at $q=4$, and hysteresis exhibits oscillatory behavior – expanding for even values of $q$ and shrinking for odd values of $q$. Yet, for $W_{0}<W_{0}^{*}(q)$, the behavior is even more exotic – phase transition switches alternately from continuous to discontinuous phase transition and back. Additionally, we have investigated two other types of dynamics, and it occurred that $W_{\exp}$ displays always a continuous phase transition while $W_{\tanh}$ switches between two types of the phase transitions. We have observed that generally higher flipping probabilities lead to discontinuous phase transitions.

We could of course propose many other dynamics which fulfill the detailed balance condition on the quenched graph, but it was not our aim. We asked the question about the role of dynamics for non-equilibrium models. On one hand, there is no excuse of referring to different dynamics or algorithms if the model is out of equilibrium, because in such a case there is no condition, such as detailed balance, which guarantees the unique steady state. On the other hand, referring to the particular dynamics within non-equilibrium models still happens Lip:Gon:Lip:15 ; Rad:Vil:Mey:07 . Such practices are not so much surprising, because it is quite common that methods of equilibrium statistical physics are transferred to the field of non-equilibrium systems and they give often reasonable results Hen:Hin:Lue:08 . Thus, one could have an intuition that the choice of a particular dynamics for non-equilibrium system will not influence significantly the outcome of the model.

Because results obtained in Jed:Chm:Szn:15 have shown that such an intuition may be very wrong, we have decided to examine several different dynamics and check how they would influence results in two limiting cases: quenched ($T_{L}=0$) and annealed ($T_{L}=\infty$). We have shown that indeed some dynamics (including heat-bath) are more resistant to the rewriring a network and display always continuous phase transition, whereas other (including Metropolis) are much more fragile and rewriring a network changes the type of the phase transition. Furthermore, our results suggest that we can generally expect discontinuous phase transitions for higher flipping rates $W_{0}$, i.e., in fact for larger noise, which is generally interesting result and coincides with observations made in Par:Hoh:17 , where discontinuity has been driven by higher temperature of a network. Of course there are more possibilities to to incorporate noise to the system. For example one could investigate the model on temporal network and introduce rewriring time $\tau$, as a parameter. Because $\tau=\infty$ corresponds to $T_{L}=0$ (i.e., quenched network) and $\tau=1/N$ to $T_{L}=\infty$ (i.e., annealed network), one could expect that there is a critical time scale $\tau^{*}$ that separates these two regimes. Indeed preliminary studies confirm these predictions Mar:Szn:17 .

It is known that switch from continuous to discontinuous phase transition can be caused by many factors. In the equilibrium statistical physics one of the best known examples is the regime-switch in the Potts model (for review see Wu:82 ): for $q>4$ the model undergoes a first-order phase transition, whereas a second-order phase transition for smaller values of $q$, with $q$ being the number of states of the spin. Moreover, it is common that systems exhibiting a discontinuous phase transition in high space dimensions may display a continuous transition below a certain critical dimension Hen:Hin:Lue:08 . For non-equilibrium systems it has been shown that by increasing the number of interacting neighbors, the fluctuations are diminished and the transitions become sharper or even change their type to discontinuous Nyc:Szn:Cis:12 ; Jed:Chm:Szn:15 ; Odo:Szo:96 . Yet, recent research has shown that for non-equilibrium systems switch from continuous to discontinuous phase transition and related tricriticality may be also caused by a high temperature of a network Par:Hoh:17 or larger flipping probabilities. These results suggest that incorporating certain type of noise to the non-equilibrium system can be responsible for tricriticality. Yet, this hypothesis requires further research.

However, the main outcome of this paper is the answer to the question posed in the introduction. Our aim was to check if being out of equilibrium, one can still speak about the Ising model with certain dynamics. As we have shown, each out of equilibrium dynamics leads to completely different behavior – some dynamics are more fragile for the non-equilibrium driving than others. In fact, it is even hard to speak about the dynamics or algorithm in case of such a non-equilibrium system because dynamics is an immanent feature of the model.

Acknowledgements.

This work was supported by funds from the National Science Centre (NCN, Poland) through grants no. 2015/18/E/ST2/00560 (to AC) and no. 2016/21/B/HS6/01256 (to AJ and KSW).