Tricritical behavior of nonequilibrium Ising spins in fluctuating environments

TL;DR

This paper explores how nonequilibrium conditions in a fluctuating network of Ising spins lead to tricritical behavior, showing transitions between continuous and discontinuous phase changes influenced by temperature differences.

Contribution

It introduces a generalized nonequilibrium Ising model on a fluctuating network, revealing tricritical points and phase transition behaviors not seen in equilibrium or simpler models.

Findings

Identification of a tricritical point in the phase diagram.

Demonstration that nonequilibrium driving can change transition order.

Analysis of heat flow between heat baths.

Abstract

We investigate the phase transitions in a coupled system of Ising spins and a fluctuating network. Each spin interacts with $q$ neighbors through links of the rewiring network. The Ising spins and the network are in thermal contact with the heat baths at temperatures $T_S$ and $T_L$, respectively, so that the whole system is driven out of equilibrium for $T_S \neq T_L$. The model is a generalization of the $q$-neighbor Ising model, which corresponds to the limiting case of $T_L=\infty$. Despite the mean field nature of the interaction, the $q$-neighbor Ising model was shown to display a discontinuous phase transition for $q\geq 4$. Setting up the rate equations for the magnetization and the energy density, we obtain the phase diagram in the $T_S$-$T_L$ parameter space. The phase diagram consists of a ferromagnetic phase and a paramagnetic phase. The two phases are separated by a continuous phase transition belonging to the mean field universality class or by a discontinuous phase transition with an intervening coexistence phase. The equilibrium system with $T_S=T_L$ falls into the former case while the $q$-neighbor Ising model falls into the latter case. At the tricritical point, the system exhibits the mean field tricritical behavior. Our model demonstrates a possibility that a continuous phase transition turns into a discontinuous transition by a nonequilibrium driving. Heat flow induced by the temperature difference between two heat baths is also studied.

Tables

Table 1. Table 1: Quartet configurations with s x = + subscript 𝑠 𝑥 s_{x}=+ and associated energy costs and realization probabilities. For s x = − subscript 𝑠 𝑥 s_{x}=- , ω α − subscript superscript 𝜔 𝛼 \omega^{-}_{\alpha} is the spin reversal of w α + subscript superscript 𝑤 𝛼 w^{+}_{\alpha} and p α − ​ ( m , e ) = p α + ​ ( − m , e ) subscript superscript 𝑝 𝛼 𝑚 𝑒 subscript superscript 𝑝 𝛼 𝑚 𝑒 p^{-}_{\alpha}(m,e)=p^{+}_{\alpha}(-m,e) .

$\alpha$	${\omega }_{\alpha }^{+}$	$e_r$	$e_f$	$e_{\mathrm{rf}}$	$p_{\alpha }^{+}(m,e)$
1	$(++++)$	0	2	2	${\left(1-e+2m\right)}^2/16$
2	$(+++-)$	0	2	2	$\left(1+e\right)\left(1-e+2m\right)/16$
3	$(++-+)$	0	2	-2	$\left(1+e\right)\left(1-e+2m\right)/16$
4	$(++--)$	4	2	-2	$\left(1-e+2m\right)\left(1-e-2m\right)/16$
5	$(+-++)$	0	-2	2	$\left(1+e\right)\left(1-e+2m\right)/16$
6	$(+-+-)$	-4	-2	2	${\left(1+e\right)}^2/16$
7	$(+--+)$	0	-2	-2	${\left(1+e\right)}^2/16$
8	$(+---)$	0	-2	-2	$\left(1+e\right)\left(1-e-2m\right)/16$

Equations

$$P_T(\Delta E)\equiv \min \left[1,e^{-\Delta E/T}\right],$$

$$H(\mathsf{A},\mathbf{s})=-\frac{J}{2}\sum _{i,j}{A}_{i,j}{s}_i{s}_j,$$

$$\sum _j{A}_{i,j}=q$$

$$m=n_{+}-n_{-}.$$

$$n_{\pm }=\frac{1\pm m}{2}.$$

$$e=-(n_{++}+n_{--})+n_{+-}.$$

$$\begin{array}{rl}{n}_{++}& =\frac{1}{4}(1-e+2m),\\ n_{--}& =\frac{1}{4}(1-e-2m),\\ n_{+-}& =\frac{1}{2}(1+e).\end{array}$$

$$f^{+}(m,e,\{n_{\alpha }\})=\frac{q!}{(n_{+}{)}^q}\prod _{\alpha =1}^8\frac{(p_{\alpha }^{+}{)}^{n_{\alpha }}}{(n_{\alpha })!}.$$

$$f^{-}(m,e,\{n_{\alpha }\})=\frac{q!}{(n_{-}{)}^q}\prod _{\alpha =1}^8\frac{(p_{\alpha }^{-}{)}^{n_{\alpha }}}{(n_{\alpha })!}.$$

$$f^{+}(m,e,\{n_{\alpha }\})=f^{-}(-m,e,\{n_{\alpha }\}).$$

$$E_{\mathrm{r}}=\sum _{\alpha }{n}_{\alpha }{e}_{\mathrm{r}}(\alpha ).$$

$$E_{\mathrm{f}}=\sum _{\alpha }{n}_{\alpha }{e}_{\mathrm{f}}(\alpha ),$$

$$E_{\mathrm{rf}}=\sum _{\alpha }{n}_{\alpha }{e}_{\mathrm{rf}}(\alpha ),$$

$$\frac{dm}{dt}=F(m,e)\text{mbox}and\frac{de}{dt}=G(m,e),$$

$$\begin{array}{rl}F(m,e)& =2\sum _{\{n_{\alpha }\}}\left(-n_{+}{f}^{+}+n_{-}{f}^{-}\right)\left[P_{T_L}(E_{\mathrm{r}})P_{T_S}(E_{\mathrm{rf}})+\left\{1-P_{T_L}(E_{\mathrm{r}})\right\}{P}_{T_S}(E_{\mathrm{f}})\right],\\ G(m,e)& =\frac{2}{q}\sum _{\{n_{\alpha }\}}\left(n_{+}{f}^{+}+n_{-}{f}^{-}\right)[E_{\mathrm{r}}{P}_{T_L}(E_{\mathrm{r}})\left\{1-P_{T_S}(E_{\mathrm{rf}})\right\}+(E_{\mathrm{r}}+E_{\mathrm{rf}})P_{T_L}(E_{\mathrm{r}})P_{T_S}(E_{\mathrm{rf}})\\ & +E_{\mathrm{f}}\left\{1-P_{T_L}(E_{\mathrm{r}})\right\}{P}_{T_S}(E_{\mathrm{f}})].\end{array}$$

$$F(m,e)=-F(-m,e)\text{mbox}andG(m,e)=G(-m,e).$$

$$I(m)\equiv F(m,\epsilon (m))=-I(-m).$$

$$L(m)=-{\int }_0^{m}I(m^{\prime })d{m}^{\prime },$$

$$L(m)=\frac{a_2}{2}{m}^2+\frac{a_4}{4}{m}^4+\frac{a_6}{6}{m}^6+O(m^8)$$

$$\frac{n_{\pm }{f}^{\pm }(m,e,\{n_{\alpha }\})}{n_{\pm }{f}^{\pm }(m,e,\{n_{\gamma (\alpha )}\}}=\frac{P_T(-\Delta E_r)}{P_T(\Delta E_r)}=e^{4{\beta }_{eq}(n_4-n_6)}$$

$$(1-e+2m)(1-e-2m)=e^{4{\beta }_{eq}}(1+e{)}^2.$$

$$\frac{n_{\pm }{f}^{\pm }(m,e,\{n_{\alpha }\})}{n_{\mp }{f}^{\mp }(m,e,\{n_{\delta (\alpha )}\})}=\frac{P_T(-\Delta E_f)}{P_T(\Delta E_f)}=e^{2\beta (n_a-n_b)},$$

$${\left(\frac{(1-m)}{(1+m)}\right)}^{2(q-1)}{\left(\frac{1-e+2m}{1-e-2m}\right)}^q=1.$$

$$T_{eq,c}(q)=\frac{2}{\ln (q)-\ln (q-2)}.$$

$$\begin{array}{rl}m& =N^{-\beta /\bar{\nu }}{O}_m((T_{eq}-T_{eq,c})N^{1/\bar{\nu }}),\\ \chi & =N^{\gamma /\bar{\nu }}{O}_{\chi }((T_{eq}-T_{eq,c})N^{1/\bar{\nu }})\end{array}$$

$$I^{\prime }(0)={\left[({\partial }_m+{\epsilon }^{\prime }{\partial }_e)F\right]}_{m=0,e=\epsilon (0)}=0,$$

$${\left[{\partial }_{m}F\right]}_{m=0,e=\epsilon (0)}=0.$$

$${\left[({\partial }_m^{3}F)({\partial }_{e}G)-3({\partial }_m{\partial }_{e}F)({\partial }_m^{2}G)\right]}_{m=0,e=\epsilon (0)}=0.$$

$${\beta }_{L,\mathrm{TCP}}=0.180954\cdots \text{mbox}and{\beta }_{S,\mathrm{TCP}}=0.408461\cdots$$

$$\dot{Q}=\frac{2}{q}\sum _{\{n_{\alpha }\}}\left(n_{+}{f}^{+}+n_{-}{f}^{-}\right)E_{\mathrm{r}}{P}_{T_L}(E_{\mathrm{r}}).$$

Full text

I Introduction

The Ising model is one of the most studied statistical physics systems for the theory of phase transitions and critical phenomena. Recently, Jędrzejewski et al. Jędrzejewski et al. (2015) studied the phase transition in the so-called $q$-neighbor Ising model. In this model, an Ising spin interacts ferromagnetically with $q$ instant neighbors which are chosen randomly among the other spins. The model was shown to undergo a phase transition from a high-temperature paramagnetic phase to a low-temperature ferromagnetic phase. Interestingly, the phase transition is of first order (discontinuous) with a discontinuous jump in the spontaneous magnetization for $q\geq 4$, while it is of second order (continuous) exceptionally at $q=3$.

The $q$-neighbor Ising model looks similar to the Ising model on an annealed network Lee et al. (2009). Suppose that Ising spins are on nodes of a network and interact with each other through links. In the annealed network, links are assumed to be rewired so fast that every spin is connected to all the others with effective coupling strengths. The equilibrium Ising model on the annealed network is described by the mean field (MF) theory and is shown to display the continuous phase transition Lee et al. (2009). In the $q$-neighbor Ising model, where spins interact with random neighbors, spatial correlations are negligible and the MF theory is also exact. Thus one might expect the continuous phase transition as the MF theory predicts. Given the MF nature of the model, the discontinuous transition in the $q$-neighbor Ising model is puzzling.

The purpose of this study is to reveal the reason why the $q$-neighbor Ising model deviates from the equilibrium MF theory prediction. We notice that not only the Ising spins but also the links connecting spins are fluctuating dynamic variables. The $q$-neighbor Ising model will be shown to be a limiting case of a nonequilibrium system driven between two heat baths $B_{L}$ and $B_{S}$ at different temperatures $T_{L}$ and $T_{S}$, respectively. The Ising spins are in thermal contact with the heat bath $B_{S}$, while the links are in thermal contact with $B_{L}$. The $q$-neighbor Ising model corresponds to the case with $T_{L}=\infty$. The nonequilibrium driving with $T_{L}\neq T_{S}$ is responsible for the deviation from the equilibrium MF theory prediction.

Phase transitions in nonequilibrium Ising models have been studied for a long time De Masi et al. (1985); Gonzalez-Miranda et al. (1987); Droz et al. (1990); Bassler and Rácz (1994); Szolnoki (2000); Pleimling et al. (2010); Borchers et al. (2014). Ising spins can be driven out of equilibrium under any dynamics breaking the detailed balance. The nature of resulting nonequilibrium phase transitions may or may not belong to the same universality class as the equilibrium counterpart. The equilibrium Ising universality class is stable against a nonequilibrium driving if the dynamics does not conserve the order parameter Grinstein et al. (1985); Blöte et al. (1990). On the other hand, nonequilibrium Ising models with order parameter conserving dynamics display different types of phase transitions Schmittmann and Zia (1991); Schmittmann (1990); Cheng et al. (1991); Praestgaard et al. (1994); Bassler and Rácz (1994); Borchers et al. (2014). Ising systems with spin-exchange dynamics are such examples. These systems can be driven out of equilibrium by introducing multiple heat baths or a directional bias in the spin exchange process. In addition to the nonequilibrium critical phenomena, energy or particle currents Borchers et al. (2014); Pleimling et al. (2010) and the entropy production Shim et al. (2016) have been attracting growing interests recently.

The Ising spins in our study are connected via fluctuating links at a different temperature. In Sec. II, we introduce a nonequilibrium Ising model involving two heat baths of temperature $T_{S}$ and $T_{L}$. This model includes the $q$-neighbor Ising model as a limiting case. The analytic theory for the model is set up in Sec. III, and the resulting phase diagram in the parameter space of $T_{S}$ and $T_{L}$ is presented in Sec. IV. We find that the ordered phase and the disordered phase are separated by the continuous phase transition line in some region of the parameter space and by the coexistence phase in the other region. The continuous phase transition line ends at the tricritical point. The equilibrium model with $T_{S}=T_{L}$ undergoes the continuous phase transition while the $q$-neighbor Ising model undergoes the discontinuous phase transition through the coexistence phase. We close the paper with summary and discussions on the heat flow in Sec. V.

II Nonequilibrium Ising model

We begin with introducing the $q$-neighbor Ising model of Ref. Jędrzejewski et al. (2015). The system consists of $N$ Ising spins $s_{n}$ ($n=1,2,\cdots,N$) in thermal contact with a heat bath at temperature $T$. The spin states are represented as $s_{n}=\pm 1$ or simply $\pm$. Spin configurations are updated following the Monte Carlo rule. Each time step, one selects a spin $s_{i}$ and $q$ other spins, denoted as $\{s_{i_{k}}|k=1,\ldots,q\}$, at random. These $q$ spins are designated as instant interacting neighbors of $s_{i}$ with the energy function $E(s_{i};\{s_{i_{k}}\})=-Js_{i}\sum_{k=1}^{q}s_{i_{k}}$ with a ferromagnetic coupling constant $J>0$. The spin $s_{i}$ is then flipped ($s_{i}\to-s_{i}$) with the probability

[TABLE]

where $\Delta E=E(-s_{i};\{s_{i_{k}}\})-E(s_{i};\{s_{i_{k}}\})$ is the energy change upon flipping $s_{i}$.

The flipping probability in (1) is taken commonly in the Metropolis algorithm simulating the thermal equilibrium states at a given temperature $T/k_{B}$ with the Boltzmann constant $k_{B}$ Binder and Heermann (2010). The Boltzmann constant will be set to unity hereafter. Thus, the $q$-neighbor Ising model appears to be a thermal equilibrium system of Ising spins interacting with random neighbors. Surprisingly, the $q$-neighbor Ising model exhibits the first-order phase transition for any $q\geq 4$ Jędrzejewski et al. (2015). The result is in sharp contrast to the equilibrium MF theory predicting the continuous phase transition Goldenfeld (1992).

In the $q$-neighbor Ising model, both the Ising spins and the links between interacting spins are fluctuating dynamic variables. The Ising spins interact with the heat bath of temperature $T$. On the other hand, the links are rewired completely randomly. This indicates that two different heat baths, one for the spins and another for the links, are involved in the $q$-neighbor Ising model.

To be more precise, we introduce the Hamiltonian for the whole system including the spins and the links as

[TABLE]

where $A_{i,j}$ is an element of an adjacency matrix $\mathsf{A}$ and $\bm{s}=(s_{1},\ldots,s_{N})$ denotes a spin configuration. The coupling constant $J$ will be set to unity. The adjacency matrix element $A_{i,j}$ takes $1$ if there is a link between $i$ and $j$ and 0 otherwise. As a convention, we set $A_{i,i}=0$ disallowing a self-loop. The adjacency matrix is constrained by the condition

[TABLE]

for all $i$ to ensure that every site has $q$ neighbors. Then, the $q$-neighbor Ising model is equivalent to the combined system of spins and links with the Hamiltonian (2) where the spins are in thermal contact with a heat bath ${B}_{S}$ of temperature $T_{S}=T$ and the links are in thermal contact with another heat bath ${B}_{L}$ of temperature $T_{L}=\infty$.

We now define the generalized model by introducing the following dynamics to the combined system with the Hamiltonian in (2). The link configuration $\mathsf{A}$ and the spin configuration $\bm{s}$ are updated as follows (see Fig. 1): (i) Select a site $x$ at random. Current neighbors of $x$ are denoted as $x_{1},\ldots,x_{q}$ ($A_{x,x_{k}}=1$). One also selects $q$ distinct sites denoted as $y_{1},\cdots,y_{q}$ among all sites but $x$. They are the potential candidates for new neighbors of $x$. For each $y_{k}$, one further selects one of its neighbor $y^{\prime}_{k}$ at random ($A_{y_{k},y^{\prime}_{k}}=1$). (ii) Try to remove existing links between $x$ and $x_{k}$ and between $y_{k}$ and $y^{\prime}_{k}$ ($A_{x,x_{k}}\to 0\mbox{ and }A_{y_{k},y^{\prime}_{k}}\to 0$), and to add new links between $x$ and $y_{k}$ and between $x_{k}$ and $y^{\prime}_{k}$ ($A_{x,y_{k}}\to 1\mbox{ and }A_{x_{k},y^{\prime}_{k}}\to 1$) for all $k=1,\ldots,q$. The link configuration after the rewiring is denoted by $\mathsf{A}^{\prime}$. The rewiring trial is accepted with the probability $P_{T_{L}}(\Delta E)$, where $\Delta E=H(\mathsf{A}^{\prime},\bm{s})-H(\mathsf{A},\bm{s})$ is the energy change upon rewiring with spin configuration $\bm{s}$ being fixed. (iii) The spin $s_{x}$ is then flipped to $-s_{x}$ with the probability $P_{T_{S}}(\Delta E)$, where $\Delta E=H(\mathsf{A}^{\prime\prime},\bm{s}_{x})-H(\mathsf{A}^{\prime\prime},\bm{s})$ is the energy change upon spin flip. Here, $\mathsf{A}^{\prime\prime}$ denotes the adjacency matrix after the rewiring trial that is $\mathsf{A}^{\prime}$ if the rewiring is accepted or $\mathsf{A}$ otherwise, and $\bm{s}_{x}$ denotes the spin configuration with $s_{x}$ being flipped from $\bm{s}$. The time $t$ is measured in unit of Monte Carlo step per site.

We adopt the so-called degree-preserving rewiring scheme in step (ii) Maslov and Sneppen (2002); Noh (2007). This method allows one to rewire the links under the constraint of (3). Trials resulting in self-loops or double-links are rejected. When $T_{L}=\infty$, rewiring trials are always accepted. Thus, our model with $T_{S}=T$ and $T_{L}=\infty$ reduces to the $q$-neighbor Ising model Jędrzejewski et al. (2015). When $T_{S}=T_{L}$, the dynamics satisfies the detailed balance and the whole system is in thermal equilibrium (see the discussion in Sec. IV.1). When $T_{L}=0$, one may think that the Ising spins will be in thermal equilibrium on the quenched network. However, the network keeps evolving even at $T_{L}=0$. Suppose that the network reaches the ground state link configuration to a given spin configuration. When spin flips at finite $T_{S}$, the links are pumped out of the ground state and rewired. Thus, the model with $T_{L}=0$ is different from the Ising model on the quenched network.

III Mean field theory

Link rewiring allows spins to interact with any other spins. Thus, spatial correlations between spins are negligible and the MF theory is a good approximation. In this section, we derive the MF rate equations for the mean magnetization density per site $m\equiv\frac{1}{N}\langle\sum_{i}s_{i}\rangle$ and the mean energy density per link $e\equiv\frac{2}{qN}\langle H(\mathsf{A},\bm{s})\rangle$ taking account of correlations up to nearest neighbors directly connected with links.

We first introduce several notations. Let $n_{+}$ and $n_{-}$ be the fractions of $+$ and $-$ spins, respectively. The magnetization density is given by

[TABLE]

The normalization $n_{+}+n_{-}=1$ yields that

[TABLE]

Let $n_{++}$, $n_{--}$, and $n_{+-}$ be the fractions of links connecting $++$, $--$, and $+-$ spin pairs, respectively, satisfying the normalization $n_{++}+n_{--}+n_{+-}=1$. The energy density per link is given by

[TABLE]

Those fractions satisfy the relations $n_{+}=n_{++}+\frac{1}{2}n_{+-}$ and $n_{-}=n_{--}+\frac{1}{2}n_{+-}$. Thus one can rewrite the fractions in terms of $m$ and $e$ as

[TABLE]

Since $n_{++}\geq 0$ and $n_{--}\geq 0$, $e$ and $m$ are restricted within the range $|m|\leq(1-e)/2$.

Rewiring a single link of a randomly selected site $x$ involves a quartet of four spins $(s_{x}s_{x_{k}}s_{y_{k}}s_{y^{\prime}_{k}})$ (see Fig. 1). A quartet can take one of the $2^{4}=16$ spin configurations. We label the configurations with $s_{x}=+1$ as $\omega^{+}_{\alpha}$ ($\alpha=1,\ldots,8)$, and the spin-reversed configurations as $\omega^{-}_{\alpha}$. These configurations are listed in Table 1. Given $m$ and $e$, the probability that a quartet is in a certain configuration $\omega^{\pm}_{\alpha}$ is given by a function of $m$ and $e$. They will be denoted as $p^{\pm}_{\alpha}(m,e)$. For example, the quartet $\omega^{-}_{8}=(-+++)$ has the probability $p^{-}_{8}=\frac{1}{2}n_{+-}n_{++}=\frac{1}{16}(1+e)(1-e+2m)$. The quartet probabilities are summarized in Table 1. It is obvious that $p^{+}_{\alpha}(m,e)=p^{-}_{\alpha}(-m,e)$ and $\sum_{\alpha}p^{\pm}_{\alpha}(m,e)=n_{\pm}=(1\pm m)/2$.

It would cost energy $e_{\rm r}(\alpha)=\left[-(s_{x}s_{y_{k}}+s_{x_{k}}s_{y^{\prime}_{k}})+(s_{x}s_{x_{k}}+s_{y_{k}}s_{y^{\prime}_{k}})\right]_{\omega^{\pm}_{\alpha}}$ for rewiring, $e_{\rm f}(\alpha)=\left[2s_{x}s_{x_{k}}\right]_{\omega^{\pm}_{\alpha}}$ for flipping of $s_{x}$ without rewiring, and $e_{\rm rf}(\alpha)=\left[2s_{x}s_{y_{k}}\right]_{\omega^{\pm}_{\alpha}}$ for flipping of $s_{x}$ after rewiring. The energy costs are summarized in Table 1. Due to the spin-reversal symmetry of the Hamiltonian, the energy costs for the configurations $\omega^{+}_{\alpha}$ and $\omega^{-}_{\alpha}$ are the same.

Monte Carlo dynamics involves $q$ quartets sharing a randomly selected site $x$. We denote the number of quartets of configuration $\omega^{\pm}_{\alpha}$ as $n_{\alpha}$ $(=0,1,\ldots,q$). Due to the spin-reversal symmetry, we do not need to count the number of quartets with $s_{x}=+1$ and $s_{x}=-1$ separately. They are constrained by the sum rule $\sum_{\alpha=1}^{8}n_{\alpha}=q$.

We are ready to set up the rate equations for $m$ and $e$. Suppose that a site $x$ is selected at random. Provided that $s_{x}=+$, the probability that $q$ associated quartets are specified by $\{n_{\alpha}\}$ is given by

[TABLE]

The term $(n_{+})^{q}=(\sum_{\alpha}p^{+}_{\alpha})^{q}$ in the denominator guarantees the normalization $\sum_{\{n_{\alpha}\}}f^{+}(\{n_{\alpha}\})=1$, where the summation $\sum_{\{n_{\alpha}\}}$ is over all sequences of non-negative integers $\{n_{\alpha}\}$ satisfying $\sum_{\alpha}n_{\alpha}=q$. Similarly, the probability for $\{n_{\alpha}\}$ with $s_{x}=-1$ is given by

[TABLE]

The symmetry property $p^{+}_{\alpha}(m,e)=p^{-}_{\alpha}(-m,e)$ yields that

[TABLE]

The updating probabilities of links and spins are determined by the associated energy changes. The link rewiring would cost

[TABLE]

The spin flip would cost

[TABLE]

without link rewiring with probability $1-P_{T_{L}}(E_{\rm r})$, or

[TABLE]

after rewiring with probability $P_{T_{L}}(E_{\rm r})$.

Combining all the quantities, we finally obtain the rate equations in the $N\to\infty$ limit as

[TABLE]

where

[TABLE]

Here, $P_{T}(E)$ is the transition probability function defined in (1) and the factor $(\frac{2}{q})$ of $G$ accounts for the link density. The dependence on $T_{L}$, $T_{S}$, and $q$ is not shown explicitly. Using the relation in (10), one finds that

[TABLE]

Note that in the $T_{L}\rightarrow\infty$ limit, the function $F(m,e)$ becomes independent of $e$ and one recovers the rate equation of Ref. Jędrzejewski et al. (2015).

IV Phase diagram

The steady-state phase diagram is determined by analyzing the fixed point solution of the rate equation in (14). Firstly, Fig. 2 demonstrates how the fixed points bifurcate as $T_{S}$ varies with fixed $T_{L}=10$ at $q=4$. When $T_{S}>T_{d1}$ with a threshold temperature $T_{d1}$, the system has a single stable fixed point at $m=0$ (see Fig. 2(a)) and is in a disordered paramagnetic phase. When $T_{d2}<T_{S}<T_{d1}$ with another threshold temperature $T_{d2}$, two pairs of stable and unstable fixed points with $|m|\neq 0$ appear additionally (see Fig. 2(b)). Hence, the system can coexist in the paramagnetic phase and in the ordered ferromagnetic phase. When $T_{S}<T_{d2}$, the fixed point at $m=0$ becomes unstable after merging with the unstable fixed points (see Fig. 2(c)). The system is in the ferromagnetic phase with nonzero spontaneous magnetization $m_{0}\equiv|m|$. In Fig. 2(d), we draw the steady state values of $e$ and $m_{0}$ against $T_{S}$. The system undergoes a first order transition with the intermediate coexistence region. This behavior is similar to that of the $q$-neighbor model with $T_{L}=\infty$ Jędrzejewski et al. (2015).

The discontinuous transition is confirmed with the Monte Carlo simulations. We have performed the simulations in two different setups. In the cooling (heating) setup, we increase (decrease) the inverse temperature $\beta_{S}$ by $0.01$ in every 2000 time steps. The Monte Carlo simulation data are presented in Fig. 2(d). The numerical data exhibit the hysteresis behavior which is characteristic of discontinuous phase transitions.

When we lower the temperature $T_{L}$, a qualitatively different behavior emerges. Figure 3 shows the evolution of the fixed points as $T_{S}$ varies with fixed $T_{L}=10/3$. When $T_{S}>T_{c}$ with a critical threshold temperature $T_{c}$, there is a single stable fixed point at $m=0$ (see Fig. 3(a)). As $T_{S}$ decreases below $T_{c}$, the fixed point at $m=0$ becomes unstable, while two stable fixed points with $m_{0}=|m|\neq 0$ appear near the unstable fixed point (see Figs. 3(b) and (c)). Hence, the spontaneous magnetization $m_{0}$ and the energy density $e$ vary continuously and the system undergoes a continuous phase transition (see Fig. 3(d)).

The nature of the phase transition can be studied systematically. Let $\varepsilon(m)$ denote the nullcline satisfying $G(m,\varepsilon(m))=0$. The symmetry property $G(-m,e)=G(m,e)$ implies that the function is even in $m$, $\varepsilon(-m)=\varepsilon(m)$. The fixed points of the rate equation are found from zeroes of

[TABLE]

It is convenient to consider

[TABLE]

which is even in $m$. The stable fixed points of the rate equations correspond to the local minima of $L(m)$. Hence, regarding $L(m)$ as the Landau free energy, we can apply the phenomenological Landau theory Goldenfeld (1992). Note, however, that $L(m)$ is not the real free energy because the system is not in thermal equilibrium.

One can expand the Landau free energy as

[TABLE]

with $T_{L}$ and $T_{S}$ dependent coefficients $a_{n}$. The paramagnetic fixed point at $m=0$ is stable when $a_{2}>0$ and unstable when $a_{2}<0$. Thus, the threshold for the paramagnetic state is determined by the condition $a_{2}=0$. If $a_{4}$ is positive near the threshold, the spontaneous magnetization scales as $m_{0}\simeq(-a_{2}/a_{4})^{1/2}$ and the system undergoes a continuous phase transition. Figure 3 exemplifies this case. On the other hand, if $a_{4}$ is negative near the threshold, the system is bistable with $m_{0}=0$ and $m_{0}\simeq[(-a_{4}+\sqrt{a_{4}^{2}-4a_{2}a_{6}})/(2a_{6})]^{1/2}$ in the region $a_{4}^{2}-4a_{2}a_{6}\gtrsim 0$. The spontaneous magnetization jumps from zero to $m_{0}\propto(-a_{4}/a_{6})^{1/2}$. Hence the system undergoes a discontinuous transition from the paramagnetic phase to the ferromagnetic phase separated by the coexistence phase, as exemplified in Fig. 2. The tricritical point is located at the point where $a_{2}=a_{4}=0$.

We present the phase diagram for the system with $q=4$ in Fig. 4. The phase diagram consists of three phases: the paramagnetic (P) phase, the ferromagnetic (F) phase, and the coexistence (C) phase. The phase diagram is constructed as follows. We first draw the lines $a_{2}=0$ and $a_{4}=0$. These lines are found numerically easily since we know the analytic expressions for $I(m)$ and $L(m)$. The two lines intersect with each other at the tricritical point (TCP). The line $a_{2}=0$ with $a_{4}>0$ is the boundary between the F and the P phases, while the line $a_{2}=0$ with $a_{4}<0$ is the boundary between the F and the C phases. The boundary between the P and the C phases, which can be approximated by the line $a_{4}^{2}=4a_{2}a_{6}$ neglecting $O(m^{8})$ term in (19), is located numerically by examining the existence of the local minimum of $L(m)$ at $m\neq 0$.

IV.1 Equilibrium case with $T_{S}=T_{L}=T_{eq}$

In order to reconcile with the results of the equilibrium Ising model on the annealed network Lee et al. (2009), we consider the equilibrium line where $T_{L}=T_{S}=T_{eq}=1/\beta_{eq}$ in detail. We can show that the transition probabilities in the rate equation satisfies the detailed balance (DB) condition.

First, consider the rewiring process which transforms each quartet configuration $\alpha=1,2,3,4,5,6,7,8$ to $\gamma=\gamma(\alpha)=1,2,5,6,3,4,7,8$, respectively. The DB requires that

[TABLE]

with $\Delta E_{r}=4(n_{4}-n_{6})$ (see Table 1). Using $p^{\pm}_{3}=p^{\pm}_{5}$ and $p^{\pm}_{4}/p^{\pm}_{6}=(1-e+2m)(1-e-2m)/(1+e)^{2}$, we find that the relation holds for all $\{n_{\alpha}\}$ if

[TABLE]

Secondly, consider the spin flip process which transforms each quartet configuration $\alpha=1,2,3,4,5,6,7,8$ to $\delta=\delta(\alpha)=8,7,6,5,4,3,2,1$, respectively. The DB requires that

[TABLE]

where $\Delta E_{f}=2(n_{a}-n_{b})$ with $n_{a}=(n_{1}+n_{2}+n_{3}+n_{4})$ and $n_{b}=(n_{5}+n_{6}+n_{7}+n_{8})$ (see Table 1). Using the expressions for $p^{+}_{\alpha}(m,e)=p^{-}_{\alpha}(-m,e)$ in Table 1 and (21), we find that the relations holds for all $\{n_{\alpha}\}$ if

[TABLE]

One can show further that the DB is also satisfied under the simultaneous rewiring and flipping by combining the calculations for each. Therefore, when $T_{L}=T_{S}$, the transition rates satisfy the DB condition and the equilibrium energy density and the magnetization are determined by (21) and (23).

We add a remark on the DB. Although the DB condition is satisfied at the rate equation level, it is not satisfied at the microscopic level of the Monte Carlo dynamics where the link rewiring and the spin flipping are tried subsequently. One can show that the rewiring and flipping do not commute with each other, which breaks the DB. Nevertheless, the preceding paragraphs show that the DB is satisfied in the average sense. Thus we will regard the model with $T_{S}=T_{L}$ as the equilibrium model.

As seen from the phase diagram in Fig. 4, the equilibrium system undergoes the continuous phase transition. The transition temperature $T_{eq,c}$ is found by analyzing (21) and (23). After a straightforward algebra, we obtain that

[TABLE]

The spontaneous magnetization behaves as $m_{0}\sim|T_{eq}-T_{eq,c}|^{\beta}$ with the MF exponent $\beta=1/2$. We have also performed the Monte Carlo simulations to measure the other critical exponents. Figure 5 shows the finite size scaling plots for the magnetization $m=\langle|\sum_{i}s_{i}|\rangle/N$ and the susceptibility $\chi=(\langle(\sum_{i}s_{i})^{2}\rangle-\langle|\sum_{i}s_{i}|\rangle^{2})/(NT_{eq})$. Near the critical point, they follow the scaling form

[TABLE]

with scaling functions $O_{m}$ and $O_{\chi}$ and the MF critical exponents $\bar{\nu}=2$ and $\gamma=1$. The data collapse confirms that the equilibrium model belongs to the MF universality class. This result suggests that the discontinuous phase transition in the $q$-neighbor model is the effect of the nonequilibrium driving.

IV.2 Tricritical point

The tricritical point TCP lies at the point where $a_{2}=a_{4}=0$ in (19). The first condition $a_{2}=0$ yields that

[TABLE]

where ′ denotes the derivative with respect to $m$ and $\partial_{m,e}$ is a shorthand notation for the partial differentiation. Note that $\varepsilon(m)$ is an even function of $m$, hence $\varepsilon^{\prime}(0)=0$. Thus, we obtain the condition

[TABLE]

The second condition $a_{4}=0$ requires that $I^{\prime\prime\prime}(0)=0$. Taking the derivatives and using $\varepsilon^{\prime}(0)=\varepsilon^{\prime\prime\prime}(0)=0$, one obtains that $(\partial_{m}^{3}+3\varepsilon^{\prime\prime}\partial_{m}\partial_{e})F(0,\varepsilon(0))=0$. The function $\varepsilon(m)$ is defined by the relation $G(m,\varepsilon(m))=0$, which yields that $\varepsilon^{\prime\prime}(0)=-(\partial_{m}^{2}G(0,\varepsilon(0)))/(\partial_{e}G(0,\varepsilon(0)))$. Thus, we obtain that

[TABLE]

By solving (27) and (28), we find the tricritical point is located at

[TABLE]

for $q=4$. The location of the TCP can be found numerically exactly for any $q\geq 4$.

The order parameter $m_{0}$ follows the tricritical scaling behavior near the TCP. In Fig. 6, we plot the spontaneous magnetization $m_{0}$ along the line $a_{4}=0$ shown in Fig. 4. It scales as $m_{0}\sim(\beta_{S}-\beta_{S,{\rm TCP}})^{1/4}$ with the MF tricritical exponent $1/4$ instead of $1/2$ Goldenfeld (1992).

At $q=3$, the lines $a_{2}=0$ and $a_{4}=0$ do not meet in the $(\beta_{S},\beta_{L})$ space, and $a_{4}>0$ along the line $a_{2}=0$. Thus the transition is always continuous and the tricritical point is absent.

V Summary and discussions

We have studied the phase transitions in the Ising spin system on the link-rewiring network. The system is in contact with two heat baths ${B}_{S}$ and ${B}_{L}$ that govern the thermal fluctuations of the spins and the links, respectively. This model is introduced in order to explain the discontinuous phase transition recently reported in the $q$-neighbor Ising model where Ising spins interact with random neighbors Jędrzejewski et al. (2015). Such a result was puzzling since the MF theory is working in the $q$-neighbor Ising model and the equilibrium Ising model in the MF theory exhibits a continuous phase transition. We have found that the $q$-neighbor Ising model is indeed a nonequilibrium system driven between two heat baths, ${B}_{S}$ for spins at finite temperatures and ${B}_{L}$ for links at the infinite temperature. We have constructed the phase diagram of the extended model in the parameter space of $\beta_{S}=1/T_{S}$ and $\beta_{L}=1/T_{L}$ with the temperatures $T_{S}$ and $T_{L}$ for spins and links, respectively. When $T_{S}=T_{L}$, the model reduces to the equilibrium model and displays the continuous phase transition belonging to the equilibrium MF Ising universality class. When $T_{L}$ is much larger than $T_{S}$, the coexistence phase emerges and the system exhibits the discontinuous phase transition. The coexistence phase terminates at the tricritical point. Our result shows that the nonequilibrium driving can change the nature of the phase transition from being continuous to being discontinuous.

A thermal system in between two heat baths at different temperatures conducts heat from a high temperature bath to a low temperature one. The steady-state heat flux, average heat flow per unit time, from the bath ${B}_{L}$ to the system will be denoted as $\dot{Q}$. The steady-state heat flux from the bath ${B}_{S}$ to the system is then equal to $-\dot{Q}$. The heat flow results in the increase of the total entropy with the rate $\dot{S}=\dot{Q}(-\beta_{L}+\beta_{S})$. Recently, the critical scaling behavior of the entropy production near the nonequilibrium phase transition has been studied Shim et al. (2016). The heat is injected into the system from the bath ${B}_{L}$ when links are rewired. Hence, by modifying (15), one finds that the heat flux per link is written as

[TABLE]

The heat flux vanishes in the equilibrium case with $T_{L}=T_{S}$ due to the detailed balance thereon.

We investigate the heat flow for the $q$-neighbor Ising model with $T_{L}=\infty$, where the expression is simplified to

[TABLE]

This expression is understood intuitively. Consider the rewiring of a single quartet. There are two links in a quartet, and the average energy of a quartet before rewiring is $2e$. After rewiring to random neighbors, the average energy becomes $-2m^{2}$. Thus, the heat flux should be given by (31).

The heat flux, evaluated from the fixed point solutions for $e$ and $m$, is presented in Fig. 7. The nonzero positive heat flux confirms that the $q$-neighbor Ising model is indeed out of equilibrium. It varies continuously at $q=3$ and discontinuously at $q=4$ as the order parameter $m_{0}$ does. It is noteworthy that the heat flux in Fig. 7 increases (decreases) as $\beta_{S}-\beta_{L}$ increase in the paramagnetic (ferromagnetic) state. The heat flux usually increases as the temperature difference becomes large. In this regard, the decrease of $\dot{Q}$ in the ferromagnetic state is odd. We also leave it for a future work to understand the peculiar behavior.

Acknowledgements.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2016R1A2B2013972).