Designing non-equilibrium states of quantum matter through stochastic resetting

We consider closed quantum many-body systems subject to stochastic resetting. This means that their unitary time evolution is interrupted by resets at randomly selected times. When a reset takes place the system is reinitialized to a state chosen from a set of reset states conditionally on the outcome of a measurement taken immediately before resetting. We construct analytically the resulting non-equilibrium stationary state, thereby establishing an explicit connection between quantum quenches in closed systems and the emergent open system dynamics induced by stochastic resetting. We discuss as an application the paradigmatic transverse-field quantum Ising chain. We show that signatures of its ground-state quantum phase transition are visible in the steady state of the reset dynamics as a sharp crossover. Our findings show that a controlled stochastic resetting dynamics allows to design non-equilibrium stationary states of quantum many-body systems, where uncontrolled dissipation and heating can be prevented. These states can thus be created on demand and exploited, e.g., as a resource for quantum enhanced sensing on quantum simulator platforms.

Introduction. Deterministic or stochastic dynamics interspersed with stochastic resetting leads to the emergence of a non-equilibrium steady state (NESS). This concept has been put forward in the classical realm in Refs. [1][2][3][4][5], where it has been established that a particle undergoing Brownian motion reaches a NESS when its diffusive time evolution is reset at constant rate γ (Poissonian reset) to its initial position. Using the renewal equation approach [1,2,[6][7][8][9][10][11][12][13], it has been further shown that stochastic resetting can improve the efficiency of search processes in terms of their mean first passage time [1,2,[14][15][16][17][18][19][20][21][22][23][24]. In contrast to the case of classical systems, much less is known about the effect of stochastic resetting on quantum systems, either closed or open. In the latter case, the Markovian Lindblad equation [25,26] describing the evolution in the presence of Poissonian resetting has been analyzed in Refs. [27][28][29][30][31][32]. Closed quantum systems which are reinitialized to their initial state after Poissonian resetting have been, instead, analyzed without reference to the Lindblad equation in Ref. [33]. Using the renewal equation approach, it was shown in the latter work that the reached NESS differs from the diagonal ensemble for any finite, non-zero, value of γ. For γ → ∞, instead, one recovers the quantum Zeno effect [86,87], where the system does not evolve in time due to frequent reset projections. Resetting in closed quantum systems has been also studied in Refs. [34][35][36][37][38][39] for first detection problems. These recent works suggest that there is a direct link between the dynamics of closed quantum systems subject to stochastic resetting and an effective open system dynamics. However, the precise nature of this relation, and of possible connections between the NESS emerging from the resetting, unitary dynamics, and non-equilibrium signatures of quantum phase transitions [40][41][42] is not fully understood.
In this work we establish such a connection by developing a general theory based on semi-Markov processes [30,43]. We prove analytically that by averaging the microscopic unitary time evolution over the reset distribution an effective non-Markovian open dynamics, governed by a generalized Lindblad equation, emerges. The ensuing NESS can be constructed analytically and its structure establishes a transparent connection to quantum quenches in closed systems, which can display signatures of equilibrium quantum phase transitions [44][45][46][47][48][49][50][51]. Our results for the NESS recover the quantum Zeno effect in the limit of infinite resetting rate. We illustrate our ideas using, as an example, the quantum Ising chain in a transverse field (TFIC). The ground-state quantum phase transition of the model is signalled in the NESS by a crossover, whose sharpness is controlled by the distribution of the reset time. Our analysis establishes stochastic resetting as a tool to generate complex stationary many-body quantum states without detrimental processes, such as heating, that typically occur when coherent dynamics is competing with dissipative processes.
Reset protocol The results we are going to develop in the following hold for a general many-body quantum system made of N particles or spin degrees of freedom, but for concreteness we will base our discussion mostly on the paradigmatic TFIC with Hamiltonian Here σ x,y,z n are the Pauli matrices acting on site n, J > 0 is the ferromagnetic coupling constant and h is the strength of a transverse magnetic field. This model exhibits a quantum phase transition in its ground state at the critical field strength h = h c = 1, which can be characterized by using the magnetization density m = N n=1 σ x n /N as order parameter. We envisage a "conditional reset" protocol, which is sketched in Fig. 1(a): Starting from a given initial state, in the Figure the state |ψ 0 = |↑↑ . . . ↑ [52]), the system is evolved unitarily under its many-body Hamiltonian H, e.g., Eq. (1) for the TFIC. After a randomly selected time, an observable is measured, in the Figure we measure the magnetization, and the many-body state is reset to a state chosen within a set of "reset states" conditionally on the outcome of the measurement. In our example the reset states are |1 = |↑↑ . . . ↑ and |2 = |↓↓ . . . ↓ (with |↑ and |↓ referring to the longitudinal-x direction in Eq. (1)). Our motivation to investigate the two aforementioned reset states is twofold: first, they are simple product states easy to prepare experimentally; second, they allow to establish a connection between the NESS of the reset dynamics and quantum quenches in the Ising chain. The two reset states can be indeed identified with the two degenerate ground states |GS(h = 0) ± of H in Eq. (1) for h = 0, namely, |1 = |GS(h = 0) + and |2 = |GS(h = 0) − . The unitary dynamics from the reset state |j (j = 1, 2), |ψ j (τ ) = exp(−iHτ ) |j (setting = 1), between consecutive resets stems therefore from a quench of the transverse field from the pre-quench h 0 = 0 to the post-quench value h [44][45][46][47][48][49][50][51][53][54][55][56][57][58][59][60][61]. The experimental motivation behind the conditional reset protocol comes from recent results on Rydberg quantum simulators [62], which allow to experimentally implement Ising models with transverse and longitudinal fields [63][64][65][66][67][68][69][70] and permit, at the same time, the spatially resolved detection of the spin state upon which the reset is conditioned. Such measurement has been indeed implemented in Rydberg atom experiments on quantum quenches of the Ising model [71][72][73][74].
After the reset, the system evolves again under the quench dynamics up to the time at which the next measurement is performed. The time between consecutive measurements is drawn from a waiting time probability density, p(τ ), which is normalized as ∞ 0 dτ p(τ ) = 1. Repeating this procedure leads at long times to a NESS, which in the case of the TFIC exhibits a sharp crossover between two stationary phases, as sketched in Fig. 1(b). Resetting can thus be interpreted as a dissipative process which leads to an effective open-system evolution.
The reset protocol also allows to suppress naturally occurring dissipative processes, e.g. spontaneous photon emission or scattering in cold-atom experiments. This is achieved by choosing the probability p(τ ) such that the reset must occur before a finite maximum reset time t max , which is shorter than the dissipation timescale. This may allow to reduce decoherence and heating and thus enables one to tailor states of matter with robustness and long-term stability. In the following we will develop the theory for deriving the form of this NESS and we discuss specifically the case of the TFIC. Note, however, that although the present discussion is based on two reset The time evolution is obtained upon alternating the quench-unitary dynamics, schematized by the "Q" boxes, and measurement-reset events, schematized by the "M-R" boxes, which happen at randomly selected times from a waiting time distribution. The average over all the ensuing "reset trajectories" generates the NESS density matrix ρness. In the inset (purple rectangle) the time evolution ensuing in a single quench dynamics and reset is shown. Conditioned on the measured value of the magnetization being positive or negative, the system is reinitialized either in a reset state completely polarized upwards (|1 = |↑↑ . . . ↑ ) or downwards (|2 = |↓↓ . . . ↓ ), respectively. (b) Sketch of the squared magnetization density m 2 ness in the NESS as a function of the transverse field h [see Eqs. (1) and (9)]. It displays remnants of the equilibrium quantum phase transition at hc = 1 through a sharp crossover between a ferromagnetic and a paramagnetic stationary phase.
states, our derivations can be generalized to the case of an arbitrary number of reset states (see Supplemental Material [75]).
Non-equilibrium stationary state When a reset takes place, the system is reinitialized in the reset state |k (k = 1, 2) dependently on the sign of the measured magnetization, see Fig. 1(a). Which reset state is chosen depends only on the time τ elapsed since the previous reset event and on its outcome |j (j = 1, 2), but, fundamentally, not on the previous history of the dynamics. The conditional reset dynamics can be therefore understood as a semi-Markov process [30,43], which is characterized by the transition probabilities P jk (τ ) (obeying k P jk (τ ) = 1) from the previous reset state |j to the following reset state |k , given that a time τ has elapsed between the two resets. The formalism of semi-Markov processes has been applied in Ref. [30] to Markovian open quantum systems. Here we consider, instead, the case of an underlying unitary time evolution, and for the conditional reset dynamics sketched in Fig. 1 we write where P k denotes the projector onto the set of states of the Hilbert space with positive (k = 1) or negative (k = 2) magnetization [76]. To describe the time evolution of the density matrix ρ(t) we exploit the fact that the resetting protocol induces a renewal structure, in the sense that the dynamics after a certain reset does not depend on what happened before that reset. This observation allows to derive "renewal equations" which express the dynamics of the state ρ(t) in the presence of resets, in terms of the state ρ free,j (t) = |ψ j (t) ψ j (t)| which evolves under a dynamics where resets are absent. We stress that in our analysis the resetting is coupled to the underlying unitary dynamics via Eq. (2). This makes the conditional reset dynamics fundamentally different with respect to resets to a single state, first analyzed in Ref. [33] for closed quantum systems, where the reset is independent from the underlying dynamics, which significantly simplifies the analysis.
We consider the observation time t to be fixed and calculate the average density matrix ρ n (t) over all the "reset trajectories" that have exactly n ≥ 0 resets within the interval [0, t]. The density matrix is then given by accounting for all the reset trajectories with an arbitrary number of resets ρ(t) = ∞ n=0 ρ n (t), as shown pictorially in Fig. 1(a). The previous expression can be resummed by taking the Laplace transform ρ(s) of ρ(t). The average density matrix ρ n (s) is, indeed, obtained starting from the initial state, similarly as in Markov chains, by considering n reset transitions R n (s), with the transition matrix R jk (s): In the previous equation, R jk (τ ) = P jk (τ )p(τ ) is the probability density that the system is reinitialized in the reset state |k in the time interval (τ, τ + dτ ) since the previous reset to the reset state |j . Note that R jk (s = 0) is a Markov matrix since k R jk (s = 0) = 1. The resulting expression for ρ(s) is (see Supplemental Material [75]) with and analogously for ρ 02 (s) in terms of ρ free,2 (t). In Eq. (5) we have introduced the survival probability q(τ ) = ∞ τ dτ p(τ ), which is the probability that no reset happens before time τ . The case of Poissonian resetting with constant rate γ corresponds to p(τ ) = p γ (τ ) = γ exp(−γτ ) and q(τ ) = q γ (τ ) = exp(−γτ ). The derivations in this manuscript apply, however, to arbitrary distributions p(τ ).
The first term in Eq. (4), involving ρ 01 (s), accounts for trajectories where no reset takes place and it is therefore given by the unitary time evolution from the initial state |ψ 0 = |1 over the interval [0, t]. The coefficients r 11 (t) ≥ 0 and r 12 (t) ≥ 0, which are obtained from the inverse Laplace transform of r 11 (s) = ∞ n=1 ( R n (s)) 11 and r 12 (s) = ∞ n=1 ( R n (s)) 12 in Eq. (4), can be interpreted as two effective time-dependent rates of jumping into the reset states |1 and |2 , respectively. Crucially, r 11 (t) and r 12 (t) are determined both by the resetting distribution p(τ ) and by the Hamiltonian dynamics, which are coupled due to the conditional protocol. The calculation of the NESS density matrix ρ ness can then be directly performed from Eq.(4) as ρ ness = lim t→∞ ρ(t) = lim s→0 s ρ(s), which gives: This equation is the first main result of this work. It expresses ρ ness as a statistical mixture of the unitary (reset-free) time evolution ensuing from the reset states |1 and |2 . Fundamentally, both the weights of the terms in the sum couple the Hamiltonian dynamics with the reset via Eqs. (2) and (3). One can further check (see Supplemental Material [75]), as noted also in Ref. [33] for resets to a single state, that ρ ness has non-vanishing offdiagonal matrix elements in the basis of the eigenstates of the Hamiltonian H. This shows that ρ ness is a genuinely non-equilibrium density matrix since it is different from the diagonal ensemble, which describes the relaxation at long times of closed quantum systems [33,77].
A complementary way to the renewal equation approach, is to construct an evolution equation for ρ(t), which for classical systems has been done in Ref. [5]. This perspective is useful as it sheds light on the effective open system dynamics which emerges upon superimposing resetting to the unitary time evolution of the system. For our conditional reset protocol we find (see Supplemental Material [75]) from Eq. (4) the generalized Lindblad equation with the time-dependent rates r 1k (t), with k = 1, 2, written in terms of the density matrix (see Supplemental Material [75]) as  Equations (7) and (8) are our second main result. They show that stochastic resetting is generating an engineered dissipative process which leads to an emergent open dynamics for the system without an actual environment being present. The dynamics in Eqs. (7) and (8) is in general non-Markovian because of the presence of the integral over the complete time history of the process. This is a consequence, as in the classical case of Ref. [5], of the fact that for non-Poissonian distributions p(τ ) one has to keep track of the time elapsed since the last reset. Only in the Poissonian case, p γ (τ ), resetting occurs at constant rate γ, independently of the time elapsed since the last reset. In this case one verifies that r(t) = γδ(t) in Eq. (7), which gives back a Markovian-Lindblad dynamics [27-29, 31, 32, 75].
The quantum Ising chain in a transverse field We particularize in this Section our theory to the TFIC introduced above [see Eq. (1)]. Periodic boundary conditions will be henceforth assumed. The quantum critical point at h = h c = 1 separates [78][79][80][81], at zero temperature and in the thermodynamic limit N → ∞, the paramagnetic (h > 1) from the ferromagnetic (h < 1) phase. In the latter, because of the spontaneous breaking of the Z 2 symmetry, the ground-state expectation value of the magnetization density m becomes non-zero Because of the Z 2 symmetry, the matrices P jk (τ ) and R jk (s) in Eqs. (2) and (3) are symmetric. This, in turn, implies that the magnetization density in the NESS in Eq. (6) is identically zero: Tr(m ρ ness ) = 0 (and the same for all the odd operators under the Z 2 symmetry; see Supplemental Material [75]). The natural order param-eter to distinguish between the ordered and the disordered phase is then the squared magnetization density m 2 ness = Tr(m 2 ρ ness ). For the latter one can derive, in the thermodynamic limit N → ∞, the expression (see Supplemental Material [75]) which directly connects the NESS expectation value of m 2 ness in the presence of reset with the magnetization density quench dynamics m(τ ) = σ x (τ ) (with the last equality coming from the choice of periodic boundary conditions). Equation (9) can be evaluated exactly in the thermodynamic limit upon extracting σ x (τ ) from the cluster decomposition principle of the two-point function σ x n σ x n+l (τ ) at large separation l (see Supplemental Material [75]), as outlined first in [82][83][84][85] for the equilibrium correlation function and in [47][48][49] for the non-equilibrium quench dynamics. In the latter case, it is known [47][48][49]61] that σ x (τ ) ∝ exp(−tΛ(h, h 0 )) relaxes exponentially to zero at long times with a characteristic time scale Λ(h, h 0 ) −1 > 0 set by the quench (see Supplemental Material [75]).
In Fig. 2 we evaluate Eq. (9) for the waiting time distribution , which represents a "chopped Poissonian distribution" with rate γ and maximum reset time t max [Θ(τ ) denotes the Heaviside step function]. The NESS reflects the equilibrium quantum phase transition through a sharp crossover due to the competition between the maximum reset time t max and the order parameter decay time Λ(h, h 0 ) −1 . Upon increasing Jt max with fixed γ/J (or, vice-versa, decreasing γ/J with Jt max fixed) the crossover becomes very sharp, as shown in Fig. 2(b) for γ = 0.1J and Jt max = 20. We emphasize that m 2 ness is always different from zero, even for large values of h deep in the paramagnetic phase, for any finite value of Jt max or γ/J. This means that ρ ness , thanks to the resetting protocol, retains part of the long-range order stored in the reset states |1 and |2 . For the reset-free quench time evolution, indeed, the diagonal ensemble predicts a vanishing expectation value of m 2 and therefore a complete loss of the order of |1 and |2 . As Jt max → 0 one encounters the quantum Zeno effect [86,87] since the initial state does not evolve in time due to very frequent resets and therefore the squared magnetization density remains close to 1. In the limit of large Jt max (small γ/J), when the crossover gets sharp and it closely resembles the equilibrium quantum phase transition, a large time, on average, elapses between consecutive resets. As a consequence, one can evaluate Eq. (9) using the asymptotics of σ x (τ ) ∝ exp(−tΛ(h, h 0 )) for long times τ → ∞, which has been analytically derived in Refs. [47][48][49]75]. The result is shown in Fig. 2(d) (see Supplemental Material [75]). We can see that the asymptotics of the order parameter is in excellent agreement with the exact numerical evaluation of Eq. (9). Only for values of h 0.9 discrepancies are visible. This is caused by the fact that when h becomes large the order parameter σ x (τ ) quickly relaxes to zero as a function of τ and the integral in Eq. (9) is therefore no longer dominated by the behavior of σ x (τ ) at large times.
Outlook. We have presented a general protocol for constructing the NESS of a unitary many-body dynamics interspersed by random resets. Using the TFIC as exemplary case, we have shown that this NESS is intimately connected to the dynamics following a quantum quench, and that it embodies the ground-state quantum phase transition of the TFIC through a sharp crossover. NESS with sharp crossovers can here be exploited to design collectively enhanced sensing protocols where a small change of the parameter (external perturbation) driving the quantum phase transition translates into an easy-to-detect macroscopic response, as, e.g., discussed in Refs. [88,89]. A possible advantage of our protocol is, however, that the establishment of the NESS relies on stochastic resetting, which diminishes the impact of uncontrolled dissipative effects such as heating.

Supplemental Material
Designing non-equilibrium states of quantum matter through stochastic resetting Here we derive in detail the results presented in the main text. The Supplemental Material is organized as follows: in Section I we show how the conditional resetting dynamics can be formulated in full generality using the formalism of semi-Markov processes S1 ; in Section II we particularize the analysis to the quantum Ising chain in a transverse field (TFIC), for which we show how to compute the phase diagram of the magnetization density squared m 2 ness in the non-equilibrium steady state (NESS).

I. CONDITIONAL RESET IN ISOLATED MANY-BODY QUANTUM SYSTEMS
In this Section we provide the general derivation of Eqs. (3)-(8) of the main text. In particular, in Subsec. I A we show how the NESS density matrix ρ ness in Eq. (6) of the main text can be analytically constructed using the renewal equation approach. In Subsec. I B we follow a complementary approach based on the corresponding evolution equation for the density matrix ρ(t) as a function of time in the presence of reset. In particular, we obtain the generalized non-Markovian Lindblad equation in Eq. (7) of the main text.

A. Renewal equation approach and construction of the NESS
We consider a many-body quantum system composed by N particles/spins, which undergoes unitary time evolution according to a certain Hamiltonian H. The Hamiltonian in this Section will be left general. In Sec. II we will consider the specific case of the TFIC Hamiltonian. The unitary time evolution is interrupted by measurements of some observables. The time between consecutive measurements is distributed according to the waiting time probability density p(τ ), i.e., after a measurement the next measurement occurs in the time interval (τ, τ + dτ ) with probability p(τ )dτ . The waiting time probability density is normalized to one ∞ 0 p(τ )dτ = 1. (S1) The following derivation applies for arbitrary p(τ ). To account, however, for the finite coherence time attained in cold-atom systems, we will be particularly interested in the "chopped exponential distribution" p γ,tmax (τ ) with rate γ and maximum reset time t max , which we define as with Θ denoting the Heaviside step function. It is useful to consider the two limiting cases of Eq. (S2) for t max → ∞ (with fixed γ) and for γ → 0 (for fixed t max ). In the former case we have that p γ,tmax (τ ) reduces to the Poisson (exponential) distribution: In the limit γ → 0, on the other hand, p γ,tmax (τ ) reduces to a uniform distribution p tmax (τ ) in the interval (0, t max ): We define for future convenience the survival probability q(τ ) which gives the probability that no measurement has been performed before time τ . In the case of p γ,tmax one has In the limit t max → ∞ from Eq. (S3) one has When a measurement is performed, the many-body wavefunction of the system is reset to a state chosen within a set of "reset states". We consider in this Section, for concreteness, the case of a set of reset states composed by two elements. As we will comment later on, the generalization to an arbitrary number of reset states is straightforward and it does not bear additional conceptual difficulties. We consider the two following reset states We denote with |↑ = 1 √ 2 (1, 1) and |↓ = 1 √ 2 (1, −1) the eigenstates of the Pauli matrix σ x . Both the states in Eq. (S9) present ferromagnetic long-range order and, indeed, they are the two degenerate ground states of the TFIC Hamiltonian in Eq. (S60) for h = 0, as we shall see in Sec. II. We then define the sets (S10a) (S10b) which are nothing but the sets of eigenstates of the magnetization operator N n=1 σ x n /N with zero, positive and negative eigenvalues, respectively. Note that for an odd number N of spins C 0 is just the empty set. The sets C 0 , C 1 and C 2 define a partition of the whole Hilbert space H with C 0 ∪ C 1 ∪ C 2 defining a complete basis of H and In addition, the orthonormality condition is satisfied by the elements of C 0 , C 1 and C 2 . Clearly |1 ∈ C 1 , while |2 ∈ C 2 . When a reset event happens, the system is re-initialized in one of the two states in Eq. (S9) according to the probability with j, k = 1, 2 and the summation in µ running over the elements of the set corresponding to the index k or 0 according to Eq. (S10). In Eq. (S12) we have defined the operator P k corresponding to the set C k as and the unitary evolution from the state |j since C 0 ∪ C 1 ∪ C 2 is a complete basis of the Hilbert space H, as already commented after Eq. (S10). On the basis of the wavefunction interpretation of quantum mechanics S2 , we, therefore, see that P jk (τ ) in Eq. (S12) represents the probability of measuring a non-negative (k = 1) or non-positive (k = 2) value of the magnetization after a time τ since the previous reset towards the reset state j. The factor 1/2 in front of the second summation of Eqs. (S12) and (S13) is chosen such that if a zero magnetization is measured, the many-body state is reset with probability 1/2 to the state |1 and with probability 1/2 to the state |2 . We also remark that the operator P k in Eq. (S13) is not a projector. This, however, is not necessary for the following derivation to hold. Indeed, only the normalization condition for P jk (τ ) in Eq. (S15) is necessary for the semi-Markov construction to be well-defined. From Eq. (S12), as a matter of fact, one realizes that the outcome of a reset transition depends on the time τ elapsed since the last reset and on the outcome j of the previous reset, but, fundamentally, not on the previous history of the dynamics. This kind of dynamics precisely defines a semi-Markov process (see, e.g., Ref. S1). The semi-Markov structure induced by Eqs. (S12)-(S15) is the essential ingredient allowing for an exact, analytical, construction of the NESS, as we now show. We combine the waiting time probability density p(τ ) and the probability P jk (τ ) in Eq. (S12) into the following probability which satisfies, for any value of j, the following normalization condition from Eqs. (S1) and (S15) Without loss of generality, we take the initial state of the system at time t = 0 to be the reset state |1 in Eq. (S9).
We emphasize, however, that the specific choice of the initial state |ψ 0 influences only the finite time, transient, dynamics. It is, indeed, straightforward to generalize the present derivation to an arbitrary initial state |ψ(0) and, thereby, to show that the steady state density matrix ρ ness , in Eq.(6) of the main text, does not depend on the choice of |ψ(0) . We have now defined all the quantities necessary for the study of the state-dependent reset dynamics. One is now interested in deriving an equation for the density matrix ρ(t) under such kind of dynamics. To do this we consider the observation time t to be fixed and we write the density matrix ρ n (t) at time t assuming that exactly n ≥ 0 resets happened in the interval (0, t). We write the first terms to provide intuition on the emerging general structure: ρ 0 (t) = q(t)ρ free,1 (t), n = 0, (S19a) with and the definition from Eq. (S14) the reset-free, unitary, evolution starting from the state |j . The physical interpretation of Eq. (S19) is clear. The first equation for ρ 0 (t) simply accounts for the case of no reset up to time t and therefore it is analogous to the reset-free evolution weighted by the survival probability q(t). The second equation for ρ 1 (t) accounts for the case of only one reset happening at time τ 1 ∈ (0, t) with probability p(τ 1 ). For the remaining time t − τ 1 the system evolves without additional resetting and therefore the weight q(t − τ 1 ) is included. The convolution structure in Eq. (S19) suggests using the Laplace transform, which is defined for an arbitrary function f (t) of time t as: Taking the Laplace transform of ρ n (t) in Eq. (S19) one obtains with the definition (cf., Eq. (S21)) and the matrix R jk (s) the Laplace transform of the matrix R jk (τ ) in Eq. (S16) (with the same notation as in Eq.(3) of the main text) The matrix R jk (s) ≥ 0 is non-negative. From Eq. (S23) in Laplace space one then has which comes from recognizing the matrix product structure in Eq. (S19). In the Laplace space the matrix R(s) satisfies the following relations for every row j From the last equality one recognizes that for s = 0 the matrix R(s) is a stochastic-Markovian matrix S1 . The geometric sum in Eq. (S26) is convergent provided the spectral radius, i.e., the absolute value of the largest eigenvalue of R(s), is smaller than 1. From Eq. (S27) S3, S4 , we see that this is true for s > 0: We can explicitly invert the matrix R(s) since it is of dimension 2 × 2 The inverse Laplace transform r 11 (t), r 12 (t), r 21 (t) and r 22 (t) of the previous equations can be performed in general numerically. The expression for r 11 (s) and r 12 (s) can then be inserted into Eq. (S26) attaining which can be readily inverted to have the desired equation for ρ(t): The previous equation can be solved to get ρ(t) once r 11 (t), r 12 (t) ≥ 0 are known from the inverse Laplace transform of Eq. (S29). Note that r 11 (t), r 12 (t) are non-negative (and the same for r 21 (t) and r 22 (t)) since they are obtained from the inverse Laplace transform of the series in Eq. (S28), with the matrix R jk (τ ) ≥ 0 from the definition in Eq. (S25). These two coefficients can be therefore considered as the effective time-dependent rates of jumping into the reset state 1 or 2, respectively, at time t having started at time 0 from the state |1 . Equation (S31) can be seen as the equivalent of the so-called "last renewal equation" of reset process, see, e.g., Refs. S5-S7 for classical systems and Refs. S8 for quantum systems. The fundamental complication with respect to the case of state-independent resets analyzed in S8 is encoded into the rates r 11 (t), r 12 (t) which depend on the conditional reset mechanism, in our case see Eq. (S12).
Note that the first term on the right hand side of Eq. (S31), as well as the rates r 11 (t) and r 12 (t), depends on the choice of the initial state that we have done in Eq. (S18). The stationary state, on the contrary, does not depend on this choice. The steady state limit of Eq. (S31) can be obtained from the "final-value theorem" of Laplace transforms S9 , which for an arbitrary function of time f (t) reads which gives the steady state limit by locating the singularity of f (s) at s = 0 (which must be a simple pole for the steady state to exist). From Eq. (S27), one immediately sees that for s = 0 one eigenvalue of R(s = 0) is equal to 1, which causes the series in Eq. (S28) to diverge. The rates r 11 (s) and r 12 (s) (and the same for r 21 (s) and r 22 (s) which are not reported for the sake of brevity) in Eq. (S29) accordingly become singular at s = 0 with a leading simple pole behavior behavior which inserted into Eq. (S26) with Eq. (S32) gives the steady state density matrix ρ ness This coincides with Eq. (6) of the main text and it is one of the main results of this manuscript. We emphasize that the steady state ρ ness is determined solely by the leading singular behavior of the rates around s = 0 in Eq. (S33), which is independent on the choice of the initial state in Eq. (S18). From Eq. (S34) one has access to the steady state average of an arbitrary observable O in terms of its free-unitary dynamics as with the notation from Eqs. (S14) and (S21) In concluding this Subsection we mention that the weights of the terms in the sum in the right hand side of Eq. (S34) have a clear interpretation in terms of the underlying semi-Markov process. As a matter of fact, one immediately realizes that the vector is a left eigenvector of R(s = 0) with eigenvalue 1. The entries P stat and P stat of the vector P stat are nothing but the steady states occupation probabilities of a two-state Markov chain S3 . On the basis of this observation it becomes also evident that the generalization of the derivation of ρ ness to an arbitrary number n of reset states, |1 , |2 . . . |n , is straightforward. Equation (S34) generalizes to this case by embodying the summation over n terms corresponding to the reset-free evolution starting from each of the n reset states. The weight of each term of the sum will be then given by the steady state occupation probability of the associated n-state Markov chain identified by the n × n matrix R jk (s = 0), j, k = 1, 2 . . . n, in the Laplace s domain.

B. Generalized master equation approach and the Lindblad equation
We now show that Eq. (S31) for ρ(t) can be seen as the solution of an evolution equation for dρ/dt. To do this, we first define the operator L: which is the generator of the reset-free unitary time evolution ρ free,i (t) of the density matrix. One, then, proceeds from Eq. (S30), which we rewrite as with Since ρ(0) = |1 1| from Eq. (S18), one can invert Eq. (S41) to the time domain eventually obtaining Note that the rates r 11 (t), r 12 (t) have to be related to r(t) in order for Eq. (S43) to be trace-preserving, in particular, one needs the constraint which can be readily checked in Laplace space as with simple algebraic manipulations starting from Eq. (S29). The first term on the right hand side of equation (S43) describes the unitary dynamics between reset events according to Eq. (S38). The second and third terms, containing the effective rates r 11 (t) and r 12 (t), can be seen as gain terms due to resetting trajectories jumping into the states |1 and |2 . The last term on the right hand side, on the other hand, can be considered as a loss term due to resetting trajectories where the systems evolves unitarily in the time interval t − t and it is then subject to resetting. Equation (S43) can be therefore considered as a generalized master equation. The derivation here adopted extends the one of Ref. S10, which regards state-independent reset processes, by including the effect of the conditional reset protocol within the rates r 11 (t), r 12 (t) (analogously to the discussion done after Eq. (S31)). The latter can be written in terms of the density matrix as and the operator P k defined in Eq. (S13). Equation (S46) can be proved by replacing ρ(t ) inside the integral with Eq. (S31) and integrating each term. In particular, the integration of the first term on the right hand side of Eq. (S31) gives where in the second equality we used the definitions in Eqs. (S21) and (S38) for the reset-free, unitary, dynamics, and Eq. (S42). The integral of the second and third term of Eq. (S31) can be similarly performed In Eq. (S52) each summation on the index µ runs over the range specified by Eq. (S53) and {. . . , . . . } denotes the anticommutator. Notice that in Eq. (S52) we used the definition of P k in Eq. (S13) and the fact that The case of Poissonian resetting at constant (time-independent) rate γ is in this perspective special. In this case, on the basis of the definition of a Poisson process, a reset occurs in the time interval (τ, τ + dτ ) with probability γdτ , without the need to know when the previous reset occurred. Particularizing Eq. (S42) to the exponential waiting-time distribution p γ (τ ) = γexp(−γτ ) one, indeed, has and therefore Eq. (S52) simplifies to which is in the Lindblad form S11,S12 with the jump operators √ γL a µ → L a µ rescaled by the resetting rate γ. Equation (S56) is, fundamentally, Markovian since it is local in time, differently from Eq. (S52). This derivation extends the analysis of Ref. S13-S16, where only state-independent (no conditional reset mechanism) reset was considered. Fundamentally, in Ref. S13-S16 the Lindblad equation is assumed as a starting point for the analysis of Poissonian resetting at rate γ. The analysis of this manuscript, on the contrary, shows exactly how the averaging of the microscopic unitary evolution over the reset trajectories leads to an effective non-Markovian generalized Lindblad equation (S52). In the specific case of the exponential waiting time probability density p γ (τ ), our analysis constitute an exact derivation of the Markovian Lindblad equation for an open quantum system subject to stochastic resetting.
We mention that for the exponential waiting-time distribution one can simulate the conditional reset dynamics with the stochastic wave function approach S17 , where the probability Prob jk of making a jump (reset) towards the state |k , being the last jump (reset) to |j , is where | ψ j (τ ) = exp(−iH eff τ ) |j , with H eff denoting the effective, non-Hermitean, Hamiltonian which dictates the time evolution | ψ j (τ ) between consecutive quantum jumps (resets) S17 and we have used the definitions in Eqs. (S12) and Eq. (S25).

II. STOCHASTIC RESETTING IN THE QUANTUM ISING CHAIN IN A TRANSVERSE FIELD
In this Section we particularize the general theory of Sec. I to the TFIC. Our motivation to investigate the TFIC model is twofold: first, it is a paradigmatic model of (quantum) statistical mechanics which displays a quantum phase transition; second, it allows for an exact determination of the NESS with reset in the thermodynamic limit. The availability of such an analytical result is particularly useful as a reference for the analysis of more complex models, such as Rydberg atoms arrays S18 , where a full analytical treatment is not possible. In Subsec. II A we briefly recapitulate some well known equilibrium properties of the quantum Ising chain which will be essential for the understanding of our results. In Subsec. II B we discuss the dynamics of the TFIC for the conditional reset protocol of Sec. I. In Subsec. II C we give details of the evaluation of the phase diagram for m 2 ness , which is reported in Fig. 2 of the main text. In Subsec. II D we eventually detail the derivation of the asymptotics of m 2 ness for large Jt max , low γ/J, which is reported in Fig. 2(d) of the main text.

A. Equilibrium properties of the Ising chain in a transverse field
The Hamiltonian of the model for a total number N of spins is given by where we assume periodic boundary conditions σ x n+N = σ x n . Various notational conventions can be adopted for the TFIC, here we follow the one of Refs. S19-S21. The ferromagnetic coupling J > 0 gives an overall energy scale. Time is accordingly measured in units of J −1 . The quantum Ising chain is a paradigmatic example of a quantum phase transition happening for the critical magnetic field value h c = 1 at zero temperature T = 0. The paramagneticdisordered phase is attained for h > 1, while the ferromagnetic-ordered phase for h < 1. In the latter phase the order parameter of the quantum phase transition becomes non-zero. The quantum phase transition comes with the spontaneous breaking of the Z 2 symmetry of rotations by 180°around the z axis The quantum Ising Hamiltonian in Eq. (S60) possesses the Z 2 symmetry, which is implemented by the unitary operator U : We define lattice fermion operators using the Jordan-Wigner transformation (see, e.g., Refs. S22 and S23) in terms of which the operator U in Eq. (S63) can be written as The unitary operator U counts, therefore, the parity of the Jordan-Wigner fermions c n . As a consequence of the Z 2 symmetry, H and U can be diagonalized simultaneously and the Hamiltonian is block diagonal where H e and H o denote the Hamiltonian restricted to the even and odd fermion sector of the Hilbert space, respectively. The exact solution of the model in terms of free fermionic quasi-particles is extremely well known, see, e.g., Refs. S22 and S23, and we therefore do not report it here for the sake of brevity. It is sufficient for our purpose to report the single-quasi-particle energy spectrum ε h (k) The ground state of H e in the even fermion sector will be denoted henceforth as |0, h N S (the subscript NS stands for Neveu-Schwarz, which is the name conventionally given to the even fermion sector S19-S21 ). The ground state of H o in the odd fermion sector will be denoted as |0, h R (the subscript R stands for Ramond which is the name given to the odd fermion sector). Note that which shows that the NS and R vacua are eigenstates of the parity operator U, i.e., they do not break the Z 2 symmetry. For h > 1, in the paramagnetic phase, it can be shown S19-S21 that |0, h N S is the ground state of the complete Hamiltonian H in Eq. (S66) since it has an energy lower than the one of |0, h R . For h < 1, in the ferromagnetic phase, the spontaneous breaking of the Z 2 symmetry is reflected by the fact that, in the thermodynamic limit N → ∞, the states |0, h N S and |0, h R become degenerate. As a consequence in the thermodynamic limit, for h < 1, either the symmetric |GS(h) + or the anti-symmetric |GS(h) − combination of |0, h N S and |0, h R is chosen as a ground state. In formulas, the ground state of the system will be denoted as |GS(h) and it is then given by Note that from Eq. (S68) Namely, the ground state for h < 1 is not an eigenstate of the fermion parity, thereby expressing the breaking of the Z 2 symmetry. Spontaneous symmetry breaking is also detected in the ground-state expectation value of the order parameter m in Eq. (S61), which is non-vanishing in the symmetry-broken phase for h < 1 S22,S24-S27 In the first equality of the previous equation we used the fact that for periodic boundary conditions the chain is translationally invariant and therefore the vacuum expectation value ± GS(h) | σ x n | GS(h) ± does not depend on the lattice site n. We conclude this subsection by defining even O e and odd O o operators under the Z 2 transformation in Eq. (S62) as with the plus sign applying to even operators, while the minus sign to odd operators. The distinction between even and odd operators is crucial for the study of the dynamics of the TFIC from the reset states |1 and |2 in Eq. (S9) according to the formalism explained in Sec. I. This is the problem that we address in the next Subsection.

B. Construction of the NESS for the quantum Ising chain
In this Subsection we construct the steady state density matrix ρ ness and the corresponding stationary expectation values according to the general theory of Sec. I (see Eqs. (S34) and (S35)). We do this by connecting our theory with quantum quenches in the Ising chain. We first note that the choice of the reset states |1 and |2 in Eq. (S9) allows for a clear connection with quantum quenches since |1 = |GS(0) + and |2 = |GS(0) − are the two degenerate ground states of the Ising Hamiltonian H(h 0 ) for h 0 = 0, according to Eq. (S70). The unitary time evolution between two consecutive reset events is therefore a quench of the transverse field S19-S21,S28-S31 from the pre-quench value h 0 = 0 to the post-quench value h. To make contact with the notation introduced in the previous Subsec. II A we write the reset states as follows We observe that, for h 0 = 0, the spin representation of the fermionic NS and R vacua is simple, as it corresponds to the superposition of the two eigenstates completely polarized in the x-longitudinal direction of the Hamiltonian in Eq. (S60). For generic h 0 = 0 the representation of the NS and R vacua in the spin language cannot be simply written. This also motivates our choice for the reset states |1 = |GS(0) + and |2 = |GS(0) − , since they have a simple representation in the spin basis as product states. This makes the experimental realization of such states simpler. We, however, observe that the following discussion can be immediately generalized , upon writing |1 = |GS(h 0 ) + and |2 = |GS(h 0 ) − in Eq. (S73), to the case where the two reset states are taken as the two degenerate ground states for a non-zero value of the pre-quench transverse field h 0 . For this reason, in all the following formulas, we denote with h 0 the pre-quench transverse field albeit in the numerical evaluation we have always used the specific value h 0 = 0. We now use Eq. (S73) to rewrite Eq. (S34) in a simpler way. First of all because of the Z 2 symmetry, the expression we have chosen for the probability P jk (τ ) in Eq. (S12) renders and therefore the matrix R(s) in Eq. (S25) has a very simple structure Using Eq. (S75) into Eq. (S34) one has where ρ free,1 (τ ) (ρ free,2 (τ )) is the quench-unitary time evolution over the time interval τ between two consecutive resets. It is important to emphasize that Eqs. (S74)-(S76) are a consequence not only of the Z 2 symmetry of the Hamiltonian (S60), but also of the chosen reset states in Eq. (S9) and of the form for the probability P jk (τ ) in Eq. (S12). We now want to draw some general statements regarding ρ ness in Eq. (S76) for the case of the waiting time probability density p γ,tmax (τ ) (with the associated survival probability q γ,tmax (τ )) in Eq. (S2). In particular, the limiting case of the exponential distribution (S3) t max → ∞ (with fixed γ) in Eq. (S3) has been already considered in Ref. S8, where it has been shown that ρ ness retains non-vanishing off-diagonal matrix elements in the basis of the eigenstates of the post-quench Hamiltonian H(h). We now show that the same applies in the more general case of the waiting time p γ,tmax (τ ) for any finite value of γ and/or t max . One can prove this by considering Eq. (S76) for γ = 0. In this limit the matrix elements of ρ ness in the basis of the eigenstates of the post-quench Hamiltonian H(h) are (ρ ness ) α,β = 1 2 In the previous equation c α = α | i denotes the overlap coefficients between the reset state |i and the eigenstates H(h) |α = E α (h) |α of the post-quench Hamiltonian. Inserting Eq. (S78) into Eq. (S77) and performing the integration one obtains It is then clear that for any finite value of t max the off-diagonal matrix elements of ρ ness are different from zero. The steady state density matrix becomes diagonal in the basis of the post-quench eigenstates only as t max → ∞ (ρ ness ) αβ = 0 α = β, Plugging Eq. (S80) into Eq. (S76) one has which coincides with the diagonal ensemble. In the last step we used that |c 1 | 2 = |c 2 | 2 because of the Z 2 symmetry. We therefore conclude that the steady state density matrix ρ ness always retains non-vanishing off diagonal matrix elements in the basis of the eigenstates of the post-quench Hamiltonian H(h) and it is therefore different from the diagonal ensemble, which describes the relaxation of closed quantum systems at long times S32 . Only in the case of γ → 0 and t max → ∞ the off-diagonal matrix elements vanish and the steady state density matrix in the presence of reset coincides with the diagonal ensemble, according to Eqs. (S80) and (S81). This is expected because in the limit γ → 0 and t max → ∞, on average, a large time between consecutive resets elapses and therefore the dynamics is dominated by the long time asymptotics of the quench-unitary dynamics. The steady state density matrix ρ ness is therefore a non-equilibrium density matrix which is substantially different from the diagonal ensemble for any finite value of γ and/or t max . The latter difference has clearly important consequences on the steady state expectation value of local observables. Inserting Eq. (S76) into Eq. (S35) one has O ness = Tr(Oρ ness ) = 1 2 It is important to emphasize that the calculation of the expectation value of an odd operator over the symmetrybroken ground state is tremendously hard since it connects the ground state |0, h 0 R in the odd fermion sector to the ground state |0, h 0 N S in the even fermion sector, see, e.g., Refs. S19-S22, S24-S27.
where we used the fact that averages of even operators O e over the ground state |0, h 0 N S are equal to the ones over the state |0, h 0 R in the thermodynamic limit N → ∞. By inserting in Eq. (S85) into Eq. (S82) one has lim N →∞ The calculation of expectation values of even operators in the symmetry broken ground state according to Eq. (S85) can be done by exploiting the Wick theorem-free fermions techniques. Exploiting the results available from Refs. S19-S22, S24-S27 for the two-point correlation function of the order parameter σ x we derive in the next Subsection the exact result for the steady state squared magnetization density in the thermodynamic limit in the presence of resetting.
C. Phase diagram for the magnetization density squared in the NESS Based on Eq. (S84) one has that the NESS expectation value m ness is identically zero according to Eq. (S84), since the operator m in Eq. (S61) is odd under the Z 2 symmetry. To get a valuable order parameter distinguishing the ferromagnetic from the paramagnetic phase, the natural choice is then the squared magnetization density which is even under the Z 2 symmetry. Inserting Eq. (S87) into Eq. (S86) one can express everything in terms of the calculation of equal-time two-points correlation functions of the longitudinal magnetization where in the second line we have used translational invariance, which causes the two-point function to depend only on the distance l between the two spins. We have also denoted, for brevity, m 2 ness = lim N →∞ m 2 ness to avoid writing henceforth the thermodynamic limit every time for the left hand side of Eq. (S88). One has to keep in mind that the following results apply in the thermodynamic limit N → ∞. The summation on the right hand side of Eq. (S88) is, indeed, dominated, as N → ∞, by the thermodynamic limit of the equal-time two-point function of the order parameter Consequently, one can evaluate the sum in Eq. (S88) directly replacing provided the latter limit exists. Replacing Eq. (S88) and then using Eq. (S90), one has where We report in the following all the expressions necessary for the evaluation of m 2 ness as in Fig. 2 of the main text. To compute σ x we follow the analysis of Refs. S24-S27, where one considers the two point function in the large distance limit l → ∞ according to Eq. (S93). The equal time two-point function of the longitudinal magnetization is non-trivial to compute because σ x is a non-local operator in terms of the Jordan-Wigner fermions because of Eq. (S64) where we have introduced the Majorana fermion operators a x n and a y n at site n, which are defined as a x n = c n + c † n , a y n = i(c † n − c n ). Equations (S105) and (S107) have been used in the manuscript to compute ρ xx (l, t) as a function of time t. Upon evaluating ρ xx (l, t) for large values of l (we have used l = 200 in our numerical analysis) one can extract the order parameter σ x (t) 2 as a function of time from the clustering decomposition formula in Eq. (S93). The resulting expression for σ x (t) 2 can be eventually inserted in Eq. (S91) to get m 2 ness . The result is shown in Fig. 2 of the main text. We emphasize that the determinant representation of the correlation function in Eq. (S107) works for arbitrary values of both the pre-quench h 0 and the post-quench h field. It can be therefore used both to study quenches within the ferromagnetic/paramagnetic phases and to study quenches across the two phases. For the case of interest in this manuscript, the pre-quench field is always fixed to zero, h 0 = 0 (cf., Eq. (S73)), while h can assume arbitrary values. With our conditional reset protocol starting from the reset states in Eq. (S73), the unitary time evolution is therefore described by either quantum quenches within the ferromagnetic phase (h 0 = 0 and h < 1), or by quenches from the ferromagnetic to the paramagnetic phase (h 0 = 0 and h > 1).
D. Asymptotics of the magnetization density squared for large Jtmax and low γ/J In this Subsection we eventually discuss the asymptotic behavior of the m 2 ness for large Jt max and low γ/J. This limit is physically interesting because in this regime the crossover displayed by m 2 ness as a function of h becomes very sharp and it closely resembles the equilibrium quantum phase transition behavior in Eq. (S71). The result of this analysis is reported in Fig. 2(d) of the main text. We provide here all the necessary expressions for the underlying numerical evaluation.
In the limit of low γ/J and large Jt max , on average, a large time elapses between one reset and the subsequent one and therefore the integral in Eq. (S91) is dominated by the long time behavior of σ x 2 (t). This is fundamental from the analytical point of view because the long time behavior of σ x (t) after a quantum quench displays universal features. In particular, for a quench within the ferromagnetic phase, h 0 , h < 1, it has been analytically proved in Refs. S19-S21 that the order parameter σ x (t) at long times t 1 decays exponentially in time as σ x (t) = (C F F ) 1/2 exp (−tΛ(h, h 0 )) , as t → ∞, (S108) with and Λ(h, h 0 ) = − π 0 dk π ε h (k)ln| cos ∆ k | > 0, where ε h (k) = d dk ε h (k). In our case h 0 = 0 and therefore we can restrict ourselves to the consideration of Λ(h, 0), which is given upon performing the integral in Eq. (S110) Inserting Eqs. (S108)-(S111) into Eq. (S91), with the survival probability in Eq. (S6), one has m 2 ness = γC F F (h, 0)e −γtmax 1 − e −γtmax − γt max e −γtmax e γtmax − 1 γ + 2Λ(h, 0) + γ(e −2tmaxΛ(h,0) − 1) 2Λ(h, 0)(γ + 2Λ(h, 0)) for h < 1. (S112) In the case h > 1 (for quenches across the quantum critical point) the formula in Eq. (S108) does not apply and no analytical formula has been so far proved. In Ref. S20 the following conjecture has been put forward for long times t 1 σ x (t) = (C F P ) 1/2 [1 + cos(2ε h (k 0 )t + α) + . . . ] 1/2 exp (−tΛ(h, h 0 )) , as t → ∞, where the . . . denotes subleading corrections. The constant α(h, h 0 ) in Eq. (S113) is an undetermined phase shift, which has to be fixed by keeping it as a fitting parameter. The momentum k 0 is determined such that cos ∆ k0 = 0. The amplitude C F P is Our aim is to use Eq. (S113) to get an estimate of the asymptotic value attained by m 2 ness at large values of h deep in the paramagnetic phase. By plugging Eq. (S113) into Eq. (S91) one obtains at leading order as h → ∞ with Λ(h, 0) = 4J π , for h > 1, (S116) and C F P (h, 0) = 1 from Eq. (S114). Equations (S112) and (S115) are reported in Fig. 2(d) of the main text as dashed lines for h < 1 and h > 1, respectively. A couple of important observations are now in order. First, m 2 ness is always different from zero, also for large values of h, for any finite value of γ and/or t max . This is a direct consequence of the fact that ρ ness is different from the diagonal ensemble, as commented after Eq. (S81). The diagonal ensemble is in the case here of interest equivalent to the generalized Gibbs ensemble, for local operators in the thermodynamic limit, since the TFIC is an integrable model S20 . The generalized Gibbs ensemble, in particular, predicts a vanishing expectation value for the magnetization density in the stationary state attained at long times after a quench. The long-range order of the pre-quench ground state is therefore entirely wiped out by the purely unitary quench time evolution. The non-zero value of m 2 ness that we computed in Fig. 2 is, therefore, an intrinsic property of the reset protocol, which shows that superimposing a controlled source of dissipation on top of the unitary time evolution eventually allows for the storage of part of the long-range order of the states |1 and |2 in the steady state. This information storage is not lost even for large values of h, deep in the disordered phase, where m 2 ness attains a finite limiting value in good agreement with the prediction of Eq. (S115). Second, we note that deviations from the asymptotics behavior predicted by (S112) and (S115) are visible only for h 0.9. In the latter case the discrepancy is caused by the fact that σ x (t) for large quenches, h h 0 = 0, decays extremely rapidly as a function of time t and it is non-zero only for short times. The time integral in Eq. (S91) is therefore not anymore dominated by the long time behavior, where equations (S108) and (S113) apply. In the case h > 1 one has a relative error of 10 − 15% between the asymptotic value in Eq. (S115) at large h and the numerically exact asymptotic value. We again remark that Eq. (S113) is a conjecture, not rigorously proven, and therefore Eq. (S115) can be taken as an estimate of the residual (minimum) magnetization attained by the system for large values of the post-quench field h.