Ferromagnetic and spin-glass like transition in the $q$-neighbor Ising model on random graphs

The $q$-neighbor Ising model is investigated on homogeneous random graphs with a fraction of edges associated randomly with antiferromagnetic exchange integrals and the remaining edges with ferromagnetic ones. It is a nonequilibrium model for the opinion formation in which the agents, represented by two-state spins, change their opinions according to a Metropolis-like algorithm taking into account interactions with only a randomly chosen subset of their $q$ neighbors. Depending on the model parameters in Monte Carlo simulations phase diagrams are observed with first-order ferromagmetic transition, both first- and second-order ferromagnetic transitions and second-order ferromagnetic and spin-glass-like transitions as the temperature and fraction of antiferromagnetic exchange integrals are varied; in the latter case the obtained phase diagrams qualitatively resemble those for the dilute spin-glass model. Homogeneous mean-field and pair approximations are extented to take into account the effect of the antiferromagnetic exchange interactions on the ferromagnetic phase transition in the model. For a broad range of parameters critical temperatures for the first- or second-order ferromagnetic transition predicted by the homogeneous pair approximation show quantitative agreement with those obtained from Monte Carlo simulations; significant differences occur mainly in the vicinity of the tricritical point in which the critical lines for the second-order ferromagnetic and spin-glass-like transitions meet.


I. INTRODUCTION
The process of opinion formation has been a rapidly growing subject of research in statistical physics in the last decades [1]. It is often investigated by means of nonequilibrium models formed of agents expressing discrete opinions on a given subject, placed in nodes and interacting via edges of a fixed network with topology reflecting the complexity and, possibly, heterogeneity of social contacts [2,3]. Widely studied examples of such models comprise, e.g., the voter model [4][5][6][7][8], the q-voter model called also nonlinear voter model [9][10][11][12], variants of the noisy q-voter model with different forms of stochasticity [13][14][15][16][17][18][19][20][21][22][23], in particular the q-voter model with independence or anticonformism [13][14][15][16]22], the q-neighbor Ising model [24][25][26][27][28] and the majority-vote model [29][30][31][32][33][34][35]. As a rule, in these models agents can express only one of two possible opinions and thus are represented by two-state spins and interactions between agents have a form of exchange interactions, either explicitly, as in the q-Ising model, or effectively, in the remaining models. Particular attention is devoted to models with ferromagnetic (FM)-like (henceforth termed FM) "friendly" interactions which tend to align the orientations of spins (agents' opinions) in parallel, and hence to transition to consensus with decreasing degree of stochasticity, similar to FM transition in magnetic systems. This transition in the noisy q-voter, q-neighbor Ising and majorityvote models can be continuous or discontinuous depending on the parameters of the model, degree distribution of the network and details of the spin-flip rate. In the case of models on complete graphs the FM-like (henceforth termed FM) transition can be quantitatively described in the framework of mean-field approximation (MFA) [13,14,16,17,20,22,24,26]. Concerning the models on networks with different degree of heterogene-ity the MFA yields quantitatively correct predictions in the case of the majority vote model [31,32,35], while in the case of the q-voter and q-neighbor Ising models better quantitative agreement is obtained between numerical simulations and predictions of a more sophisticated homogeneous or heterogeneous pair approximation (PA) [7,8,15,18,19,21,23,28,33,[35][36][37].
Apart from FM interactions also antiferromagnetic (AFM)-like (henceforth termed AFM) ones, which prefer antiparallel orientation (opposite opinions) of interacting spins (agents) can be introduced in the above-mentioned models. For example, anticonformist agents can be included, preferring to express opinion opposite to that suggested by their neighborhood, which corresponds to asymmetric AFM interactions. Increasing the fraction of anticonformists in the q-voter model [13,22] and the majority vote model [34] leads to disappearance of the FM transition so that the models remain in the disordered paramagnetic (PM) phase even for small degree of stochasticity. Besides, in the majority vote model with large fraction of anticonformists a spin glass (SG)-like (henceforth termed SG) phase can appear [34] characterized by the absence of long-range ordering (hence, absence of global consensus) and presence of only shortrange ordering [38][39][40][41]. Symmetric AFM interactions can also be included, corresponding to "hostile" interactions between agents who thus tend to express opposite opinions. Then, for the majority vote model on random graphs phase diagram is obtained resembling qualitatively that for the equilibrium Ising model on random graphs which is a model for a dilute SG [42], with secondorder transitions from the PM to the FM and SG phases and a tricritical point (TCP) in which critical lines corresponding to the two above-mentioned transitions meet [35].
The aim of this paper is to investigate the effect of sym-metric AFM interactions between agents on phase transitions in a nonequilibrium model for the opinion formation in which both first-and second-order FM transitions are possible. For this purpose the q-neighbor Ising model on random graphs is studied which belongs to a family of nonequilibrium Ising models in which the spins in nodes and the edges of the network of interactions are in contact with thermal baths with different temperatures [25]. The q-neighbor Ising model is more complicated than noisy q-voter models but has an advantage that its phase diagram in the presence of AFM interactions can be roughly predicted following physical intuitions from the equilibrium Ising model for dilute SG [42]. In this paper the model is investigated only on homogeneous or weakly heterogeneous networks in order to focus on the effect of the AFM interactions and avoid complications related to heterogeneity of the network of interactions (e.g., the need to use high-dimensional system of differential equations resulting from heterogeneous PA [36,37] for theoretical description of the model). The paper is organized as follows. In Sec. II the model is introduced and the main results concerning the observed phase transitions are outlined. In Sec. III the MFA and PA for the q-neighbor Ising model [27,28] are extended to take into account the presence of the AFM symmetric interactions between agents; since it is assumed that the network of interactions is not strongly heterogeneous only homogeneous MFA and PA are considered. In Sec. IV the FM transition in the model is investigated in the framework of the above-mentioned approximations and conditions are discussed for the occurrence of the first-and secondorder transition and of the TCP in which the critical lines for these transitions meet. Theoretical predictions for the FM transition, in particular those from the PA, are shown to agree quantitatively with results of Monte Carlo (MC) simulations. Besides, also in Sec. IV numerical evidence for the second-order SG transition in the model is provided and conditions are discussed for its occurrence and for the presence of the TCP in which the critical lines for the second-order FM and SG transitions meet. Finally, Sec. V is devoted to summary and conclusions.

II. THE MODEL
A starting point for the definition of the q-neighbor Ising model on networks is the usual Ising model with two-state spins s i = ±1 located in nodes i = 1, 2, . . . N and with exchange integrals J ij associated with edges of a network of interactions, with the Hamiltonian In this paper the networks of interaction under study are limited to a class of homogeneous and weakly heterogeneous random graphs, such as random regular graphs (RRGs) with a degree distribution P (k) = δ k,K and mean degree of nodes k = K or Erdös-Rényi graphs (ERGs) with binomial degree distribution P (k) = N −1 k p k (1 − p) N −1−k and mean degree of nodes k = (N −1)p, p 1 [43]. As an extension with respect to previous studies of the q-neighbor Ising model [24][25][26][27][28] it is assumed that the exchange integrals J jj are drawn form a probability distribution (J = 1 is assumed in simulations) and randomly associated with the edges; if the nodes j, j are not connected by an edge J j,j = 0 by definition. In this way interactions between spins are either FM (J jj = J > 0, with probability 1 − r) or AFM (J jj = −J < 0, with probability r) and the Hamiltonian (1) is that for a model for dilute SG [42]. Each simulation of the q-neighbor Ising model described below starts with random realization of the above-mentioned network of interactions with a given degree distribution P (k) and random association of the exchange integrals with edges according to the distribution (2); both the network and the exchange integrals remain fixed in the course of simulation (quenched disorder).
The q-neighbor Ising model is a nonequilibrium variant of the above-mentioned Ising model in which, at each time step, each spin interacts only with q randomly chosen neighbors. The dynamics of the q-neighbor Ising model on networks is a modification of that of the kinetic Ising model with Metropolis spin-flip rate. MC simulations of the model are performed using random sequential updating of spins, with each MC simulation step (MCSS) corresponding to updating all N spins; the possible spin flips correspond to changes of opinions of the agents. Each MCSS is performed as follows: 1. Randomly choose a node j, 2. From the set of k j neighbors of the node j choose randomly and without repetitions a subset nn j,q of its q neighbors (q-neighborhood), 3. Calculate first the local energy-like quantity E j s j , s j :j ∈nnj,q by summing only these terms in the Hamiltonian (1) which account for the interactions of the spin s j with spins belonging to its q-neighborhood, then the same quantity corresponding to the model with the spin s j flipped and the remaining spins unchanged E j −s j , s j :j ∈nnj,q , and the change of the local energy-like quantity ∆E j related to the potential flip of the spin s j , 4. Flip the spin s j with probability given by a Metropolis-like formula where β = 1/T and T is the effective temperature which measures the level of internal noise (uncertainty in agents' decision making), 5. Repeat steps (1)(2)(3)(4) until all N nodes are updated (without repetitions).
The q-neighbor Ising model, defined as above, with uniform FM exchange integrals (all J ij = J > 0) was investigated in detail on a fully connected graph [24][25][26][27] and networks with finite mean degree k [28]. In the former case the model exhibits first-order FM transition for q = 4 and q ≥ 6 with clearly visible hysteresis loop. Width of the hysteresis loop in general increases with q, though there are oscillations superimposed on this trend such that loops for consecutive odd values of q are narrower than for the neighboring even values of q [24]. The same is true for the model on networks provided that q k . However, as q is increased and becomes comparable with k the hysteresis loop becomes narrower and eventually disappears, and the FM transition becomes second-order [28].
In this paper the effect of AFM interactions (with J ij = −J < 0) on phase transitions in the q-neighbor Ising model is studied. Results of both numerical and theoretical investigations can be briefly summarized as follows. In general, as the fraction r of the AFM interactions increases the critical temperature (or temperatures in the case of the first-order transition) for the FM transition decreases. For q k and small r the FM transition is first-order for q = 4 and q ≥ 6, but for q ≥ 6 the width of the hysteresis loop decreases to zero with r. Eventually a TCP occurs, the FM transition becomes second-order and then disappears as the critical temperature reaches zero, and for large r the PM phase is stable for any temperature. For q comparable with k the FM transition is continuous for any r; besides, for large r continuous SG transition occurs and the phase borders between the PM and FM phases and between the PM and SG phases meet in a TCP. Thus, in the latter case the obtained phase diagram for the q-neighbor Ising model qualitatively resembles that for dilute SG [42].

III. THEORY
In this section two theoretical approaches to the FM transition in the q-neighbor Ising model with FM and AFM interactions are presented, the (homogeneous) MFA which is exact for the model on fully connected graphs and provides good approximation for the model on random graphs with large mean degree k as well as homogeneous PA which provides better approximation for a model on random graphs with arbitrary k . These approaches are extensions of the MFA and PA for the q-neighbor Ising model with purely FM exchange integrals [24,28]. In particular, predictions of the homogeneous PA show good quantitative agreement with results of MC simulations for a broad range of the model parameters q, k , r, as discussed in Sec. IV. Unfortunately, the SG transition which is also observed in MC simulations of the q-neighbor Ising model cannot be analyzed within these approaches: the order parameter for the SG phase is based on two-spin correlation function [39][40][41], while such correlations are neglected in the MFA and PA.

A. Mean field approximation
Before introducing homogeneous PA for the q-neighbor Ising model on random graphs, in this subsection simple MFA for the q-neighbor Ising model on a fully connected graph with FM and AFM exchange integrals drawn from the distribution (2) is presented. It is a straightforward extension of the MFA for the q-neighbor Ising model on a fully connected graph with uniform FM exchange inegrals [24], corresponding to r = 0 in the distribution (2). In this model the neighborhood of a given spin consists of all remaining N − 1 spins. In the MFA the macroscopic quantity characterizing the model is the concentration c of spins with direction up, related to the order parameter, the usual magnetization m, by c = (1 + m)/2. A dynamical equation for the concentration c has a form of the rate equation, where γ + (γ − ) are rates of spin flips in the direction up (down) averaged over all spins. In the thermodynamic limit N → ∞ these rates can be evaluated as follows. Let us define inconsistency of opinions of a pair of interacting agents as orientations of the corresponding spins such that their interaction increases energy of the associated equilibrium Ising model with the Hamiltonian (1). Then, spins in nodes j, j have inconsistent orientations if they have opposite orientations and interact via FM exchange integral J jj = J > 0 or they have the same orientations and interact via AFM exchange integral J jj = −J < 0. Hence, probability that a spin s j = −1 has a neighbor with inconsistent orientation is p = (1 − r)c + r(1 − c) = (1 − 2r)c + r, and the number i of neighbors with inconsitent orientations among its N − 1 neighbors obeys a binomial distribution If in a single time step such spin attempts to flip and among its q selected neighbors there are l ones with inconsitent orientation (l ≤ i) change of the local energy (4) caused by flipping this spin is ∆E j = −2Jl + 2J(q − l) = 2J(q − 2l). Then from Eq. (5) the average spin flip rate for spins possessing i neighbors with inconsistent orientation is where Since all nodes are statistically equivalent the spin flip rate averaged over all spins γ + (c, T ) in Eq. (6) can be obtained by averaging the rate f (i, T ), Eq. (7), over the binomial distribution of the number of neighbors of a representative spin s j = −1 with inconsistent orientation and multiplying the result by the concentrations of spins with orientation down. The rate γ − (c, T ) in Eq. (6) can be obtained similarly by repeating the above reasoning for a representative spin s j = +1. The results are Rewriting Eq. (6) in terms of the magnetization m = 2c − 1 it is obtained that Expanding the right-hand side of Eq. (11) in powers of m it is possible to write it as a derivative of an effective potential V (m, T, r, q), V (m, T, r, q) = C 2 (T, r, q)m 2 + C 4 (T, r, q)m 4 + C 6 (T, r, q)m 6 + . . . etc.
It can be seen that Eq. (11) and (12) have a fixed point m = 0 corresponding to the PM phase. In general, this fixed point is stable for high temperatures T and high fractions of AFM exchange interactions r and for fixed q can lose stability as T or r are decreased with the other parameter kept constant which corresponds to the transition from the PM to the FM phase. Let us assume that r is fixed and focus on the FM transi-tion occurring with decreasing temperature at a critical point T c which can be determined numerically from the condition C 2 (T c , r, q) = 0. This condition can be fulfilled only if r < r M F A (q) for which T c = 0; otherwise, the PM phase is stable in the whole range of T > 0. Character of the FM transition can be deduced using the phenomenological Landau formalism for phase transitions: the transition is second-order if C 4 (T c , r, q) > 0 and first-order if C 4 (T c , r, q) < 0 and simultaneously C 6 (T c , r, q) > 0 (if the latter condition is not fulfilled, higher-order coefficients in the expansion (13) should be positive). The critical temperature for the second-order FM transition is T c,M F A a pair of symmetric stable fixed points of Eq. (11) with |m| > 0 appears corresponding to the FM phase. In the case of the first-order FM transition the temperature T c1,M F A = T c is the lower critical temperature below which the PM phase loses stability and the only stable fixed points of Eq. (11) are the two symmetric fixed points with |m| > 0 corresponding to the FM phase. The latter fixed points are in turn stable for T < T c1,M F A ) which is the higher critical temperature above which the only stable fixed point is that with m = 0 corresponding to the PM phase. T The number of equations may be reduced since it is possible to evaluate r form the condition C 2 (T, r, q) = 0 using Eq. (14), , substitute in the right-hand side of Eq. (15) and solve numerically the resulting equation For the sake of completeness should be mentioned that the above version of the MFA is usually called homogeneous MFA. A more advanced heterogeneous MFA usually provides better description of models on networks with different degree of heterogeneity. In the latter approximation dynamical variables may be concentrations of nodes with degree k in which spins have orientation up. However, for the q-neighbor Ising model (as well as for the q-voter model) it can be shown that in the heterogeneous MFA these concentrations do not depend on the degree of nodes, thus heterogeneous MFA is equivalent to the homogeneous PA presented above.

B. Pair approximation
In the framework of the homogeneous PA the qneighbor Ising model on networks is described in terms of concentrations c k,↑ of nodes with degree k in which spins have orientation up (normalized to the number of nodes with degree k which is N P (k); the respective concentration of nodes with spins with orientation down is c k,↓ = 1 − c k,↑ ) and concentration b of active links (normalized to the total number of edges N k /2). In this section dynamical equations for these concentrations are derived, which represent an extension of Eq. (6) and (11) obtained in the MFA.
Let us start with the definition of active links. In the case of model with only FM interactions active links are associated with edges connecting nodes containing spins with opposite orientations, i.e., they correspond to interactions between pairs of spins increasing the energy (1) in the related equilibrium Ising model. Hence, in the model under study with FM and AFM interactions a natural generalization of the concept of active links is to assume that they correspond to interactions between pairs of spins which increase the energy (1), i.e., that active links connect pairs of neighbors with inconsistent opinions defined in Sec. III A. Thus, since J j,j = ±1 if an edge connecting nodes j, j exists and J j, Assumption of homogeneity in the PA consists precisely in the assumption that the dynamics of active links can be described by a single concentration b rather than a set of concentrations of active links connecting pairs of nodes with given degrees, as in the case of more advanced heterogeneous pair approximation [8,36,37]. This assumption is valid for models on homogeneous and weakly heterogeneous networks such as RRGs and ERGs investigated in this paper. Nevertheless, in the homogeneous PA the possible heterogeneity of the network is partly reflected since it is allowed that the concentrations c k,↑ depend on the degree of nodes. Hence, homogeneous PA differs from the (possibly heterogeneous) MFA in that apart from the degree-dependent concentrations of nodes with spins directed up also the concentration b of active links is treated as a separate dynamical variable. Another basic assumption in the PA is that orientations of spins in the neighboring nodes are not mutually correlated, thus the number of active links i (i ≤ k j ) attached to the node j with degree k j containing spin with orientation ν, (ν ∈ {↑, ↓}) obeys a binomial distri- Here, θ ν are conditional probabilities that a link is active provided that it is attached to a randomly chosen node containing spin with orientation ν. Due to the above-mentioned homogeneous approximation these probabilities can be expressed in terms of the macroscopic concentrations c k , b. In the presence of FM and AFM interactions this can be done in the following way [35]. Since the number of nodes in the graph is N , the total number of links is N k /2, the number of FM links with J ij = J > 0 is (1 − r)N k /2, the number of AFM links with J ij = −J < 0 is rN k /2, the number of active links is N k b/2, each link has two tips (ends) attached to two different nodes, thus the total number of tips is N k , the number of tips of FM links attached to the nodes is (1−r)N k , the number of tips of AFM links attached to the nodes is rN k and the number of tips of active links attached to the nodes is N k b. Let us denote by P (ν, ν ), where ν, ν ∈ {↑, ↓}, concentration of tips of links attached to the nodes with spins with orientation ν such that the other tip of the link is attached to a node with spin with orientation ν , normalized to the total number of tips. Hence, the corresponding number of above-mentioned tips is N k P (ν, ν ). In order to proceed with calculation it should be assumed that signs of the exchange integrals associated with subsequent links are not correlated with orientations of spins in the nodes to which these links are attached. Then, using the definition of an active link it is obtained that the number of tips of active links can be expressed as Obviously, P (↓, ↑) = P (↑, ↓) and ν,ν ∈{↑,↓} P (ν, ν ) = 1, thus Besides, the number of tips of links attached to nodes with spin with orientation ν is N k ν ∈{↑,↓} P (ν, ν ); this number can be also expressed in terms of the concentrations c k,ν as k N P (k)kc k,ν = N k C ν , where C ν = k −1 k P (k)kc k,ν is a "weighted", or "link", concentration of nodes with spins with orientation ν, such that C ↓ = 1 − C ↑ . Hence, for ν ∈ {↑, ↓}. The conditional probability θ ν can be evaluated as the ratio of the number of tips of active links attached to nodes with spins with orientation ν to the total number of tips of links attached to such nodes N k C ν . Using Eq. (18,19) the conditional probabilities can be eventually expressed as Taking into account the aforementioned assumptions and results, dynamical equations for the concentrations c k , b in the PA can be eventually obtained. The equations for the concentrations c k,↑ of nodes with spins with orientation up have again a form of rate equations, where γ + (γ − ) are rates of spin flips in the direction up (down) averaged over all spins located in nodes with degree k. These rates can be evaluated in the same way as the corresponding rates in the MFA by replacing the number of neighbors N − 1 by the degree k, the number of neighbors with inconsistent orientation i by the number of attached active links and the probability p that a node has a neighbor with inconsistent opinion with appropriate probabilities θ ν , ν ∈ {↓, ↑}, Eq. (20). Thus the spin-flip rate for spins in nodes with i active links attached, provided that they have degree k, is and the average rates in Eq. (21) are Substituting Eq. (23) in (21) can be easily seen that for any combination of degrees c k,↑ − c k ,↑ → 0 with increasing time, thus the system of equations (21) has a stable stationary solution c k,↑ = c ↑ for any degree k. Hence there is also C ↑ = c ↑ , thus c ↑ = c, where c is the usual (unweighted) concentration of nodes with spins up, occurring also in the MFA in Sec. III A and related to the usual magnetization via c = (1 + m)/2. It follows that in order to characterize the stationary states (corresponding to thermodynamic phases) of the q-neighbor Ising model in the framework of the homogeneous PA it is sufficient to replace all concentrations c k,↑ with a single concentration c ↑ = c and c k,↓ with c ↓ = 1 − c. This situation is identical as in the case of the q-voter model [18]. The dynamical equation for the concentration of active links b can be obtained by observing that each flip of a spin located in the node with degree k and i active links attached (i ≤ k) changes b by ∆ b (i|k) = 2 N k (k−2i), and such changes happen at the rate f (i, T |k) (22). Averaging over the binomial distribution of the number of active links attached to nodes with degree k containing spin with a given orientation, over the concentrations c k,ν , ν ∈ {↓, ↑} of nodes with degree k containing spin with orientation ν and over the two possible orientations of spins, approximating again c k,↑ = c ↑ = c, c k,↓ = c ↓ = 1 − c, taking into account that that the elementary time step is ∆t = 1/N and going to the thermodynamic limit N → ∞ yields a general dynamical equation for the concentration of active links [15], (24) The last two sums in Eq. (24) can be evaluated as in Ref. [28]. Eventually, the following system of dynamical equations for the concentrations c of nodes with spin with orientation up and b of active links is obtained, where For r = 0, when there are only FM interactions associated with edges, θ ↓ = b/2(1−c), θ ↑ = b/2c and equations for the q-neighbor Ising model from Ref. [28] are recovered.
thus the eigenvalues of the Jacobian matrix are Let us again assume that r is fixed and focus on the possible FM transition occurring with decreasing temperature. For the parameters used in the MC simulations below numerical analysis of Eq. (28), (30), (31) reveals that for θ being a solution of Eq. (28) the eigenvalue λ 2 < 0 in the whole range of T , while λ 1 can change sign with varying T provided that r < r P A ( k , q), where r P A ( k , q) is also determined numerically. Thus, for fixed r, the critical temperature T c at which the PM solution with c = 1/2 loses stability as well as the corresponding value θ c are determined from simultaneous solution of equations B(c = 1/2, θ) = 0, Eq. (28), and λ 1 = 0, Eq. (30). Depending on the order, the FM transition occurring at, or in the vicinity of, T = T c , corresponds to different bifurcations of the fixed point or points (including the PM fixed point) of the two-dimensional system of equations (25)(26). For k → ∞ predictions concerning the occurrence and order of the transition from the PM to the FM phase obtained using the PA and MFA coincide for any q. For the q-neighbor Ising model on networks with finite mean degree k , provided that q k predictions of the PA and MFA are still qualitatively similar, with quantitative differences becoming more pronounced with decreasing k or increasing q. In particular, for r = r P A ( k , q) the above-mentioned T c = 0 and for r < r P A ( k , q) the FM transition occurs with decreasing T which can be secondor first-order, depending on the parameters k , q, r. In the case of the second-order transition the PM fixed point loses stability via a supercritical pitchfork bifurcation at T (F M ) c,P A = T c and for T < T (F M ) c,P A a pair of stable equilibria with c > 1/2 (m > 0), b < 1/2, or c < 1/2 (m < 0), b < 1/2, exists, corresponding to the FM phase with positive or negative magnetization, respectively. In the case of the first-order transition as T is decreased two pairs of stable and unstable equilibria appear via two saddle-node bifurcations taking place simultaneously at temperature c2,P A > T c , which can be determined only numerically. For T c < T < T (F M ) c2,P A the two above-mentioned stable equilibria, one with c > 1/2 (m > 0), b < 1/2, and the other with c < 1/2 (m < 0), b < 1/2, corresponding again to the FM phase with positive or negative magnetization, respectively, coexist with the stable equilibrium with c = 1/2 (m = 0), b ≤ 1/2 corresponding to the PM phase; the basins of attraction of the three stable equilibria are separated by stable manifolds of the two unstable equilibria. Eventually at T = T If q and k are comparable predictions concerning the FM transition based on the PA are qualitatively different from those based on the MFA. For q > k /2 from the PA follows that the FM transition with decreasing T is always second-order and occurs for 0 < r < r P A ( k , q) at the critical temperature T c = T (F M ) c,P A . This prediction is reasonable since, in particular, the PA predicts that the FM transition in the q-neighbor Ising model on a RRG with P (k) = δ k,K is continuous in the limiting case q = K corresponding to the equilibrium Ising model on a RRG, as expected. Moreover, at r = r P A ( k , q) there is T c > 0; for fixed k the value of r P A ( k , q) weakly depends on q, but the critical temperature T c at r = r P A ( k , q) increases with q. Besides, for a range of r below r P A ( k , q) as temperature is further reduced the two symmetric stable fixed points with m > 0 or m < 0 and b < 1/2 which exist for T < T fork bifurcation. This means that for a range of r below r P A ( k , q) the PA predicts occurrence of another critical line T (F M ) c,P A (T, r) corresponding to a continuous transition from the FM to the PM phase with decreasing temperature. This line merges with that for the usual transition from the PM to the FM phase T (F M ) c,P A (T, r) at r = r P A ( k , q), and the two critical lines form a characteristic cusp marking the borders of stability of the FM phase (Fig. 1). The range of r for which the additional critical line occurs increases with q and for q = k it comprises a whole interval 0 ≤ r ≤ r P A ( k , q). Anticipating results of MC simulations (Sec. IV) it can be shown that the critical temperatures T c,P A at r = r P A ( k , q) are very close, or even equal, to the critical temperature for the SG transition T (SG) c,M C for given q obtained from simulations (Fig. 1). Thus, the abovementioned cusp at the borderline of the FM phase in the framework of the PA might be interpreted as a TCP in which the critical lines for the FM and SG transitions from the PM phase meet. It is then tempting to speculate that the predicted critical line T (F M ) c,P A (T, r) is related to the de Almeida -Thouless line determining the lower border of stability of the FM phase in the replicasymmetric solution [39][40][41]45], which exists also in the model for dilute SG [42]. However, the de Almeida -Thouless instability of the FM phase with decreasing temperature leads to the occurrence of the re-entrant SG phase characterized by non-zero magnetization rather than the PM phase with zero magnetization predicted by the PA in the q-neighbor Ising model. It is interesting to note that a similar additional critical line marking the lower border of stability of the FM phase for decreasing temperature-like model parameter is predicted by the PA in the majority-vote model with FM and AFM interactions [35], even with a more complex structure (firstorder transition from the FM to the PA phase is possible). Hence, appearance of such line can be typical for the PA in models with FM and AFM interactions. However, it should be emphasised that results of MC simulations of the q-neighbor Ising model differ significantly from the predictions of the PA for r in the vicinity of r P A ( k , q) for any q, and the numerically obtained location of the TCP in which the critical lines for the FM and SG transitions from the PM phase meet is far from the above-mentioned cusp at the border of stability of the FM phase resulting from the PA. Similar discrepancy is also observed in the majority-vote model [35].

IV. RESULTS AND DISCUSSION
In order to verify the occurrence of the FM or SG phase transition MC simulations of the q-neighbor Ising model under study were performed and in the case of the FM transition their results were compared with predictions of the MFA and homogeneous PA from Sec. III. In this section results of simulations of the model on RRGs are only presented; in most cases results for the model on ERGs with the same parameters k , q, r are quantitatively similar. Simulations were performed on networks with the number of nodes 10 3 ≤ N ≤ 10 4 using simulated annealing algorithm with random sequential updating of the agents' opinions. For each realization of the network and of the distribution of the exchange integrals simulation is started in the disordered PM phase at high temperature with random initial orientations of spins. Then the temperature is decreased in small steps toward zero, and at each intermediate value of T , after a sufficiently long transient, the order parameters for the FM and the possible SG transitions are evaluated as averages over the time series of the opinion configurations. Alternatively, to check for the presence of the hysteresis loop in the first-order FM transition, the above algorithm can be applied with FM initial conditions and temperature increased. The results are then averaged over 100 − 500 (depending on N ) realizations of the network and of the distribution the exchange integrals.
The possible FM and SG transitions in the q-neighbor Ising model are investigated in the same way as in the corresponding equilibrium Ising model. The order parameter for the FM transition is the absolute value of the magnetization wherem denotes a momentary value of the magnetization at a given MCSS, · t denotes the time average for a model with given realization of the network according to P (k) and of the associated distribution of the exchange integrals according to P (J ij ), and [·] av denotes average over different realizations of the network and of the distribution of exchange integrals. The order parameter for the SG transition (henceforth called the SG order parameter) is the absolute value of the overlap parameter [39][40][41], where α, β denote two copies (replicas) of the system simulated independently with different random initial orientations of spins andq is a momentary value of the overlap of their spin configurations at a given MCSS. In the PM phase both M and Q are close to zero. In the case of the FM transition both M and Q increase as T is decreased.
In the case of the SG transition the SG order parameter Q increases as T is decreased while the magnetization M remains close to zero. The order of the FM or SG transition and the critical values of the temperature can be conveniently determined using respective Binder cumulants U (M ) vs. T and U (Q) vs. T [44], In the case of the second-order FM or SG transition the respective cumulants are monotonically decreasing functions of temperature: for T → 0 there is U (M ) → 1 in the FM phase and U (Q) → 1 in the SG phase, and for c,M C for the FM and SG transitions can be determined from the intersection point of the respective Binder cumulants for models with different numbers of agents N [44]. In the case of the first-order FM transition it is sometimes possible to observe directly the hysteresis loop, by measuring magnetization M as a function of decreasing temperature for a model started in the PM phase as well as as a function of increasing temperature for a model started in the FM phase with all spins directed up (or down). This requires a model with large enough number of nodes N and large enough width of the hysteresis loop. If the hysteresis loop is narrow the Binder cumulants again become useful. In the case of the first-order transition behavior of the cumulant U (M ) for T → 0 and T → ∞ is similar as in the case of the second-order transition, but close to the critical temperature the cumulant exhibits negative minimum which deepens and becomes sharper with increasing number of nodes N . The critical temperature for the first-order FM transition again can be determined from the intersection point of the cumulants U (M ) for models with different numbers of agents N , started, e.g., in the PM phase.
Since evidence for the second-order SG transition in the nonequilibrium q-neighbor Ising model is purely numerical, care was taken a discarded transient after each change of temperature in the simulated annealing algorithm was long enough, as well as averaging over time and over different realizations of networks and of the distribution of exchange integrals in MC simulations were performed over long enough time intervals and large enough number of realizations, respectively, to obtain reliable values of Q, U (Q) as functions of T for fixed r. This was further verified by evaluating the critical temperature for the SG transition using the Binder cumulants for the model with q = 50 and r = 1 on a RRG with k = K = 50. This model is equivalent to the equilibrium Ising model on a RRG with purely AFM exchange integrals. The obtained critical temperature T q. It is known that for r = 0 (purely FM interactions) this model exhibits first-order FM transition characterized by a particularly wide hysteresis loop [24,28]. Both the MFA and homogeneous PA predict that the FM transition occurs for a small range of r > 0 and it remains first-order, still with a broad hysteresis loop (Fig. 2); the upper and lower critical temperatures determining borders of the hysteresis loop predicted by the MFA are slightly lower than the corresponding critical temperatures predicted by the PA. Moreover, since the lower critical temperature T c1,P A ) at which the PM phase loses stability with decreasing temperature reaches zero at smaller value of r = r M F A (r = r P A ) than the upper critical temperature T c2,P A ), both theories predict that for a certain range of r there is no transition from the PM to the FM phase with decreasing temperature and the PM phase remains stable for T → 0, and simultaneously the FM phase with increasing temperature remains stable up to T = T c2,M C (Fig. 2). It seems that the width of the hysteresis loop obtained from simulations is signif- c2,M C reach zero could not be accurately estimated due to rapid decrease of these critical temperatures at the end of the critical lines, in accordance with predictions of the MFA and PA. Finally, it should be mentioned that for q = 4, K = 50 SG transition characterized by increase of the SG order parameter Q with decreasing temperature was not observed, and for large r the PM phase remains stable as T → 0.
More complex phase diagrams are predicted theoretically and observed in MC simulations for the q-neighbor Ising model on random graphs with q ≥ 6 and q k . As an example, the phase diagram for the model with q = 8 on RRG with k = K = 50 is shown in  c,M C and exhibit negative minima as functions of T which become deeper with increasing number of nodes, and the magnetization M changes discontinuously at the critical point ( Fig. 4(a)). For r ≥ 0.15 the FM transition becomes second-order: up to N = 10 4 the Binder cumulants U (M ) for different N cross at one point corresponding to the critical temperature T (F M ) c,M C and monotonically decrease with T , and the magnetization M changes smoothly at the critical point ( Fig.  4(b)). The critical temperatures as well as the extent of the bistability region associated with the first-order FM transition and the location of the TCP separating it from the second-order FM transition are quantitatively well predicted by the homogeneous PA, while predictions of the MFA are noticeably worse. Hovever, as r approaches r P A the critical temperature T c,P A before decreasing sharply, and the range of r for which the FM transition occurs turns out to be slightly wider than predicted by the PA. Again, SG transition is not observed and for large r the PM phase remains stable as T → 0.
Another kind of complex phase diagrams is predicted theoretically and observed in MC simulations for the qneighbor Ising model on random graphs with q > k /2. As an example, the phase diagram for the model with q = 46 on a RRG with k = K = 50 is shown in  (Fig. 5). It is remarkable that in the case of the q-neighbor Ising model this agreement is much better than, e.g., in the case of the q-voter model on random graphs, where predictions of the PA deviate much from results of MC simulations as q approaches k [24]. Only in the vicinity of r P A = 0.36 . . . for which T c,M C remained stable as T → 0. In MC simulations of the model with q > k /2, apart from continuous FM transition, for larger values of r also second-order SG transition is observed. It is characterized by increase of the SG order parameter Q with decreasing temperature, and the critical temperature for this transition T c,M C = 5.75 ± 0.05. It is noteworthy that the SG transition occurs in the nonequilibrium q-neighbor Ising model on random graphs with such combinations of parameters q, k that the resulting phase diagram in Fig. 5 qualitatively resembles that for the model for dilute SG [42]: it contains only critical lines for the second-order FM and SG transitions, and the first-order FM transition is not observed. Similar phase diagrams were observed also in another nonequilibrium counterpart of the Ising model on random graphs, the majority-vote model [35].
The fact that the SG transition in the nonequilibrium q-neighbor Ising model on random graphs occurs only if q > k /2 may be, perhaps naively, understood as follows. It is known that the energetic landscape of the corresponding equilibrium Ising model with the Hamiltonian (1) and mixture of FM and AFM exchange integrals is littered with many local minima corresponding to metastable spin configurations [38][39][40][41][42]. In the nonequilibrium model there is no energetic landscape, however, it may be speculated that in the case of homogeneous graphs if the condition q > k /2 is fulfilled, then in each consecutive MCSS the network of interactions (varying in time due to random and in general not symmetric choices of the q-neighborhoods for different spins) reproduces with enough accuracy the underlying random graph. Then approximate shape of the energetic landscape is recognized by the nonequilibrium model as different spin configurations are explored during the time evolution. As a result, at low enough temperatures the nonequilibrium model can be trapped in a spin configuration, or in a set of similar configurations, close to one of the metastable configurations of the corresponding equilibrium Ising model, which results in the appearance of the SG phase. For small q even weak thermal noise is able to destabilize such configurations, thus the critical temperature for the SG transition T (SG) c,M C is low. As q is increased fine details of the energetic landscape exert effect on the evolution of the nonequilibrium model and T (SG) c,M C increases and approaches that for the corresponding equilibrium Ising model. This interpretation poses a question if the SG transition can be observed in the q-neighbor Ising model on heterogeneous networks, where it is rather impossible to reproduce the structure of the underlying network containing nodes with arbitrarily high degree (hubs) by interactions of each spin only with its finite q-neighborhood. Verification of this possibility requires further extensive MC simulations and is beyond the scope of this paper.

V. SUMMARY AND CONCLUSIONS
In this paper the q-neighbor Ising model on homogeneous random graphs with quenched disorder of FM and AFM exchange interactions associated with the edges of the network was considered. In comparison with the original q-neighbor Ising model with purely FM exchange integrals [15,[24][25][26] the model under study for fixed topology of connections (e.g., the mean degree of nodes k ) and size of the q-neighborhood shows richer critical behavior with varying temperature as the fraction of AFM exchange integrals r is varied. For example, first-and second-order FM phase transitions or second-order FM and SG transitions can occur for different r, and the corresponding critical lines on the T vs. r phase plane meet in TCP. Concerning the FM transition, MFA and homogeneous PA were extended to take into account the effect of AFM exchange interactions on the transition. Quantitative agreement between predictions of the homogeneous PA with results of MC simulations was obtained for a broad range of the model parameters, with noticeable discrepancies observed in the vicinity of the abovementioned TCPs or in the vicinity of r for which the critical temperature approaches zero; in particular, for q comparable with k destabilization of the FM phase with T → 0 predicted by the PA for a certain range of r was not confirmed in MC simulations. Concerning the SG transition, evidence for its occurrence is based solely on MC simulations; this kind of phase transition occurs in the models with q > k /2 for larger values of r than the FM transition.
In the context of modelling opinion formation it should be mentioned that in the q-neighbor Ising model the SG phase can occur in the range of parameters, in particular of the fraction r of the AFM exchange integrals, which is realistic in the models for social interactions; e.g., in the case of the political stage divided between two parties each person can easily interact both with the followers of the same or another party and thus tend to follow or object their opinions. The same is true also in the majority vote model [35]. This suggests that in real societies apart from the spectacular and widely studied FM transition to consensus also a more subtle transition to the SG phase may occur, characterized by local rather than global ordering of agents' opinions. In the context of nonequilibrium models it is interesting to note that phase diagrams obtained from MC simulations of two nonequilibrium counterparts of the Ising model on random graphs with a fraction r of the AFM exchange integrals, the q-neighbor Ising model and the majority vote model [35], are quatitatively similar to each other and resemble those for the equilibrium model for dilute SG [42]: on the T vs. r phase plane there are critical lines for the second-order FM and SG transitions merging in a TCP. Moreover, predictions of the PA in the two above-mentioned nonequilibrium models are also qualitatively similar; in particular, in both cases the PA suggests the possibility of destabilization of the FM phase and occurrence of the PM phase as T → 0. Hence, it seems justified to search for such similarities in other related nonequilibrium models. Taking into account that in the noisy q-voter model on random graphs many results in the MFA and PA can be obtained analytically [13][14][15][16][17][18][19], this model, with a sort of AFM interactions included, could be a good candidate for further studies in the above-mentioned direction.