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 ferromagnetic 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 extended 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.


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-neighbor 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 a e-mail: andrzej.krawiecki@pw.edu.pl (corresponding author) 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 majority-vote 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 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 heterogeneity 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 abovementioned 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 disappear- ance 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 longrange ordering (hence, absence of global consensus) and presence of only short-range 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 second-order 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 symmetric 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 qneighbor 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 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 Sect. 2, the model is introduced and the main results concerning the observed phase transitions are outlined. In Sect. 3, the MFA and PA for the qneighbor 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 Sect. 4, 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 firstand second-order 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 Sect. 4, 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. Unfortunately, the MFA and PA introduced in this paper do not predict SG transition in the q-neighbor Ising model under study, thus evidence for this transition is based solely on MC simulations. Finally, Sect. 5 is devoted to summary and conclusions.

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 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 (2) (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 use of the spin-flip rate given by Eq. (5) distinguishes the q-neighbor Ising model from other related models, e.g., from the (noisy) q-voter models in which rate of flips of the spin s j is affected by the qneighborhood only if all its members are unanimous and have orientation opposite to s j . 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].

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 Sect. 4. Unfortunately, the SG transition which is also observed in MC simulations of the q-neighbor Ising model cannot be analyzed within these approaches, as briefly discussed below.

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 qneighbor Ising model on a fully connected graph with uniform FM exchange integrals [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 inconsistent 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 inconsistent 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 Eqs. (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 transition 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 MFA (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 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 C 4 (T, r(T ), q) = 0 for T =T MFA . As pointed out above, the MFA cannot be used to describe the SG transition in the model under study. This is mainly because 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 where the only order parameter is the magnetization m.
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 MFA presented above.

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 NP (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 Eqs. (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 Sect. 3.1. 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 different spins in the q-neighborhood of a given node 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 ν, 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 abovementioned tips is N k P (ν, ν ). 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 NP (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 ν, 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 Eqs. (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 (22) 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) − 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 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], 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.

Stability and bifurcations of fixed points of dynamical equations in the pair approximation
The fixed point(s) of the system of equations (25,26) Stability of the PM fixed point can be determined from the eigenvalues of the Jacobian matrix of the right-hand sides of Eqs. (25), (26) After some calculations it can be found that ∂A ∂b c=1/2,b=θ = ∂B ∂c c=1/2,b=θ = 0; (29) 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 Eqs. (28), (30) and (31)  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 PA ( k , q) the above-mentioned T c = 0 and for r < r PA ( k , q) the FM transition occurs with decreasing T which can be second-or firstorder, 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 c,PA 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 T = T (FM) c2,PA > T c , which can be determined only numerically. For T c < T < T (FM) c2,PA 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 PA ( k , q) at the critical temperature T c = T (FM) c,PA . 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 PA ( k , q) there is T c > 0; for fixed k the value of r PA ( k , q) weakly depends on q, but the critical temperature T c at r = r PA ( k , q) increases with q. Besides, for a range of r below r PA ( 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  c,MC for given q obtained from simulations (Fig. 1). Thus, the above-mentioned 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 (FM) c,PA (T, r) is related to the de Almeida-Thouless line determining the lower border of stability of the FM phase in the replica-symmetric 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 temperaturelike model parameter is predicted by the PA in the majority-vote model with FM and AFM interactions [35], even with a more complex structure (first-order 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 PA ( k , q) for any q, and the numerically obtained loca-tion 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].
As pointed out above the homogeneous PA introduced in Sect. 3.2 does not predict occurrence of the SG transition in the q-neighbor Ising model. Such transition is second-order and is expected to appear for q comparable with k , fixed r > r PA ( k , q) and decreasing temperature. The SG phase is characterized by zero magnetization, i.e., by c = 1/2 as in the case of the PM phase, and occurrence of local ordering of spins, e.g., for r = 1 (purely AFM interactions) the nearest-neighbor spins should have mostly opposite orientations. Hence, in the framework of the PA the SG transition should correspond to the appearance of a stable fixed point of the system of equations (25,26) with c = 1/2 and the concentration of active links b = θ smaller than that in the PM phase. Taking into account the diagonal structure of the Jacobian in the PM fixed point of the system of equations (25,26) at the SG transition temperature the PM fixed point should lose stability due to the eigenvalue λ 2 (31) crossing zero, with λ 1 (30) negative. However, then the PM fixed point would become a hyperbolic unstable fixed point and just below the transition temperature a pair of stable fixed points with c = 1/2 would occur via a supercritical pitchfork bifurcation, one with b = θ smaller and the other one with b = θ larger than in the PM phase, the latter with no clear physical interpretation. In fact, numerical analysis of Eqs. (28,30,31) reveals that for r > r PA ( k , q) both eigenvalues λ 1 , λ 2 remain negative with temperature decreasing to zero, thus the PM point remains stable and no SG transition occurs. Nevertheless, the concentration of active links b = θ corresponding to the PM fixed point decreases significantly with T → 0 which suggests that some correlation between orientations of interacting spins appear.

Results and discussion
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 Sect. 3. 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 c,MC 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 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) for fixed r as functions of T in the vicinity of the expected SG transition. 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   (Fig. 2). It seems that the width of the hysteresis loop obtained from simulations is significantly smaller than that predicted by the MFA or PA. However, the width increases with the number of nodes in the network, with T  c2,MC 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 Fig. 3. According to the MFA and homogeneous PA with increasing r the width of the hysteresis loop associated with the first-order FM transition decreases to zero and eventually the FM transition becomes second-order. and monotonically decrease with T , and the magnetization M changes smoothly at the critical point (Fig. 4b).
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. However, as r approaches r PA the critical temperature T c,PA 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  different N (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].    [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,MC 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,MC 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.
Failure of the PA to describe quantitatively the FM transition for large r as well to predict the SG transition is probably due to violation of its assumptions mentioned in Sect. 3.2. For example, orientations of spins belonging to the q-neighborhood of a given node can be correlated, e.g., if they are connected by an edge; however, in the case of ERGs and RRGs probability that two neighbors of a given node are mutual neighbors (clustering coefficient) is small [2]. Another possibility is that the activity of the links, i.e., the relative orientations of the connected spins, is correlated with the signs of the exchange integrals associated with them. Moreover, the latter correlation may depend on r. Verification of these hypotheses also requires further numerical studies.

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 above-mentioned 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, it is not predicted by the proposed homogeneous PA and 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. Explanation of the origin of the above-mentioned discrepancies between results of MC simulations and predictions of the PA requires further investigation.
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 quantitatively 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 secondorder 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 abovementioned direction.

Author contributions
The author (A.K.) is responsible for the whole content of the paper.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This paper is theoretical and has no associated experimental data, and all numerical results are shown in the figures.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.