Stochastic Two-temperature Nonequilibrium Ising model

TL;DR

This study explores the nonequilibrium stationary state of a 2D Ising model with stochastic temperature modulation, revealing non-monotonic responses and a Boltzmann-like state with an effective temperature, highlighting intrinsic nonequilibrium features.

Contribution

It introduces a stochastic dichotomous temperature modulation in the Ising model and analyzes its effects on NESS, providing a renewal approach and dynamical response theory for the first time.

Findings

Magnetization and energy show non-monotonic dependence on switching rate

In the slow-switching limit, a renewal approach explains the dynamics

A finite energy current confirms the nonequilibrium nature

Abstract

We investigate the nonequilibrium stationary state (NESS) of the two-dimensional Ising model under a stochastic dichotomous modulation of temperature, which alternates between $T_c \pm \delta$ around the critical temperature $T_c$ at a rate $\gamma$. Both magnetization and energy exhibit non-monotonic dependence on $\gamma$, explained by a renewal approach in the slow-switching limit, while for small $\delta$ dynamical response theory quantitatively captures the $\gamma$-dependence of the observables. In the fast-switching regime, the NESS appears Boltzmann-like with a $\gamma$-dependent effective temperature. However, a finite energy current flowing through the system from hot to cold reservoir confirms the intrinsic nonequilibrium nature of the dynamics.

Equations

$$H(C)=-\sum _{⟨ij⟩}{s}_i{s}_j.$$

$$w_{\sigma }(\Delta E)=\frac{{\psi }_0}{2}[(1+\sigma )e^{-{\beta }_{+}\Delta E/2}+(1-\sigma )e^{-{\beta }_{-}\Delta E/2}],$$

$$M=\frac{1}{N}\left|\sum _{i=1}^N{s}_i\right|,E=\frac{1}{2N}\sum _{⟨ij⟩}{s}_i{s}_j.$$

$$⟨M{⟩}_{\delta }-M_c={\chi }_2^M(\gamma ){\delta }^2+O({\delta }^4),\text{and}⟨E{⟩}_{\delta }-E_c={\chi }_2^E(\gamma ){\delta }^2+O({\delta }^4).$$

$$T_{\text{eff}}=T_c(1-\frac{{\delta }^2}{T_c^2}).$$

$$P_{\gamma }(C,t|C_0)$$

$$P_{\gamma }(C,t)\simeq \frac{1}{2}{e}^{-\gamma t}[{\mathbb{P}}_{T_{+}}(C,t)+{\mathbb{P}}_{T_{+}}(C,t)]+\frac{\gamma }{2}\underset{0}{\overset{t}{\int }}d{t}^{\prime }{e}^{-\gamma t^{\prime }}[{\mathbb{P}}_{T_{+}}(C,t^{\prime })+{\mathbb{P}}_{T_{-}}(C,t^{\prime })],$$

$$⟨O{⟩}_{\gamma }$$

$$⟨O(t){⟩}_{\mathrm{eq}}^{\sigma }=⟨O{⟩}_{\mathrm{eq}}^{\sigma }+b_{\sigma }{e}^{-t/{\tau }_{\sigma }},$$

$$⟨O{⟩}_{\gamma }$$

$$g_O^{\sigma }(t)=⟨O(t){⟩}_{\mathrm{eq}}^{\sigma }-⟨O{⟩}_{\mathrm{eq}}^{\sigma }$$

$$⟨O{⟩}_{\delta }=⟨O{⟩}_{\mathrm{eq}}+{\delta }^2{\chi }_2+\mathcal{O}({\delta }^4),$$

$${\chi }_2(\omega )=\frac{1}{2}{⟨O(\omega )\left[-D^{\prime \prime }(\omega )+\frac{1}{2}{S}^{\prime \prime }+(D^{\prime }{)}^2(\omega )+\frac{1}{4}(S^{\prime }{)}^2(\omega )-D^{\prime }(\omega )S^{\prime }(\omega )\right]⟩}_{\mathrm{eq}},$$

$$S(\omega )$$

$$D(\omega )$$

$$⟨O(t){⟩}_{\delta }-⟨O{⟩}_{\mathrm{eq}}={\delta }^2[\frac{1}{2}⟨S^{\prime \prime }O(t){⟩}_{\mathrm{eq}}-⟨S^{\prime }{D}^{\prime }O(t){⟩}_{\mathrm{eq}}].$$

$$k(\{C,\sigma \},\{C^{\prime },\sigma \})$$

$$k(\{C,\sigma \},\{C,{\sigma }^{\prime }\})$$

$$S(\omega )$$

$$D(\omega )$$

$$S^{\prime }(t)$$

$$S^{\prime \prime }(t)$$

$$D^{\prime }(t)$$

$${\chi }_2(t)$$

$${\chi }_2^{(1)}(t)=-\frac{1}{T_c^3}⟨U(t)O(t){⟩}_{\mathrm{eq}},$$

$${\chi }_2^{(2)}(t)=-\frac{{\psi }_0}{2T_c^4}⟨O(t){\mathcal{J}}_E(t)\underset{0}{\overset{t}{\int }}d{t}^{\prime }\sigma (t^{\prime })\sum _{k=1}^{L^2}\Delta {\epsilon }_k(t^{\prime })e^{-{\beta }_c\Delta {\epsilon }_k(t^{\prime })/2}{⟩}_{\mathrm{eq}}$$

$${\chi }_2^{(1)}(t)=-\frac{1}{T_c^3}[⟨H(t)O(t){⟩}_{\mathrm{eq}}-⟨H(0)O(t){⟩}_{\mathrm{eq}}].$$

$${\chi }_2^{(1)}=\frac{1}{T_c}⟨V;O{⟩}_{\mathrm{eq}},$$

$${\chi }_2^{(2)}({\gamma }_O^{\ast })=-{\chi }_2^{(1)}.$$

$${\chi }_2^{(2)}(t)\stackrel{\gamma \to 0}{\to }-\frac{{\psi }_0}{2T_c^4}\sum _{{\sigma }_0}\rho ({\sigma }_0){\sigma }_0^2⟨O(t)U(t){\int }_0^{t}d{t}^{\prime }\sum _{k=1}^{L^2}\Delta {\epsilon }_k(t^{\prime })e^{-{\beta }_c{\epsilon }_k(t^{\prime })/2}{⟩}_{\mathrm{eq}}.$$

$$⟨O{⟩}_{\delta }\stackrel{\gamma \to 0}{\to }⟨O{⟩}_{\mathrm{eq}}-\frac{1}{T_c^3}⟨O(t)U(t){⟩}_{\mathrm{eq}}-\frac{{\psi }_0}{2T_c^4}⟨O(t)U(t)\underset{0}{\overset{t}{\int }}d{t}^{\prime }\sum _{k=1}^{L^2}\Delta {\epsilon }_k(t^{\prime })e^{-{\beta }_c{\epsilon }_k(t^{\prime })/2}{⟩}_{\mathrm{eq}},$$

$$w_{\mathrm{eff}}(\Delta E)$$

Full text

1 Introduction

The Ising model has served as a cornerstone of statistical physics for nearly a century [1, 2]. Originally introduced to understand critical phenomena [3], it has been extensively studied using various dynamical rules [4, 5, 6] that preserve detailed balance with respect to the equilibrium Boltzmann distribution. Despite its apparent simplicity, the Ising model exhibits a wide range of nontrivial phenomena, including phase transitions, domain coarsening, phase separation, and giant fluctuations. Owing to this rich behavior, its influence extends far beyond statistical physics, with applications in diverse fields such as biology, sociology, and economics [7, 8, 9].

The Ising model also provides a natural framework for exploring out-of-equilibrium behaviour. Over the years, several nonequilibrium versions of the model have been explored in different contexts. A commonly studied setup involves coupling different regions of the system to thermal baths at different temperatures. This naturally drives currents and leads to nontrivial steady-state profiles of magnetization and energy. Such arrangements have been investigated in both one [10, 11, 12, 13, 14, 15, 16] and two dimensions [17]. Generalisations involving spatially varying couplings or temperature fields have also been investigated [18, 19]. Another direction focuses on models that break detailed balance at the microscopic level—typically by introducing spin-flip rates that depend on two or more temperatures [20]. These systems have been examined using superstatistics [21], and comparisons across different update schemes (Metropolis, Glauber, Swendsen-Wang) suggesting that certain signatures of equilibrium criticality may persist even far from equilibrium [22]. In related studies, the interplay between finite and infinite temperature baths has been shown to induce interesting steady-state features, such as nonmonotonic variation of surface tension [23]. More recently, time-dependent temperature protocols have attracted attention, where the Ising system is alternately coupled to thermal baths at different temperatures. These studies range from exploring calorimetric response [24] and entropy production rate [25] in mean-field models to exploring applications such as protein sequence inference [26]. However, the nonequilibrium stationary behaviour of the Ising model under such stochastic temporal driving remains mostly unexplored.

In this paper we investigate a scenario where the two-dimensional Ising model is subject to a stochastic dichotomous temperature. In particular, we consider the case where the temperature intermittently switches between two values $T_{c}\pm\delta$, with a constant rate $\gamma$. We characterize the resulting NESS by measuring the mean and fluctuations of magnetization and energy. Surprisingly, for large $\delta$, the average magnetization and energy both show a non-monotonic behaviour as $\gamma$ is increased. We show that this non-monotonicity originates from the difference in the typical relaxation times of the equilibrium Ising model below and above the critical temperature $T_{c}.$ On the other hand for small $\delta$, we characterise the NESS using dynamical response formalism and show that the second-order response coefficient changes its sign as $\gamma$ is increased. Numerical evidence further suggests that, in the large $\gamma$ regime, the stationary weight of the spin configurations admits a Boltzmann-like form with a $\gamma$ and $\delta$ dependent effective temperature $T_{\text{eff}}<T_{c}$. However the dynamics remains nonequilibrium, which is manifest from a persistent flow of energy current through the system.

The paper is organised as follows. In the next section we define the model and present a brief summary of our main results. Section 3 focuses on the behaviour of magnetization and energy in the stationary state. The effective temperature picture in the large $\gamma$ regime is discussed in Sec. 4. We conclude with some general remarks in Sec. 6

2 Model and Results

We consider a periodic square lattice of size $N=L\times L$ where each site $i$ is associated with a classical spin $s_{i}=\pm 1$. The spins interact via the nearest neighbour Ising Hamiltonian, where the energy of a configuration $C=\{s_{1},s_{2},\cdots s_{N}\}$ is given by,

[TABLE]

The system is alternately coupled to two thermal reservoirs at temperatures $T_{c}-\delta$ and $T_{c}+\delta$, and the switching between them is governed by a dichotomous stochastic process with constant rate $\gamma$ [see Fig. 1]. The configuration of the system evolves via a Glauber-like local detail balanced single spin-flip dynamics where a randomly chosen spin flips with a rate,

[TABLE]

where $\Delta E$ is the change in the energy due to the proposed spin flip and $\beta_{\sigma}\equiv 1/T_{\sigma}=1/(T_{c}+\sigma\delta)$ denotes the inverse temperature with $\sigma=\pm 1$. The intermittent change in the temperature is incorporated via the dichotomous process $\sigma\to-\sigma$ with rate $\gamma$. Note that, $\psi_{0}$ is taken to be a constant, independent of temperature, and characterises the inverse of the typical time-scale of spin flips in the system. Thus for a given $\sigma$ the system evolves via ordinary Glauber dynamics at a temperature $T_{\sigma}$.

For $\delta=0$ the system is equivalent to the ordinary Ising model at a constant temperature $T_{c}$. In that case, for any finite size $L$, the system reaches an equilibrium state with average magnetization $M_{c}$ and average energy $E_{c}$, that are independent of $\gamma$. In the thermodynamic limit, $L\to\infty$, the aforementioned equilibrium state becomes critical. For $\delta>0$, the system is expected to reach a nonequilibrium stationary state (NESS) the properties of which depends explicitly on $\gamma$ and $\delta$. Clearly, this NESS remains invariant under the transformation $\delta\to-\delta$. Our goal in this work is to characterize this NESS. In particular, we focus on the behaviour of magnetization $M$ per site and energy $E$ per bond, defined as,

[TABLE]

Typical time evolutions of $M(t)$ and $E(t)$ are shown in Fig. 1. Using numerical simulations, we measure the mean and variance of these observables in the NESS, and find that these quantities exhibit certain unusual features, including non-monotonic behaviour as a function of $\gamma$ for a finite sized system. We analytically characterise the behaviour in the limiting scenarios of small and large temperature switching rate $\gamma$. Before going into the details, we provide a brief summary of our main results.

• •

For large $\delta\lesssim T_{c}$, both average magnetization $\langle M\rangle$ and average energy $\langle E\rangle$ show non-monotonic behaviour with $\gamma$ — while $\langle M\rangle_{\gamma}$ first decreases and then increases, $\langle E(\gamma)\rangle$ first shows an increase, eventually decreasing as $\gamma$ is increased [see Figs. 2(a) and 4(a)]. For small $\delta$, the non-monotonicity disappears, where $\langle M\rangle_{\gamma}$ $\big(\langle E\rangle_{\gamma}\big)$ monotonically increases (decreases) as a function of $\gamma$. Moreover, the variance of magnetization $\sigma_{M}^{2}$ shows a highly non-trivial behaviour with $\gamma$, whereas the variance of energy $\sigma_{E}^{2}$ decreases monotonically with $\gamma$ [see Figs. 2(b) and 4(b)].

• •

In the small $\gamma$ regime the typical duration between two successive temperature switching events is much larger compared to the typical relaxation time of the Ising model at a given temperature. Using this time-scale separation, we obtain the behaviour of $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$ in the small $\gamma$ regime. This also explains the emergence of the non-monotonicity for large $\delta\lesssim T_{c}$ [see Fig. 6].

• •

Investigating the behaviour of $\langle M\rangle$ and $\langle E\rangle$ as a functions of $\delta$, we find that for small $\delta$,

[TABLE]

where $M_{c}$ and $E_{c}$ denote the equilibrium values of average magnetization and energy, respectively, at $T=T_{c}$. Using dynamical response theory, we investigate the $\gamma$ dependence of the second order susceptibilities $\chi_{2}^{M}$ and $\chi_{2}^{E}$ [see Fig. 7]. Interestingly, we find that both $\chi_{2}^{M}$ and $\chi_{2}^{E}$ change sign as $\gamma$ is increased beyond certain threshold values.

• •

We show that, in the $\gamma\to\infty$ limit, the spin-flip rate $w(\Delta E)$ can be expressed as that of an equilibrium Ising model with an effective temperature

[TABLE]

In fact, we find that for large $\gamma\gg\psi_{0}$, an effective equilibrium picture emerges, with a $\gamma$ and $\delta$ dependent effective temperature [see Fig. 8]. However, the inherent nonequilibrium nature of the process is apparent from a non zero energy current flowing through the system [see Fig. 10].

In the next section, we study the behavior of the NESS in detail.

3 Characterization of the NESS

To characterise the stationary state of this stochastic two-temperature nonequilibrium Ising model, we first focus on the behaviour of mean and variance of magnetisation and energy. For such a nonequilibrium setting it is generally hard to compute these quantities analytically. Hence we take recourse to Monte-Carlo simulations to explore the behaviour of these quantities for a finite system in the $\delta-\gamma$ plane. For the sake of completeness, we briefly describe the details of the simulation.

We initialize the system in a random configuration, where each spin $s_{i}$ is chosen to be $\pm 1$ with equal probability $1/2$. Additionally, we initialize $\sigma=\pm 1$ with equal probability $1/2$, which corresponds to a temperature $T_{\sigma}=T_{c}+\sigma\delta$. The temperature remains constant for a time $\tau$ which is drawn from an exponential distribution $P(\tau)=\gamma e^{-\gamma\tau}$, during which the spins evolve via Glauber dynamics—a randomly chosen spin flips with the rate $w(\Delta E)$ defined in (2). After time $\tau$ is elapsed, $\sigma\to-\sigma$ resulting in a temperature switch. The temperature now remains constant at $T_{-\sigma}$ for a time drawn from $P(\tau)$. This dynamics eventually leads to an NESS where we measure the relevant observables.

Let us first investigate the behaviour of magnetization and energy as functions of $\gamma$. Figure 2(a) shows a plot of $\langle M\rangle$ versus $\gamma$ for different values of $\delta$. Evidently, for $\delta=0$, the magnetization remains constant at $\langle M\rangle=M_{c}$, where $M_{c}$ is the average value of magnetization for the equilibrium $L\times L$ Ising model at temperature $T_{c}$. For $\delta>0$, average magnetization shows a non-trivial dependence on $\gamma$—for small values of $\gamma$, $\langle M\rangle<M_{c}$, while for large $\gamma$, $\langle M\rangle>M_{c}$. In fact, it appears that for some intermediate $\gamma=\gamma^{\star}_{M}$, $\langle M\rangle=M_{c}$ irrespective of $\delta$, which corresponds to the crossing point of the different curves in Fig. 2(a). Remarkably, for large values of $\delta$, i.e., as the system gets further away from equilibrium, the dependence of $\langle M\rangle$ on $\gamma$ becomes non-monotonic in the $\gamma<\gamma^{\star}_{M}$ regime. However, for $\gamma>\gamma^{\star}_{M}$, the magnetization monotonically increases for all values of $\delta$ indicating an overall increase in the spin alignment. Unusual signatures are also exhibited by the variance of the magnetization, $\sigma_{M}^{2}=\langle M^{2}\rangle-\langle M\rangle^{2}$, which shows a peak near $\gamma=\gamma^{\star}_{M}$ that becomes sharper with increasing $\delta$ [see Fig. 2(b)].

Thus from Fig. 2 it is apparent that the stochastic temperature switching leads to some intriguing features in the behaviour of the magnetization. These features are most prominent when the switching rate is smaller than the threshold value $\gamma^{\star}_{M}$. We also investigate how this threshold changes with the typical intrinsic spin-flip time-scale. Figure 3 shows a plot of $\gamma^{\star}_{M}$ as a function of $\psi_{0}$ which shows a linear increase. This indicates that the effect of stochastic evolution of the temperature is most striking when the typical time-scale of temperature switching is smaller than that of the spin flip.

We also explore the behaviour of mean and variance of the energy per bond as functions of $\gamma$. Figure 4(a) shows plots of $\langle E\rangle$ versus $\gamma$—similar to magnetization, the curves corresponding to different values of $\delta$ show a crossing a threshold value $\gamma^{\star}_{E}$. Note that, $\gamma^{\star}_{E}$ is different from $\gamma^{\star}_{M}$ although they are of the same order. Contrary to the magnetization, the average energy $\langle E\rangle>E_{c}$, the value of average energy at $\delta=0$, for $\gamma<\gamma^{\star}_{E}$, while $\langle E\rangle<E_{c}$ for $\gamma>\gamma^{\star}_{E}$. It is clear from Fig. 4(a) that $\langle E\rangle$ also shows a nonmonotonicity for large values of $\delta$ in the $\gamma<\gamma^{\star}_{E}$ regime. On the other hand, the variance of energy $\sigma^{2}_{E}=\langle E^{2}\rangle-\langle E\rangle^{2}$ decreases monotonically with $\gamma$ [see Fig. 4(b)].

We also investigate how average magnetization and energy change with $\delta$. As shown in Fig. 5 both $\langle M\rangle$ and $\langle E\rangle$ are symmetric functions of $\delta$. The change in the behaviour of the system across the threshold values of $\gamma$ is manifest in the change of curvature of the $\langle M\rangle$ and $\langle E\rangle$ curves. The average magnetization $\langle M\rangle$ is concave around $\delta=0$ when $\gamma<\gamma^{\star}_{M}$ and becomes convex for $\gamma>\gamma^{\star}_{M}$. In contrast, $\langle E\rangle$ is convex for $\gamma<\gamma^{\star}_{E}$ and turns concave for $\gamma>\gamma^{\star}_{E}$.

In the following, we first focus on the two limiting scenarios, $\gamma\ll\psi_{0}$ and $\delta\ll T_{c}$, to understand the physical origin of the non-trivial $\gamma$ and $\delta$ dependence of magnetization and energy. Finally we focus on the effective temperature picture emerging in the regime $\gamma\gg\psi_{0}$.

3.1 Small $\gamma$ Regime

To understand the behaviour of the NESS in the small $\gamma$ regime, it is useful to consider a renewal approach [27, 28]. Since the time evolution of the system between two consecutive temperature switches follows ordinary Ising dynamics at a constant temperature, we can write a last renewal equation,

[TABLE]

where, $P_{\gamma}(C,t|C_{0})$ denotes the probability that the system is in spin configuration $C$ at time $t$ starting from $C_{0}$. This equation can be understood by considering contributions from trajectories with no temperature switch, indicated by the first term and atleast one temperature switch, indicated by the second term. In the second term, $t^{\prime}$ denotes the time elapsed since the last temperature switch to the value $T_{\sigma}$, $\mathbb{P}_{T}(C,t^{\prime}|C^{\prime})$ denotes the probability that ordinary Ising model at constant temperature $T$ is in configuration $C$ at time $t^{\prime}$ starting from $C^{\prime}$. Moreover, $P_{\gamma}(C^{\prime},t-t^{\prime}|C_{0})$ denotes the probability that the system has configuration $C^{\prime}$ at the time of last temperature switch, $t-t^{\prime}$, starting from $C_{0}$. Here we have assumed that at $t=0$, the system can be at temperature $T_{\pm}$ with equal probability $1/2$, i.e., $\sigma_{0}=\pm 1$ with $\rho(\sigma_{0})=1/2$.

In general, it is hard to solve the integral (6) to obtain $P_{\gamma}(C,t)$ explicitly. However, in the small $\gamma$ limit, this equation can be used to obtain the stationary value of the average of any observable $\langle O\rangle$ as discussed below. For ordinary 2D Ising model with Glauber dynamics, the relaxation time $\tau_{r}\sim\psi_{0}^{-1}L^{2}$. For small values of $\gamma\ll\tau_{r}^{-1}$, the typical duration between consecutive temperature flips $\gamma^{-1}$ is much larger compared to the typical relaxation time of ordinary Ising model. Consequently, the system approaches close to equilibrium at temperatures $T_{+}$ and $T_{-}$ between successive temperature switches. In this regime, $\mathbb{P}_{T_{\sigma}}(C,t^{\prime}|C^{\prime})$ becomes independent of the initial spin configuration $C^{\prime}$ and (6) reduces to,

[TABLE]

The expectation values of any single-time observable $O(C)$, like magnetization and energy, in the stationary state, can thus be expressed as,

[TABLE]

The above expression allows us to compute the stationary average $\langle O\rangle_{\gamma}$ of any observable using its equilibrium relaxation behaviour $\langle O(t)\rangle_{\rm eq}^{\pm}$ for ordinary Ising model at temperature $T_{\pm}$. For any finite-sized system at temperature $T_{\sigma}$, the observables typically show an exponential relaxation to a equilibrium value $\langle O\rangle_{\rm eq}^{\sigma}$. Thus we can assume,

[TABLE]

where $\tau_{\sigma}$ denotes the relaxation time-scale at temperature $T_{\sigma}$, and $b_{\sigma}$ is an empirical constant. Substituting (9) in (8) and performing the integration over $t^{\prime}$, we get,

[TABLE]

Clearly, in the limit $\gamma\to 0$, we have $\langle O\rangle_{\gamma}=\Big[\langle O\rangle_{\rm eq}^{+}+\langle O\rangle_{\rm eq}^{-}\Big]/2$, as expected. Equation (10) allows us to determine the small $\gamma$ behaviour of $\langle M\rangle$ and $\langle E\rangle$ in the NESS, using the corresponding $\langle O\rangle_{\rm eq}^{\sigma}$, $b_{\sigma}$ and $\tau_{\sigma}$ values extracted from the equilibrium relaxation dynamics at $T_{\sigma}$. This can be done by looking at,

[TABLE]

for ordinary Ising model at $T_{\sigma}$ and extracting the parameters $\tau_{\sigma}$ and $b_{\sigma}$, which depend on the observable, by fitting it to the form given in (9). We estimate these parameters for a given $\delta$, for both energy and magnetization, from their equilibrium relaxation dynamics and compute their $\gamma$-dependence. Figure  6 compares these predicted behaviours with $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$ obtained from numerical simulations, which show an excellent match. This validates our prediction (10) and the assumptions leading to it.

This analysis also explains the non-monotonic behaviour of $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$. In the small $\gamma$ regime, at each temperature switching event $\sigma\to-\sigma$, the expectation value $\langle O(t)\rangle_{\gamma}$ starts from $\langle O\rangle_{\rm eq}^{\sigma}$, eventually reaching $\langle O\rangle_{\rm eq}^{-\sigma}$ before the next switch [see Fig. 1]. Thus, $b_{+}$ and $b_{-}$ must have opposite signs. For example, for magnetisation $b_{+}>0$ and $b_{-}<0$, while for energy $b_{+}<0$ and $b_{-}>0$. It is then easy to show that, for $\tau_{+}\neq\tau_{-}$, (10) has an extremum at some intermediate value of $\gamma$. Since, the relaxation in the ferromagnetic phase is much slower compared to that in the paramagntic phase, $\tau_{-}>\tau_{+}$, with their difference becoming larger as $\delta$ is increased. This difference in the relaxation time-scales leads to the non-monotonic variation of $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$ in the small $\gamma$ regime, as shown in Fig. 6.

3.2 Small $\delta$ regime

In this section, we study the behaviour of average energy and average magnetization for small values of $\delta$, i.e., when $\delta/T_{c}\ll 1$. The $\delta=0$ limit corresponds to the case of the equilibrium Ising model at temperature $T_{c}$. Therefore, the small-$\delta$ behaviour of the system can be explored by applying the dynamical response theory [29, 30, 31, 32] around the equilibrium limit, treating $\delta$ as the perturbation parameter. Since for $\delta=0$, the temperature flip dynamics is decoupled from the spin dynamics, the unperturbed system corresponds to two independent equilibrium processes, namely, the Ising dynamics at temperature $T_{c}$ and a dichotomous process, $\sigma\to-\sigma$ with rate $\gamma$. For non-zero values of $\delta$ the coupled dynamics reaches an NESS which remains invariant under the transformation $\delta\to-\delta$ and the expectation value of any observable $\langle O\rangle_{\delta}$ must be an even function of $\delta$ [see Fig. 5]. Therefore, for small $\delta$, we can write,

[TABLE]

where $\langle O\rangle_{\rm eq}$ denotes the equilibrium value of the observable and $\chi_{2}$ represents the second order response coefficient. In the following, we compute this coefficient using the non-linear dynamical response theory, where frenetic contributions play a significant role [30, 31]. For a generic observable $O(\omega)$, which depends on the trajectory $\omega$ of the system over a time-interval $[0,t]$, the second order response coefficient is given by,

[TABLE]

where $S^{\prime}(\omega)$ and $D^{\prime}(\omega)$ denote the first-order excesses in entropy and frenesy, generated due to the perturbation. The corresponding second-order contributions are indicated by $S^{\prime\prime}(\omega)$ and $D^{\prime\prime}(\omega)$. For Markov jump-processes, $S(\omega)$ and $D(\omega)$ along the trajectory $\omega$, have the general form [32],

[TABLE]

where $k(x_{t^{\prime}_{\ell}},x_{t^{\prime}_{\ell+1}})$ is the transition rate for the $\ell$-th jump from the configuration $x_{t^{\prime}_{\ell}}$ to $x_{t^{\prime}_{\ell+1}}$, $\lambda(x)=\sum_{y\neq x}k(x,y)$ denotes the escape rate from the state $x$ and $\sum_{\ell}$ denotes a sum over all transitions $\ell$ along the trajectory $\omega$.

Since we are interested in the response coefficient of single-time observable like $M(t)$ and $E(t)$, it is useful to consider a time-antisymmetric observable of the form $O(\omega,t)=O(t)-O(0)$ in (13). In this case, the the correlations $\langle(D^{\prime})^{2}O\rangle_{\rm eq}$, $\langle(S^{\prime})^{2}O\rangle_{\rm eq}$, and $\langle D^{\prime\prime}O\rangle_{\rm eq}$ vanish due to the time-reversal invariant nature of equilibrium dynamics. Consequently, to second order in $\delta$, we have,

[TABLE]

Clearly, the second order response coefficient is a sum of two terms—the first one, $\langle S^{\prime\prime}O(t)\rangle_{\rm eq}$ is purely entropic while the second one, $\langle S^{\prime}D^{\prime}O(t)\rangle_{\rm eq}$ involves the excess dynamical activity.

To proceed further, we note that for our two-temperature Ising model, the trajectory is characterised by $\omega=\{C(t^{\prime}),\sigma(t^{\prime});0<s<t\}$, where $C(t^{\prime})$ denotes the spin configuration and $\sigma(t^{\prime})$ indicates the instantaneous temperature $T(t^{\prime})=T_{c}+\sigma(t^{\prime})\delta$ at time $t^{\prime}$. The transition rates are given by,

[TABLE]

where $H$ is the Hamiltonian of the Ising system given in (1) and $\beta_{\sigma}=1/(T_{c}+\sigma\delta)$. Since we need to compute excesses of $S(\omega)$ and $D(\omega)$ with respect to $\delta$, it suffices to consider only the $\delta$ dependent terms,

[TABLE]

Here, $\Delta E^{\sigma}(t)$ in (19) is the total energy gained by the system from the reservoir at temperature $T_{\sigma}$ during the interval $[0,t]$. Moroever, $\Delta\varepsilon_{k}(t^{\prime})$ in (20) denotes the energy cost associated with flip of the $k$-th spin in the configuration $\{C(t^{\prime}),\sigma(t^{\prime})\}$.

The excess entropy can be computed by taking derivative of (19) with respect to $\delta$ and then taking $\delta=0$.

[TABLE]

where, $\mathscr{J}_{E}(t)$ denotes the time-integrated energy current flowing from the hot reservoir at temperature $T_{+}$ to the cold reservoir at temperature $T_{-}$. Moreover, $U(t)$ denotes the net energy gain of the system from both the reservoirs during the time interval $[0,t]$.

The excess frenesy can be similarly obtained from (20) and is given by,

[TABLE]

Substituting the above equations in (16), we explicitly obtain the second-order response coefficient, which can be conveniently expressed as,

[TABLE]

where,

[TABLE]

denotes the purely entropic contribution, and

[TABLE]

involves the contribution from frenesy, explicitly depending on the dynamics. In the limit $t\to\infty$, $\chi_{2}(t)$ reaches the stationary response coefficient.

It is noteworthy that the net energy gained $U(t)=H(t)-H(0)$, where $H(t)=-\sum_{\langle ij\rangle}s_{i}(t)s_{j}(t)$ denotes the total energy of the system at time $t$. Thus, the entropic component of the susceptibility (25) can be expressed as,

[TABLE]

In the limit $t\to\infty$, the correlation $\big\langle H(0)O(t)\big\rangle_{\rm eq}$ factorizes and $\chi_{2}^{(1)}$ can be expressed in the more familiar form of equilibrium response coefficient,

[TABLE]

where $V=-\delta^{2}H/T_{c}^{2}$ can be interpreted as the perturbation to the original Hamilonian (1). Evidently the entropic term $\chi_{2}^{(1)}$ is independent of $\gamma$ and the $\gamma$ dependence of the susceptibility $\chi_{2}$ originates solely from $\chi_{2}^{(2)}$.

Using Eqs. (24)-(26) we compute the stationary response coefficients $\chi_{2,M}$ and $\chi_{2,E}$ for magnetization and energy, respectively, from numerical simulations of the equilibrium Ising model at $T_{c}$. This is compared with the response measured directly from the perturbed system $\tilde{\chi}_{2,O}=\big(\langle O\rangle_{\delta}-\langle O\rangle_{\rm eq}\big)/\delta^{2}$ in Fig. 7, which shows an excellent match with the prediction (24). From Fig. 7, it is clear that the response coefficient vanishes at an intermediate value of $\gamma$, where,

[TABLE]

Thus, there exists a $\gamma^{*}_{O}$ where the expected value of the observable $\langle O\rangle$ in the NESS is equal to its equilibrium value. This explains the crossing of $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$ for different values of $\delta$ shown in Figs. 2(a) and 4(a).

It is now useful to consider what happens to the susceptibility in the limiting situations $\gamma\to 0$ and $\gamma\to\infty$. In the limit of vanishing switching rate $\gamma\to 0$, the temperature remains constant at its initial value $T_{\sigma_{0}}$, allowing us to substitute $\sigma(t^{\prime})=\sigma_{0}$ in (26), which leads to,

[TABLE]

Performing the sum over $\sigma_{0}$ and using (30) in (12), we get

[TABLE]

which infact is the small $\delta$ expansion of (10) in the $\gamma\to 0$ limit.

On the other hand, for large switching rate, $\gamma\gg\psi_{0}$, we find that the $\gamma$ dependent part of the response function $\chi_{2}^{(2)}\sim 1/\gamma$ [see Appendix A for details]. Hence, in the limit $\gamma\to\infty$, $\chi_{2}^{(2)}\to 0$ and the susceptibility becomes solely entropic in nature [see (28)]. This is suggestive of an equililibrium-like picture in the large $\gamma$ scenario, which we investigate in the following section.

4 Large $\gamma$ regime and effective temperature description

In this section, we focus on the behavior of the NESS in the large $\gamma$ regime, i.e., when $\gamma\gg\psi_{0}$. In this regime, the typical interval between two successive temperature switching events is much smaller than the typical time-scale of spin-flips. Hence, in the limit of infinite temperature switching rate, we can replace $\sigma$ by $\langle\sigma\rangle=0$ in (2) to write the effective spin-flip rate as,

[TABLE]

Clearly, $w_{\rm eff}(\Delta E)/w_{\rm eff}(-\Delta E)=\exp(-\beta_{\infty}\Delta E)$ and (32) corresponds to a detailed balanced Glauber dynamics at an effective temperature $T_{\infty}$. Thus, in the limiting scenario of $\gamma\to\infty$, we expect the spin configuration weight to have the Boltzmann form $P(C)\sim e^{-\beta_{\infty}H(C)}$. Physically, this corresponds to a scenario where the temprature switches happens so fast that the Ising spins are not able to distinguish between the two reservoirs.

Nonequilibrium stationary states can sometimes be partially characterised by an effective temperature picture [33]. For example, it has been shown that an Ising model in contact with two or more heat reservoirs can sometimes be described using an equilibrium-like picture with an effective temperature [34, 22]. An effective temperature picture also describes a generalised FDT for nonequilibrium critical dynamics of Ising model [35].

It is thus interesting to ask whether such an effective temperature picture also holds for the stochastic two-temperature Ising model for finite $\gamma$. In this case no such effective temperature can be directly obtained from the flip-rate. Instead, we extract an effective temperature $T_{\rm eff}(\gamma)$ by comparing $\langle M\rangle$ vs $\gamma$ curve with the equilibrium $\langle M\rangle$ vs T curve such that,

[TABLE]

for a given value of $\delta$. This process is illustrated in the inset of Fig. 8(a). Similarly, a $\gamma$ and $\delta$ dependent effective temperature can also be extracted from $\langle E\rangle_{\gamma}$. Figure. 8(a) shows the plots of these two effective temperatures as functions of $\gamma$ for a particular value of $\delta$—clearly for $\gamma\gg\psi_{0}$, $T_{\rm eff}(\gamma)$ becomes independent of the choice of the observable, indicating the existence of a unique effective temperature. Infact, the effective temperatures extracted from $\langle M^{2}\rangle$ and $\langle E^{2}\rangle$ also correspond to the same values of $T_{\rm eff}(\gamma)$, although these are not explicitly shown here. Note that in the limit $\gamma\to\infty$, $T_{\rm eff}(\gamma)\to T_{\infty}$. Figure 8(b) illustrates that $T_{\rm eff}(\gamma)-T_{\infty}$ decays algebraically as $\sim 1/\gamma$ for different values of $\delta$.

Such an effective temperature description is not expected to hold when the time-scales of temperature switch and spin-flip become comparable. This is also apparent from Fig. 8(a), which illustrates that the $T_{\rm eff}(\gamma)$ curves obtained from $\langle M\rangle_{\gamma}$ and $\langle E\rangle_{\gamma}$ no longer match, in the small $\gamma$ regime.

The existence of a unique effective temperature suggests that the weight of the spin configuration $C$ in the NESS admits a Boltzmann-like form $P(C)\sim\exp\big[-H(C)/T_{\text{eff}}\big]$ in the large $\gamma$ regime. Such a Boltzmann-like form for the staionary state weight however does not necessarily imply that the system is in equilibrium [36, 37, 38]. To understand whether the spin degrees indeed achieves equilibrium, we measure the effective spin-flip rates.

For fixed $\gamma$ and $\delta$, the effective spin-flip rate $w_{\rm eff}(\Delta E)$ can be estimated from the ratio of number of successful spin-flips to attempted spin-flips with energy change $\Delta E$. From these rates, an effective temperature can be extracted as,

[TABLE]

From numerical simulations, we extract $T_{\rm eff}$ for the allowed values of $\Delta E=\pm 4$ and $\pm 8$. These are compared with the $T_{\rm eff}(\gamma)$ obtained from the stationary values of the observables using (33) [see Fig. 9]. Remarkably, for $\gamma\gg\psi_{0}$, all the curves converge, indicating that the dynamics effectively does satisfy detailed balance. Thus it appears that the spin degrees, when considered independently of the bath degrees of freedom, exhibit an effective equilibrium-like behaviour in this large switching rate regime.

However, the full system consisting of the spin degrees of freedom and the two reservoirs is inherently nonequilibrium in nature. To illustrate this, we measure dynamical observables like the energy current flowing through the system, which should vanish in equilibrium. The average instantaneous current flowing to the system from the reservoir at temperature $T_{\sigma}$ in the stationary state is given by,

[TABLE]

where $\Delta E^{\sigma}$, introduced in Eqs. (21)-(22), is the total energy gained by system from the reservoir during the time interval $[0,t]$. Then the average instantaneous energy current flowing through the system from the hot reservoir to the cold reservoir is given by,

[TABLE]

Figure 10(a) shows the plot for $j_{E}$, obtained from numerical simulations, as a function of the switching rate $\gamma$ for different value of $\delta$. Evidently, $j_{E}$ increases with increasing $\gamma$ and eventually saturates. The average current also increases with increasing $\delta$ as shown in Fig. 10(b). Thus a finite energy current always flows from the hot reservoir to the cold reservoir, notwithstanding the effective temperature description in the large $\gamma$ regime. This is a signature of the inherent nonequilibrium nature of the stochcastic two-temperature Ising model dynamics.

5 Finite-size Effects

In the previous sections, we have analysed the behaviour of the NESS using various theoretical techniques and limiting scenarios that are independent of the system size. However, so far we have relied on a fixed system size to validate these results. It is therefore natural to examine how the observables considered depend on the system size. To this end, we investigate the system-size dependence of the average magnetization and the energy in the stationary state. Figure 11 shows plots of $\langle M\rangle$ and $\langle E\rangle$ as functions of $\gamma$ for different values of system size $L$. Clearly, the non-monotonic nature of both the average magnetization and energy remain robust across system size variation. In fact, the qualitative behaviour of both the observables remain same as the system size is varied, although the magnetization shows a pronounced quantitative dependence on $L$.

The inset in Figure 11(a) shows the behaviour of $\langle M\rangle$ as a function of the system size $L$ for different values of $\gamma$. Clearly, for small $\gamma\ll\psi_{0}$, the magnetization decreases with increasing $L$, while for large $\gamma$, it becomes largely independent of $L$. This can be understood as follows. As already discussed, in the $\gamma\to 0$ limit [see (10)], the magnetization is given by, $\langle M\rangle_{\gamma}=[\langle M\rangle_{\text{eq}}^{+}+\langle M\rangle_{\text{eq}}^{-}]/2$, where $\langle M\rangle_{\text{eq}}^{\sigma}$ denotes average magnetization of equilibrium Ising model at temperature $T_{\sigma}$. In the limit of thermodynamically large system size, $\langle M\rangle_{\text{eq}}^{+}\to 0$ while $\langle M\rangle_{\text{eq}}^{-}$ attains a finite value, leading to $\langle M\rangle_{\gamma}\simeq\langle M\rangle_{\text{eq}}^{-}/2$. On the other hand, as discussed in Sec 4, for $\gamma\to\infty$, the spin-configurations in the NESS has an effective equilibrium description at a temperature $T_{\text{eff}}<T_{c}$. This implies that for large $\gamma$, the system attains an ordered state, yielding a large value of the magnetization. The average energy shows an opposite trend—for small $\gamma$, $\langle E\rangle$ increases as $L$ is increased. Moreover, the overall system size dependence of the average energy is less pronounced, as can be seen from Fig. 11(b).

It should be noted that both $\langle M\rangle$ and $\langle E\rangle$ exhibit little dependence on the system size in the large-$\gamma$ regime, implying that the effective temperature at large $\gamma$ is essentially independent of the system size.

6 Conclusion

In this work, we have studied the two-dimensional Ising model subjected to a stochastic dichotomous temperature. The temperature alternates between $T_{c}\pm\delta$ with a rate $\gamma$ which leads the system to a nonequilibrium stationary state characterised by unusual behaviour of observables like magnetization and energy. For example, the average magnetization and average energy change non-monotonically as functions of flip-rate $\gamma$. We have characterised this NESS through two complementary approaches. In the small $\gamma$ regime, we use a renewal approach to explain the origin of the non-monotonic variation of the obseravbles. On the other hand, for small $\delta$, we employ the dynamical response theory to understand the $\gamma$ dependence of the observables.

We find that the NESS admits an effective temperature description for fast temperature switches, i.e., for large values of $\gamma$. We explicitly compute this effective temperature in the $\gamma\to\infty$ limit, which turns out to be below the critical temperature of the ordinary Ising model. Using extensive numerical simulation, we show that the effective temperature description remains valid for a finite range of $\gamma$ and breaks only in the small $\gamma$ regime. We further show that such an effective temperature description does not imply that the system reaches an equilibrium state—the inherent nonequilibrium nature of the dynamics leads to a non-zero energy current flowing through the system from the hot reservoir to the cold reservoir. We also show that the qualitative behaviour of the observables in the NESS remains robust across system size variation.

A natural open question is whether similar effects are observed in other systems that belong to the Ising universality class. In a more general sense, our study can be thought of as a special case of a driving protocol where a control parameter of a system is randomly switched between two distinct values according to a dichotomous process. In this context, it would be worthwhile to investigate systems with underlying nonequilibrium dynamics like the Katz-Lebowitz-Spohn model [36] or the driven lattice gas [39]. It would also be of interest to explore whether our results can be tested experimentally using thin-film set-ups similar to the ones used in Ref. [40, 41].

Acknowledgements

We acknowledge the National Supercomputing Mission (NSM) for providing computing resources of ‘PARAM RUDRA’ at S.N. Bose National Centre for Basic Sciences, which is implemented by C-DAC and supported by the Ministry of Electronics and Information Technology (MeitY) and Department of Science and Technology (DST), Government of India.

Funding information: RS acknowledges the CSIR grant no. 09/0575(11358)/2021-EMR-I. UB acknowledges the support from the Anusandhan National Research Foundation (ANRF), India, under a MATRICS grant [no. MTR/2023/000392].

Data Avaiability

The data that support the findings of this article are obtained from numerical simulations. These are available from the authors upon reasonable request.

Appendix A Response coefficient in the limiting case of infinite switching rate

In this appendix, we derive the expressions for second-order response coefficient $\chi_{2}(t)$ in the limiting scenario of infinite switching rate, i.e., for $\gamma\to\infty$. We start by considering the $\gamma$ dependent part of the response coefficient quoted in (26),

[TABLE]

In writing the above expression we have used,

[TABLE]

in (26), where $\mathscr{I}(t)$ denotes the instantaneous energy current flowing into the system. In the unperturbed dynamics, i.e., for $\delta=0$, the dichotomous process $\sigma\to-\sigma$ and the Ising spin-flip dynamics at temperature $T_{c}$ become statistically independent, allowing the averages over $\sigma$ and over the spin configurations to be factorised. Consequently, (A.1) now simplifies to

[TABLE]

For large $\gamma$, $\sigma$ effectively emulates a white noise with autocorrelation,

[TABLE]

Substituting (A.4) in (A.3) and subsequently performing the integral over $t^{\prime\prime}$, yields

[TABLE]

It is thus evident from the above expresion that in the limit $\gamma\to\infty$, $\chi_{2}^{(2)}$ vanishes like $1/\gamma$.