Ising model with stochastic resetting

We study the stationary properties of the Ising model that, while evolving towards its equilibrium state at temperature $T$ according to the Glauber dynamics, is stochastically reset to its fixed initial configuration with magnetisation $m_0$ at a constant rate $r$. Resetting breaks detailed balance and drives the system to a non-equilibrium stationary state where the magnetisation acquires a nontrivial distribution, leading to a rich phase diagram in the $(T,r)$ plane. We establish these results exactly in one-dimension and present scaling arguments supported by numerical simulations in two-dimensions. We show that resetting gives rise to a novel"pseudo-ferro"phase in the $(T,r)$ plane for $r>r^*(T)$ and $T>T_c$ where $r^*(T)$ is a crossover line separating the pseudo-ferro phase from a paramagnetic phase. This pseudo-ferro phase is characterised by a non-zero typical magnetisation and a vanishing gap near $m=0$ of the magnetisation distribution.

Stochastic resetting has seen enormous activities during the last few years [1], notably in the context of search processes which are ubiquitous in nature [2].The simple intuition behind resetting is as follows.If one is searching for a target via a stochastic process such as simple diffusion, it may take a long time due to trajectories that run off from the target.It is then advantageous to restart the search process with a certain resetting rate r from the same initial condition [3,4].The idea is that one may explore new pathways leading to the target, thereby reducing the search time.Recently, a large number of studies have shown that there is an optimal resetting rate r * that makes the search process most efficient [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].There is yet another interesting aspect of resetting dynamics in addition to optimising the search process.The resetting move breaks detailed balance and hence drives the system into a nontrivial non-equilibrium stationary state (NESS).Characterising such a NESS has recently become a problem of central interest in statistical physics [3,[22][23][24][25][26][27].
The phenomenon of resetting has found a large number of applications across disciplines.For example, in biology, the process of RNA polymerisation, which is responsible for the synthesis of RNA from a DNA template, is stochastically interrupted by backtracking [28,29].Similar notions are also found in the ecological context: for instance, animals such as rhesus monkeys [30,31], during the foraging period, are known to perform stochastic resetting to the previously visited sites and the effects of such memory induced resetting have been studied in several models [18,24,[32][33][34][35].In computer science, stochastic restarts can be used to reduce the running time of randomised search algorithms [36][37][38][39][40][41].Stochastic resetting has also been studied in the context of complex chemical processes as in the Michaelis-Menton reaction scheme [9,12], active run-and-tumble particles [42][43][44], biological traffic models [45], and also recently in quantum systems [46,47].(m) in the stationary state, in the presence of resetting, is shown schematically in three different regions in the (T, r) plane.For T > Tc, there is a crossover line r * (T ) (shown schematically by the dashed blue line) that separates the para phase for r < r * (T ) and the "pseudo-ferro" phase for r > r * (T ).In the para phase, the P stat r (m) has a divergent peak at m = mtyp = 0 (and no gap at m = 0), while in the pseudo-ferro phase, P stat r (m) vanishes at m = 0 with a peak at mtyp > 0 (and still no gap at m = 0).For T < Tc (ferro phase), a nonzero gap opens up at m = 0 in P stat r (m) and moreover the distribution peaks at mtyp > 0.
Most of the systems discussed above concern the effect of resetting on the statics and dynamics of a single particle (or equivalently for noninteracting systems).It is natural to ask how resetting affects the stationary state of a many-body interacting system.This question has been addressed in a number of interacting systems in one dimension, such as reaction-diffusion systems [48], fluctuating interfaces [49], exclusion and zero-range processes subjected to resetting [50][51][52].However, none of these systems, in equilibrium, exhibit a thermodynamic phase transition.It is then interesting to know how resetting affects a system that, in the absence of resetting, displays a phase transition in its equilibrium state.The simplest paradigmatic model that exhibits an equilibrium phase transition is the Ising model in d-dimensions, with d ≥ 2. It is then natural to ask the question: what kind of stationary state is reached if an Ising model, evolving under the natural Glauber dynamics, is subject to resetting (to its initial configuration) at a constant rate r?
It is useful to first recall the properties of the nearest neighbour Glauber Ising model in the absence of resetting (r = 0) [53].Let s i = ±1 denote the spin at site i of the Ising model.Starting from an initial configuration where the spins are independently chosen to have value ±1 with probability (1±m 0 )/2 (where m 0 ∈ [0, 1]), the individual spins flip according to the Glauber rate that satisfies detailed balance (see later for details).Let P ({s i }, t) denote the probability distribution of a spin configuration {s i } at time t.The most natural observable is the order parameter, i.e., the average magnetisation m(t) = 1 N i s i (t) where • • • denotes averaging over the probability measure P ({s i }, t) at time t.At late times, the system approaches the thermal equilibrium state, where the magnetisation m(t) approaches, irrespective of the initial value m 0 , the final value m eq > 0 for T < T c (ferromagnetic phase), and 0 for T ≥ T c (paramagnetic phase and at the critical point).Of course, What happens to the magnetisation m(t) at long times when a finite resetting rate r is switched on?This means that the system still evolves under the Glauber dynamics but, with a rate r, it now goes back to the initial configuration and the Glauber dynamics restarts.In this Letter we show that this nonzero resetting drives the system into a new non-equilibrium stationary state where the magnetisation has a nontrivial distribution and not a single value as in thermal equilibrium.This nontrivial stationary state arises from the fact that, even though after each resetting event the system starts the Glauber dynamics from the same initial configuration, the state of the system at time t is governed by the time of evolution since the last resetting event and this time is itself a random variable drawn from an exponential distribution with mean 1/r.Consequently, the measured thermally averaged magnetisation m(t) fluctuates from one resetting history to another.In particular, in the stationary state it acquires a nontrivial distribution.This is thus markedly different from the equilibrium case (r = 0) where the magnetisation distribution is trivially a delta function centred either at m = 0 (for T > T c ) or at m = m eq > 0 (for T < T c ).Thus the knowledge of the full distribution of the magnetisation is necessary to characterise the steady state of the system in the presence of resetting.
It is useful to first summarise our main results.We show that a nonzero resetting leads to a rich phase diagram in the (T, r) plane as displayed in Fig. 1 with the emergence of a new phase which we call "pseudo-ferro" phase.At all temperatures, the stationary magnetisa-tion distribution has now a finite support.For T > T c , there is a new crossover line r * (T ) that separates the paramagnetic phase r < r * (T ) from a pseudo-ferro phase for r > r * (T ).In the para phase, the stationary magnetisation distribution P stat r (m) diverges as m → 0, as P stat r (m) ∼ m ζ with ζ = r/r * (T ) − 1 < 0 and hence the typical magnetisation m typ = 0 (m typ denotes the value of m at which P stat r (m) reaches its maximum).In addition, there is no gap at m = 0, i.e. g = 0.In contrast, in the "pseudo-ferro" phase (T > T c and r > r * (T )), while the gap g still remains zero, P stat r (m) now vanishes as m → 0 as P stat r (m) ∼ m ζ with ζ = r/r * (T )−1 > 0. Consequently, the maximum of P stat r (m) occurs at a nonzero value m typ > 0 (see Fig. 1).For T < T c (ferro phase), the distribution P stat r (m) has a finite support [m eq , m 0 ] (for m eq < m 0 ) or over [m 0 , m eq ] (if m eq > m 0 ).Thus, in this phase, there is a finite nonzero gap g = min(m eq , m 0 ).In addition, m typ > 0 in the ferro phase.Exactly at T = T c , the distribution P stat r (m) vanishes extremely rapidly, P stat r (m) ∼ e −Am −κ as m → 0. We show that the exponent κ is related to the equilibrium critical exponents via the relation κ = νz/β, where ν and β are respectively the correlation length and the order parameter critical exponents, while z is the dynamical critical exponent associated to the Ising Glauber dynamics at T = T c .We establish these results from an exact solution in d = 1 and, for d = 2, we provide scaling arguments supported by numerical simulations.
We start with an Ising model with ferromagnetic nearest-neighbour interactions H = −J i,j s i s j on a d-dimensional lattice with N sites and periodic boundary conditions.Starting from an initial condition where the spins are independently ±1 with probability (1 ± m 0 )/2, the Glauber dynamics, in the absence of resetting, consists in flipping a single spin with rate [53] where β = 1/(k B T ) is the inverse temperature and ∆E = 2Js i j∈n.n.s j is the change of energy in flipping the i-th spin.In d = 1, this rate simplifies to w( ) where γ = tanh (2βJ).This property makes the 1d Glauber dynamics exactly solvable as the evolution equation for the n-point correlation functions only involve n-point functions, i.e. it satisfies a closure property [53].This closure property however does not hold for d > 1.Note that under the Glauber dynamics, the magnetisation m(t) = (1/N ) i s i (t) evolves deterministically with time t.Now imagine that we switch on the resetting mechanism, whereby the system goes back randomly in time to the initial configuration with a nonzero rate r.This means that, between two successive resetting events, the system evolves by the standard Glauber dynamics mentioned above (1).If we now observe the system at a fixed time t, what matters is the time τ elapsed since the last resetting before t.This is because the system has evolved by the pure Glauber dynamics during the interval [t−τ, t].But since the resettings happen stochastically, the time τ itself is a random variable.As a result, any observable, such as the average magnetisation, measured at time t also becomes a random variable.One can express the distribution P r (m, t) of the average magnetisation m in the presence of resetting with rate r by the simple renewal equation P r (m, t) = r t 0 dτ e −rτ P 0 (m, τ ) + e −rt P 0 (m, t) , (2) where P 0 (m, τ ) = δ(m − m(τ )) denotes the magnetisation distribution in the absence of resetting (r = 0) since it evolves deterministically as m(τ ).The second term in ( 2) is easy to explain: it corresponds to having no resetting up to time t and the system evolves by the standard Glauber dynamics during [0, t] and ends up with a magnetisation m at time t.The first term in (2) corresponds to the event that there is a resetting event at time t − τ which happens with probability r dτ , followed by no-resetting in the interval [t − τ, t] which occurs with probability e −rτ .During this interval of length τ , the system evolves via the standard Glauber dynamics, hence at time t, the magnetisation is just m(τ ).In the large time limit, the second term in (2) drops out and hence the stationary magnetisation distribution is given by Given this simple renewal property (3), we need just to know the deterministic Glauber evolution m(τ ) for all τ in the absence of resetting to determine P stat r (m).Resetting in the 1d Ising model.We start with the exactly solvable case on a 1d-lattice with N sites and periodic boundary conditions in the absence of resetting.The one-point average s i (t) evolves via [53] with initial condition s i (0) = m 0 for all i.Thus the average magnetisation m(t) = (1/N ) i s i (t) evolves via dm(t)/dt = −(1 − γ)m(t) whose solution is trivially Substituting this solution (5) in Eq. ( 3), one obtains exactly where In this case T c = 0 and we have only the part of the phase diagram in Fig. 1 with T ≥ T c .The result in Eq. ( 6) clearly shows that, near m = 0, P stat r (m) ∼ m ζ where the exponent ζ = r/r * (T ) − 1 varies continuously with temperature.Thus P stat r (m) either diverges (for r < r * (T )) or vanishes (for r > r * (T )) as m → 0. In the former case, m typ = 0 -this is the para phase.In contrast, for r > r * (T ), m typ > 0: this is the new "pseudo-ferro" phase induced by resetting.We have also done numerical simulations in d = 1 to verify our analytical prediction in Eq. ( 6) and found excellent agreement (see Fig. 1 in the Supplementary Material [59]).Resetting in the 2d Ising model.In 2d, the Ising model at equilibrium has a finite T c ≈ 2.269 with the choice J = k B = 1.Unlike in d = 1, the Glauber dynamics is not exactly solvable in d = 2.However, using the well established phenomenological behaviour of m(t), in particular at late times, in the renewal equation ( 3), we can make some predictions for the stationary magnetisation distribution P stat r (m) in various parts of the phase diagram in the (T, r) plane in Fig. 1.We then verify these predictions with numerical simulations and find a very good agreement.Below we consider the three cases T > T c , T < T c and T = T c separately.
T > T c .We start with the paramagnetic phase T > T c .In this case, the average magnetisation, for the pure Glauber dynamics is expected to decay at late times as m(t) ∼ a 1 e −λ1 t where the amplitude a 1 and the leading decay rate λ 1 both depend on temperature [55] and can be estimated very precisely from Monte Carlo simulations.This pure exponential decay of m(t) holds only when t 1/∆λ where ∆λ is the first gap in the relaxation spectrum.Substituting this functional form in Eq.
(3) we get, for r ∆λ This is a good approximation at high temperature T T c where ∆λ is large.In this case a 1 ≈ m 0 and thus we recover qualitatively a similar distribution for P stat r (m) (6) as in the d = 1 case.In the case of 2d, for T T c , we then have r * (T ) = λ 1 .Thus, as in the d = 1 case, we have a crossover from the usual para phase for r < r * (T ) to the "pseudo-ferro" phase for r > r * (T ) across the crossover line r * (T ) in the (T, r) plane, for T T c .Our numerical simulations are completely consistent with this scenario.In Fig. 2 a), we plot the cumulative stationary distribution F stat r (m) = m 0 P stat r (m )dm vs m for three different resetting rates: one in the para phase (the top curve), one in the "pseudo-ferro" phase (bottom curve) and finally at the crossover line r = r * (T ) (middle curve).In the last case, the probability distribution function (PDF) of the magnetisation P stat r * (m) is uniform [from Eq. ( 8)] and hence the cumulative distribution function (CDF) F stat r * (m) increases linearly with m.T < T c .In the ferro phase and in the absence of resetting, the magnetisation density of the 2d Ising model reaches a nonzero equilibrium value m eq with a stretched exponential decay at late times [55] where the parameters a, b and 0 < c < 1, not known analytically, need to be determined from simulations.In (9), the + and − signs are used in the case m 0 > m eq and m 0 < m eq respectively.When a constant resetting rate r is introduced in the system, the stationary PDF P stat r (m) is obtained from the general formula in Eq. ( 3), where for m(τ ) we now use Eq. ( 9).The resulting P stat r (m) is non-trivial and its detailed form is discussed in [59].For instance, a plot of the associated CDF is given in Fig. 2 c) for the case m 0 > m eq and compared to simulations, showing an excellent agreement.In this case, the PDF is supported over the finite interval [m eq , m eq + a].It has non-trivial asymptotic behaviours at the edges.For example, near the lower edge, where m → m + eq , P stat r (m) vanishes faster than a power law as Exactly at the critical point, the magnetisation m(t), without resetting, has a non-monotonic decay with time [56][57][58], which is not known analytically, except at short and long times.Hence it is difficult to evaluate P stat r (m) exactly from Eq. ( 3) for all m.However, when r and m are both small, one can use in Eq. ( 3) the late time form of m(t) ≈ b c t −φ where the exponent φ = β/(νz) is related to the standard critical exponents defined earlier.We then find that there is a scaling regime as m → 0, r → 0 but with m r −φ fixed where the distribution P stat r (m) takes the scaling form where the scaling function G(y) = A y −1−1/φ e −(bc/y) 1/φ , with A = b (m) has still a finite support m ∈ [0, Γ] where the upper cut-off Γ depends on system parameters.While there is no strict gap at m = 0 (as in the ferro phase), P stat r (m) vanishes extremely rapidly as m → 0. Thus, even in this resetting induced NESS, there is a remnant signature of the equilibrium critical point T c that is manifest in this essential singularity near m = 0.
To summarise, in this Letter, we have addressed a general question: how does resetting affect a many-body interacting system that, in equilibrium, exhibits a thermodynamic phase transition at T = T c ?A natural candidate to study this question is the paradigmatic Ising model evolving under the Glauber dynamics.We have shown, both analytically and numerically, that resetting (with a constant rate r) leads to a nontrivial phase digram of the Glauber-Ising model in the (T, r) plane (see Fig. 1).In particular, we have shown that resetting leads to the emergence of a new pseudo-ferro phase for T > T c and r > r * (T ) where the system has a non-zero typical magnetisation in the stationary state and yet there is no gap in the magnetisation distribution near m = 0.
The qualitative features of the phase diagram, established here for the d = 1 and d = 2 Ising model with Glauber dynamics, are also expected to hold in higher dimensions, as well as for other single spin-flip dynamics, such as the Metropolis dynamics.Finally, going beyond the magnetisation distribution, it would be interesting to investigate the structure of the distribution of the twopoint correlation functions in this resetting induced nonequilibrium stationary state.

I. RESETTING IN THE 1d ISING MODEL: MONTE CARLO SIMULATIONS
For the 1d Ising model, T c = 0 and hence, there is no ferromagnetic phase for T > 0. Hence, in the phase diagram of Fig. 1 of the main text, we have only the phase T ≥ T c = 0.The stationary probability distribution of the average magnetisation in the 1d Ising model in the presence of resetting is given by Eq. ( 6) in the main text and reads where We have performed Monte-Carlo simulations on a one-dimensional lattice of length N = 10000 at temperature T = 3.5 (setting J = k B = 1) to verify our analytical results.Simulation data are collected by averaging over 3000 independent realisations.The fixed initial configuration to which the system is reset at a constant rate r has magnetisation m 0 = 0.992.We perform the Glauber dynamics of the 1d Ising model at temperature T = 3.5, for which r * (T = 3.5) = 1 − tanh(2/3.5)= 0.4836.Hence, for r < r * (T = 3.5), P stat r (m) diverges as m → 0. This is clearly shown in Fig. 1 a), where, in the presence of a resetting rate r = 0.2, the typical value of the average magnetisation is m typ = 0: this regime is the para phase.In contrast, for r > r * (T = 3.5), a new "pseudo-ferro" phase emerges, where m typ > 0 and P stat r (0) = 0, as shown in Fig. 1 b).Exactly on the crossover line r = r * (T ), Eq. ( 1) predicts a uniform distribution for P stat r (m) in [0, m 0 ]: this is numerically verified in Fig. 1 c).It is also interesting to know how the first moment of P stat r (m) depends on the resetting rate r.This can be computed easily as where r * (T ) = 1 − γ.We have measured this quantity m(r) for different values of r from our numerical simulations.They are plotted in blue in Fig. 2, showing an excellent agreement with Eq. (3).However, as discussed in the main text, the average magnetisation is not enough to characterise the stationary state under resetting and one needs the full distribution P stat r (m).For the Ising model in d ≥ 2, the critical temperature T c is finite.Hence, in the phase diagram of Fig. 1 of the main text, we will have all the three phases.Unlike in 1d, the Glauber dynamics is not exactly solvable in d ≥ 2, even in the absence of resetting.In this section, we provide some general scaling arguments to predict scaling forms for P stat r (m) for T > T c , T < T c and T = T c , using the renewal equation (3) in the main text which reads A. The high temperature phase T > Tc In the high temperature phase, the average magnetisation, quite generally, decays with time as m(t) = i a i e −λit where λ 1 < λ 2 < • • • denote the eigenvalues associated to the relaxation spectrum.At late times, this is dominated by the lowest eigenvalue λ 1 , i.e. m(t) ≈ a 1 e −λ1 t .This pure exponential decay is thus valid when t 1/∆λ where ∆λ = λ 2 − λ 1 is the first gap in the relaxation spectrum.Substituting this pure exponential form in Eq. ( 4) we get, for r ∆λ This is Eq. ( 8) of the main text and is valid for T T c , such that the first gap ∆λ is large.

B. The low temperature phase T < Tc
In the low temperature phase, the average magnetisation m approaches the equilibrium magnetisation m eq > 0 as a stretched exponential at late times [2], i.e., m(t) ≈ m eq ± a e −bt c , where the + and − signs are used in the cases m 0 > m eq and m 0 < m eq respectively.
Inserting this expression for m(t) in Eq. ( 4) of the main text, we obtain for r ∆λ where Its support is (m eq , m eq + a] if m 0 > m eq or [m eq − a, m eq ) if m 0 < m eq .In the second equality we have used x = a e −bτ c .

C. At the critical point T = Tc
At T = T c , the late time evolution of m(t) in the absence of resetting follows a power law decay to 0 since the spectrum is gapless.In this case, at late times, m(t) ≈ b c t −φ where the exponent is φ = β/(νz) [3].More precisely, this power law decay holds for t < L z , where L is the linear size of the system.This condition is fulfilled in our simulations since 1 rmin L z , where r min is the smallest resetting rate used in our simulations.Using Eq. ( 4) of the main text, we obtain for r 1 where in the second equality we made the change of variable x = b c t −φ .In Eq. ( 7), if we set the scaling variable y = mr −φ , then P stat r (m) takes the scaling form given in Eq. ( 10) of the main text.

III. NUMERICAL ESTIMATE OF THE PARAMETERS IN THE TIME EVOLUTION m(t) OF THE 2d ISING MODEL
In order to use Eq. ( 3) of the main text and compute P stat r (m), we need to estimate the parameters in the time evolution m(t).To this end, we have done Monte Carlo simulations of the Glauber dynamics in the three temperature regimes on a 256 × 256 square lattice with the choice J = k B = 1 and with m 0 = 0.9905 as the magnetisation of the starting configuration at time t = 0. Fig. 3 shows the fit of m(t) in the paramagnetic phase (Fig. 3 a)), ferromagnetic phase (Fig. 3 b)) and at the critical temperature (Fig. 3 c 8 and ν = 1, one gets z = 1/(8φ) ≈ 2.17, which is in excellent agreement with previous large-scale simulations [2].
FIG. 3. Fit of m(t) in the three temperature regimes.For T < Tc and T = Tc, it is necessary to consider much more time steps to obtain an accurate fit, since the decay to the equilibrium magnetisation is slower than the one for T > Tc.
IV. FIRST MOMENT OF P stat r

IN THE 2d ISING MODEL
A quantity of interest that clearly depends on the resetting rate r is the first moment of P stat r (m) in the 2d Ising model.Similarly to the approach followed for the 1d case, we can easily compute it for T > T c , T < T c and T = T c as where m(t) is the time evolution of the average magnetisation in the Glauber dynamics which has been estimated in the previous section.
For T > T c , using m(t) ≈ a 1 e −λ1t at late times, the average magnetisation in the stationary state is given, for small r (so that we can use the pure exponential late time decay of m(t)), by similarly to the one dimensional case Eq. ( 3).
For T < T c , using m(t) ≈ m eq ± a e −bt c , m(r) takes the form, for small r, where the last integral can be estimated numerically very accurately.At the critical temperature T c , using m(t) ≈ b c t −φ , the average magnetisation in the stationary state reads, for small r, where Γ(z) is the gamma function.
We have performed Monte Carlo simulations with different values of r to verify the dependence of m(r) on r.The agreement is excellent in all the three phases (Fig. 4 a), Fig. 4 b) and Fig. 4 c)), especially for small values of r.This is consistent with the fact that we have taken the late time evolution of m(t) in Eq. ( 8), which holds only for small r.  9), Eq. ( 10) and Eq. ( 11)), while the blue dots come from Monte Carlo simulations.

V. GENERAL EXPRESSION OF THE CUMULATIVE DISTRIBUTION
In the main text, using renewal theory, we have obtained the expression of the non-equilibrium stationary PDF of the average magnetisation in presence of a resetting rate r: where m(τ ) denotes the deterministic Glauber evolution of the average magnetisation.It is easy to compute exactly the CDF that we plot in Fig. 2 of the main text.We can write P stat r (m) as a slightly more general expression where of course in our case f (τ ) = r e −rτ and g(τ ) = m(τ ).Then, making the change of variable z = g(τ ), we can write the previous integral as Since f (t) = −r f (t), the corresponding CDF is given by where the sign is + and const = 0 if g is a decreasing function of time, while the sign is − and const = 1 if g is an increasing function of time.Since in our case f (t) = r e −rt , Eq. ( 16) reduces to This is the CDF that we plot in the main text, with the use of different functions g in the three temperature regimes.The derivation of Eq. ( 17) is based on two properties of the functions f and g.The first one is that f and its time derivative f are proportional, since f is an exponential function.The second property is that the function g, which describes the (deterministic) time evolution of the average magnetisation, is a monotonic function and, therefore, is invertible.In particular, this property of the function g is valid in the time window we considered for all the three different temperature regimes (paramagnetic and ferromagnetic phases and at the critical temperature).

VI. NUMERICAL ESTIMATE OF r * (T ) IN THE PARAMAGNETIC PHASE OF THE 2d ISING MODEL
The exact expression of the crossover line r * (T ) in the 1d Ising model is given by Eq. ( 7) of the main text as r * (T ) = 1 − tanh( 2J k B T ).It is then interesting to know if a similar curve can be obtained, at least numerically, also in the two dimensional case.For several temperatures T > T c , we perform numerical simulations of the 2d Glauber dynamics in presence of several resetting rates r and detect the value r * (T ) that makes the resulting PDF P stat r * (m) more similar to the uniform distribution in [0, m 0 ].We thus expect to obtain a curve r * (T ) that divides the para phase from the new "pseudo-ferro" phase as depicted in Fig. 1 of the main text.
If the time evolution of m(t) was given by a simple exponential decay to 0, then we would find that r * (T ) = λ 1 (T ), as in the one dimensional case.However, this is not the case, since this time dependence holds only at late times.Indeed, the time evolution of the average magnetisation in the paramagnetic phase is given by where λ 1 < λ 2 < • • • < λ N constitute the whole relaxation spectrum, a i are coefficients and N = 2 L 2 is the total number of spin configurations in a L × L lattice.Eq. ( 18) is valid for any time t and is therefore a generalisation of the simple exponential decay.In the long-time limit the two time evolutions coincide, because only the minimum value among λ i 's (i.e.λ 1 ) survives as t → ∞.
In principle, N should be equal to the total number of spin configurations.Of course, to make a reasonable fit, this number is too big: therefore we choose to take N = 3, that is we approximate m(t) to be the sum of the three exponentials with the smallest decay rates.As a consequence, from the fit of m(t) we estimate six parameters, which are a i and λ i for i = 1, 2, 3.
It is therefore interesting to compare the curve r * (T ) obtained as illustrated above with the value λ 1 (T ) obtained by fitting this m(t) to numerical data, for different temperatures T > T c .The inset of Fig. 5 shows the values of r * and λ 1 for different temperatures above T c .Note that r * (T ) > λ 1 (T ) ∀T > T c .In particular, if we imagine to draw a line connecting all the blue dots, we obtain the red dashed line in the phase diagram in Fig. 1 of the main text.In the main plot of Fig. 5 the ratio r * /λ 1 is plotted for different temperatures above T c , using the same data points of the inset.We see that this ratio converges to 1 as the temperature T increases, proving that, when T T c , the approximation r * = λ 1 (which is the exact relation found in the 1d case) is correct.

FIG. 1 .
FIG. 1. (Color online) Phase diagram in the (T, r) plane.The magnetisation distribution P stat r

FIG. 2 .
FIG. 2. (Color online) Plot of the CDF F stat r (m) = m 0 P stat r (m ) dm vs m, for different values of the resetting rate r, in a) the para phase T > Tc, b) at the critical point T = Tc ≈ 2.269 (for J = kB = 1) and c) in the ferro phase T < Tc. a): the solid curves show the simulation results for T = 3.5, m0 = 0.9905 and for three different values of r, respectively for r < r * (top blue curve), r = r * ≈ 0.117 (middle magenta curve) and r > r * (bottom green curve).The red dashed lines correspond to our theoretical prediction obtained by integrating Eq. (8).The agreement with the theoretical predictions gets better for smaller r as explained in the text.b): scaling collapse at T = Tc of the CDF F stat r (m) ≈ H(m r −β/(νz) ), compared with the theoretical function H(y) = e −(bc/y) νz/β (red solid line) for three different small values of r.Here we used β = 1/8, ν = 1 and z ≈ 2.17 together with the estimated parameter bc = 0.9576.c): CDF for T = 2.24 < Tc with meq = 0.70732 and r = 0.00459 (blue solid line), compared with the theoretical prediction (see the text and Eq.(23) in the Supp.Mat.[59]) shown by the dashed red line.All the simulations were performed on a 256 × 256 square lattice.One observes finite size effects at the edges of the support.

c
/φ, vanishes extremely rapidly as y → 0. Consequently the CDF F stat r dm ≈ H(m r −φ ) where H (y) = G(y) = e −(bc/y) 1/φ .This scaling behaviour is verified numerically in Fig. 2 b) where F stat r (m) shows a beautiful scaling collapse for three different values of r.The full distribution P stat r

Fig. 1
presents the plot of the PDF P stat r (m) for three different values of r, showing the excellent agreement between theory (red line) and Monte-Carlo simulations (plotted in blue).

FIG. 1 .
FIG. 1. Plot of the PDF P stat r (m) for three different values of r: (a) para phase with r = 0.2 < r * (T = 3.5) = 0.4836 . .., (b) pseudo-ferro phase with r = 3 > r * (T = 3.5) and (c) exactly on the crossover line r = r * (T = 3.5).The red solid lines correspond to the theoretical result in Eq. (1), while the blue curves correspond to numerical simulations.The almost blue vertical lines at the edges of the distribution are due to finite size effects.
II. SCALING PREDICTIONS FOR P stat r (m) IN THE THREE TEMPERATURE REGIMES FOR THE ISING MODEL WITH RESETTING IN d ≥ 2

FIG. 2 .
FIG. 2. Plot of the first moment m(r) as a function of r.Numerical data and Eq.(3) are plotted in blue and red respectively.
)).The blue dots (simulation data) are obtained by averaging over 4770, 490 and 70 independent realisations at temperatures T = 3.5, T = 2.24 and T = T c = 2.269 respectively.The red dashed lines are obtained by fitting the well established phenomenological behaviour of m(t) in the three temperature regimes to the numerical data.As mentioned above, for T > T c , the average magnetisation decays at late times as m(t) ≈ a 1 e −λ1t .The fit shown in Fig. 3 a) provides the values a 1 = 0.889 and λ 1 = 0.117 for T = 3.5.For the low temperature phase, we set T = 2.24 < T c = 2.269 and in Fig. 3 b) we show the fit of the functional form m(t) ≈ m eq + a e −bt c with estimated parameter values a = 0.283, b = 0.37 and c = 0.316.At the critical temperature T c = 2.269 the fit, m(t) = b c t −φ , shown in Fig. 3 c) provides an estimation of the parameters b c = 0.9576 and φ = 0.0576.Using the relation φ = β/(νz), and the known exact values of β = 1

FIG. 5 .
FIG. 5.The ratio r * λ 1 and, in the inset, the values of λ1 and r * are estimated with Monte Carlo simulations at different temperatures above Tc.