On time crystallinity in dissipative Floquet systems

TL;DR

This paper explores how dissipative Floquet systems can maintain stable subharmonic responses, identifying conditions that preserve time-crystallinity despite noise and dissipation effects, and contrasting with disordered many-body localized systems.

Contribution

It provides a framework for understanding stability of time-crystalline order in dissipative quantum systems, including examples of stable oscillations and distinctions from disordered systems.

Findings

Short-ranged 2D system exhibits persistent subharmonic oscillations.

Dissipation can stabilize time-crystalline order under certain conditions.

Fully connected models show stability of subharmonics in the thermodynamic limit.

Abstract

We investigate the conditions under which periodically driven quantum systems subject to dissipation exhibit a stable subharmonic response. Noting that coupling to a bath introduces not only cooling but also noise, we point out that a system subject to the latter for the entire cycle tends to lose coherence of the subharmonic oscillations, and thereby the long-time temporal symmetry breaking. We provide an example of a short-ranged two-dimensional system which does not suffer from this and therefore displays persistent subharmonic oscillations stabilised by the dissipation. We also show that this is fundamentally different from the disordered DTC previously found in closed systems, both conceptually and in its phenomenology. The framework we develop here clarifies how fully connected models constitute a special case where subharmonic oscillations are stable in the thermodynamic limit.

Figures

Equations

$${\hat{H}}_{\mathrm{TFIM}}=-\sum _{\ell }{\sigma }_{\ell }^z{\sigma }_{\ell +1}^z+g\sum _{\ell }{\sigma }_{\ell }^x,$$

$${\partial }_t\rho =-i[\hat{H}(t),\rho ]+\sum _i\left({\hat{L}}_i\rho {\hat{L}}_i^{†}-\frac{1}{2}\{{\hat{L}}_i^{†}{\hat{L}}_i,\rho \}\right)$$

$$\hat{H}(t)=\{\begin{array}{l}\frac{\theta }{t_R}{\sum }_{\ell }{\sigma }_{\ell }^x;0\le t<t_R\\ {\hat{H}}_{\mathrm{TFIM}};t_R\le t<t_R+t_C.\end{array}$$

$$L_{\alpha }=\sqrt{\gamma }\left[\Theta (m_{\alpha })|\Uparrow ⟩⟨\alpha |+\Theta (-m_{\alpha })|\Downarrow ⟩⟨\alpha |\right],$$

$$C(t)=\frac{1}{N^2}\sum _{\ell \ne r}\mathrm{Tr}\left[\rho (t){\sigma }_{\ell }^z{\sigma }_r^z\right].$$

$$L_{\alpha ,\ell }=\sqrt{\gamma }\left(\prod _{i\ne \ell }{\mathcal{P}}_{i,s_i^{\alpha }}\right)\times \{\begin{array}{l}{\sigma }_{\ell }^x,M_{\{r_{\ell }\},\alpha }=0\\ {\sigma }_{\ell }^{\pm },\mathrm{sgn}\left(M_{\{r_{\ell }\},\alpha }\right)=\pm 1\end{array}$$

$${\partial }_t{\rho }_{\alpha \alpha }=\sum _{\beta }{\gamma }_{\alpha \beta }{\rho }_{\beta \beta }-\left(\sum _{\beta }{\gamma }_{\beta \alpha }\right){\rho }_{\alpha \alpha }$$

$${\rho }_{\alpha \beta }={\varrho }_{m{m}^{\prime }}\forall (\alpha ,\beta )s.t.m_{\alpha }=m\mathrm{and}{m}_{\beta }=m^{\prime }.$$

$$\begin{array}{rl}⟨m⟩& =\sum _{m}m{\varrho }_{mm}{\mathcal{N}}_m,\\ C& =\frac{1}{2}\sum _m\left(m^2-\frac{1}{N}\right){\varrho }_{mm}{\mathcal{N}}_m,\end{array}$$

$$\begin{array}{rl}{\partial }_t{\varrho }_{m,m^{\prime }}=-\frac{i}{2}N\theta [& (1-m^{\prime }){\varrho }_{m,m^{\prime }+\frac{2}{N}}+(1+m^{\prime }){\varrho }_{m,m^{\prime }-\frac{2}{N}}\\ & -(1-m){\varrho }_{m+\frac{2}{N},m^{\prime }}+(1+m){\varrho }_{m-\frac{2}{N},m^{\prime }}].\end{array}$$

$${\partial }_t{\varrho }_{m,m^{\prime }}=\{\begin{array}{l}\gamma \sum _{m^{\prime \prime }=1/N}^{1-1/N}{\mathcal{N}}_{m^{\prime \prime }}{\varrho }_{m^{\prime \prime },m^{\prime \prime }};m=m^{\prime }=1\\ \gamma \sum _{m^{\prime \prime }=-1/N}^{-1+1/N}{\mathcal{N}}_{m^{\prime \prime }}{\varrho }_{m^{\prime \prime },m^{\prime \prime }};m=m^{\prime }=-1\\ -\gamma {\varrho }_{m,m^{\prime }};\mathrm{otherwise}\end{array},$$

$${\varphi }_m(t)={\delta }_{m,m(0)}{e}^{-\gamma t}+{\delta }_{m,\pm 1}(1-e^{-\gamma t}),$$

Full text

Introduction. Understanding how statistical mechanics emerges in closed quantum many-body systems undergoing coherent dynamics with time-independent Hamiltonians has been one of the major themes of physics research over the last few decades. More recently, attention has been focussed on closed systems with time-periodic (“Floquet”) Hamiltonians, where fundamentally novel out-of-equilibrium phases describable in macroscopic terms have been discovered; none more prominent than the $\pi$-spin glass also termed the discrete time crystal (DTC) Khemani et al. (2016); Else et al. (2016); Moessner and Sondhi (2017); Sacha (2015); Yao et al. (2017); Zhang et al. (2017); Choi et al. (2017); Pal et al. (2018); Rovny et al. (2018); O’Sullivan et al. (2018).

Generically, a major obstacle to working with Floquet systems is that they suffer heat death due to unbounded increase of entropy, approaching an infinite-temperature state at long times D’Alessio and Rigol (2014); Lazarides et al. (2014); Ponte et al. (2015a). This heating can be avoided by introducing disorder-induced localisation Lazarides et al. (2015); Ponte et al. (2015b), or by coupling the system to an external environment which drains energy from the system. The former approach, used in Ref. Khemani et al. (2016), additionally endows the Floquet eigenstates with a discrete-symmetry broken spatial glassy order. Crucially, the eigenstates connected by the symmetry (and having the same spatial ordering pattern) are separated in quasienergy by $\pi/T$ with $T$ the driving period, leading to the subharmonic oscillation of an appropriate local observable. It was shown in Ref. Lazarides and Moessner (2017) that an external Markovian environment, unless explicitly fine-tuned, destroys such a delicate coherence required for the subharmonic oscillations; the system is driven towards a mixture of various Floquet eigenstates all with uncorrelated patterns of the spatial glassy order.

Here we analyse general dissipative Floquet systems in a different setting where the ordered phases, time-crystalline or otherwise, are stabilised by the dissipation, and which would be entirely absent without it. The mechanism, manifestly different from that of the $\pi$-SG, involves a periodic rotation between two “sectors” of a Hilbert space followed by dissipative “cooling” to states distinguishable by a measure such as magnetisation.

This intuitively appealing picture ignores the possibility that the dissipation in addition to the desired cooling also generates noise and eventual loss of phase coherence in the oscillations. In terms of the density matrix of the system, the noise potentially translates to the probability distribution of the observable over individual realisations being broad, which we argue destroys the DTC. We show that while this is indeed the case for 1-$d$ spin chains with local interactions and dissipation, in 2-$d$ the dissipative processes can naturally lead to a narrow distribution and hence persistent time-crystalline order unlike in 1-$d$. We note that this broadening mechanism can be absent altogether in mean-field dynamics such as for fully-connected models Russomanno et al. (2017); Gong et al. (2018); Gambetta et al. (2019); Yao et al. (2018).

How a dissipative DTC can exist$\ldots$

First, let us describe the general arguments for the stability, or lack thereof, of DTCs to dissipation, before demonstrating them concretely using specific quantum spin-models. For a visual schematic, see Fig. 1.

• •

The Hilbert space is divided into different (for concreteness: two) sectors, $\mathcal{H}_{\pm}$, which could be symmetry sectors or simply based on an empirical criterion based on the expectation value of an observable, and

• •

each of these sectors has a manifold of states $\{\ket{G_{\pm}}\}$ which posses quantum order characterised by an observable, $\hat{M}$, which may for example be due to symmetry-broken order. Crucially,

• •

the two ordered manifolds, $\{\ket{G_{\pm}}\}$, cannot be connected to each other via local operators, and

• •

the expectation value $m$ of the order parameter, $\hat{M}$, is sufficiently narrowly distributed over the states within each of the manifolds that the two distributions for the two manifolds do not overlap (up to exponentially small corrections).

In Fig. 1, the two sectors are the positive- and negative-$m$ halves of the vertical axis. Within this setting consider a two-step Floquet protocol in the presence of dissipative processes such that,

• •

In the first step (“rotation phase”), the system evolves under the simultaneous action of a Hermitian rotation operator $H_{R}$ and the dissipative terms dissipation. In the absence of dissipation the evolution over the step is unitary and given by $U_{R}(\theta)$, which maps neighbourhoods of the ground state manifold of one sector to states of the other sector and vice-versa. For a particular $\theta_{\ast}$, $U_{R}(\theta_{\ast})$ maps states from $\{\ket{G_{+}}\}$ exactly onto states from $\{\ket{G_{-}}\}$ and vice versa.

• •

In the second step (“cooling phase”), the system is governed by a Hamiltonian for which $\{\ket{G_{\pm}}\}$ are ground state manifolds as well as by the same dissipation processes as in the first step.

The dissipative processes during both steps cool the system down so that under their sole influence all states in the $\mathcal{H}_{+}$ sector would be driven to $\{\ket{G_{+}}\}$ and likewise for $\mathcal{H}_{-}$. Under the combined action of the unitary and dissipative terms, they take the system towards the $\ket{G_{\pm}}$ of the sector it instantaneously finds itself in.

A quantum system driven with such a dissipative Floquet protocol, initialised in either of the ground state manifolds, shows a time crystalline response trivially if $\theta=\theta_{\ast}$, as the expectation value of $\hat{M}$ oscillates stroboscopically between that in $\{\ket{G_{+}}\}$ and $\{\ket{G_{-}}\}$ with a period twice that of the Floquet drive, provided the dissipation is inactive in the rotation cycle. A certain robustness of the temporal order to deviations of $\theta$ from $\theta_{\ast}$ is expected to be induced via the dissipation: if the unitary rotation does not take states from the ground state manifold of one sector (say $\{\ket{G_{+}}\}$) entirely to that of the other sector ($\{\ket{G_{-}}\}$), but admixes nearby excited states in the other sector, the cooling step of the drive can push the weight back onto the $\{\ket{G_{-}}\}$ manifold. Simply put, cooling kills off the excitations left behind by the imperfect rotation, stabilising a DTC.

$\ldots$* and how it can die:* In a broad sense, the dissipative processes have three effects:

• •

Hindering rotation during rotation phase Recall that $\{\ket{G_{+}}\}$ and $\{\ket{G_{-}}\}$ are not connected via local operators. $U$ then naturally has the form of a global rotation of the degrees of freedom. If the rotation phase is not instantaneous, the state goes through excited states at intermediate times. However, the dissipation cools the system, opposing this creation of excitations hence making the rotation process less effective. Therefore the overall rotation with dissipation is less than without, trapping (part of) the weight in the wrong sector. This is unfavourable to the presence of a stable DTC.

• •

Correcting error caused by imperfect rotation during cooling phase Imperfect rotation potentially leaves the state in the correct sector, but not in the $\ket{G_{\pm}}$ manifolds; dissipation corrects this, favouring the DTC.

• •

Broadening the distribution during both phases: The rotation and the cooling acting in conjunction can, for short times, increase the width of the distribution of the state’s overlaps with the excited states such that the resulting state is spread over both the sectors. In the following cooling cycle of the Floquet drive, the weights in each sector can get pushed to their respective ground state manifolds, resulting in a finite weight in the wrong sector (see later discussion and Fig. 1). This is generally fatal to the DTC.

Of the three, the third (broadening) invariably causes the DTC signal to decay eventually. In its absence, when the probability distribution of the observable remains sharp over time, a stable DTC phase is possible with the first two mechanisms determining the parameter regime of the stability. Let us also note that while dissipation is favourable for the temporal order in the cooling cycle, it is detrimental in the rotation cycle, and it is a priori not obvious whether increasing the strength of the dissipation from some finite value favours or disfavours the temporal order.

In what follows, we introduce an explicit microscopic model and give three examples of dissipative processes. First we show that dissipation that cleanly separates the two sectors but broadens the magnetization distribution leads to a decay of the oscillation. We then introduce spatially local dissipation processes and show that: 1) in 1-$d$ they fail to separate the two sectors, cause broadening, and lead to a decaying oscillation. In 2-$d$, they may cleanly separate the sectors and in addition do not result in broadening, so that in this case a stable DTC appears.

Quantum spin systems.

To analyse the above ideas in a concrete setting, we consider a system of spins-1/2, first in 1-$d$. Using a basis constituted by the products states of $\sigma^{z}$ (which we henceforth denote as $\{\ket{\alpha}\}$), the two sectors $\mathcal{H}_{\pm}$ can be taken to be the set of product states, $\{\ket{\alpha_{\pm}}\}$ which satisfy $\braket{\alpha_{\pm}}{\hat{M}}{\alpha_{\pm}}\gtrless 0$ respectively with $\hat{M}=\sum_{\ell}\sigma^{z}_{\ell}$ 111Modulo the ambiguity for the product states with zero magnetisation; we choose to put half of them in the first sector and half of them in the other. The choice has no bearing on the subsequent dynamics.. This is a natural choice for a system described by a ferromagnetic Ising Hamiltonian

[TABLE]

as in the limit of $g\to 0$, $\{\ket{\alpha}\}$ is a possible set of eigenstates. Moreover, for $g\neq 0$, $\ket{G_{\pm}}$ are adiabatically connected to the $\ket{\Uparrow}$ (all-up) and $\ket{\Downarrow}$ (all-down) states as long as the Hamiltonian is in the ferromagnetic phase, $|g|<1$. Note that, defining the two sectors and the corresponding ground state manifolds in the fashion we do, also allows us to label the basis states and the sectors with the magnetisation density $m=\braket{\hat{M}}/N$ ($N$ being the system size).

The unitary operator $U_{R}$ which in the thermodynamic limit maps states $\{\ket{G_{+}}\}\leftrightarrow\{\ket{\alpha_{-}}\}$ is given by $U_{R}(\theta)=\exp\left[-i\theta\sum_{\ell}\sigma^{x}_{\ell}\right]$ with $\theta\in(\pi/4,\pi/2]$ and is produced by the action of the Hamiltonian $H_{R}=\frac{\theta}{t_{R}}\sum_{j}\sigma^{x}_{j}$ over time $t_{R}$. It easily follows that for $\theta_{\ast}=\pi/2$, $U_{R}(\theta_{\ast})$ precisely maps the all-up state to the all-down state. In fact, since the ground state of the system breaks the $\mathbb{Z}_{2}$ symmetry of the Hamiltonian spontaneously, $U_{R}(\theta_{\ast})$ connects the two ground states exactly throughout the ferromagnetic phase. As anticipated, the rotation $U_{R}(\theta)$ is manifestly a non-local operation. In the thermodynamic limit, a product state with definite magnetization $m$ is mapped to a (in the $\sigma^{z}$-basis, non-product) state with definite magnetization $m\cos(2\theta)$, so that the rotation operation does not result in broadening.

Depending on the dissipative processes involved, this model can show decaying (bottom of Fig. 1) or persistent (top) subharmonic oscillations depending on whether broadening occurs or not. In the following, we use three explicit Markovian dissipative processes to show how the absence (presence) of broadening due to them favours (disfavours) the persistence of the temporal order.

Lindblad dynamics.

Focussing on Markovian dissipative processes, the equation of motion for the density matrix of the system is governed by the Lindblad equation

[TABLE]

where $H(t)$ is the time-dependent (in our case, time-periodic) Hamiltonian and $\{L_{i}\}$ is the set of time-independent quantum jump operators which arise due to the coupling to the dissipative environment. Our binary Floquet protocol with period $T=t_{C}+t_{R}$ is

[TABLE]

We focus mostly on the $g\to 0$ limit of $\hat{H}_{\mathrm{TFIM}}$.

Direct jump operators. To demonstrate the deleterious effects of broadening we begin by considering a set of jump operators $\{L_{\alpha}\}$

[TABLE]

where $m_{\alpha}$ denotes the magnetisation of the product state $\ket{\alpha}$. These jump operators take the weight from any product state and transfer it directly to the ground state of the corresponding sector as well as causing exponential decay of offdiagonal elements of the density matrix in the product state basis. They therefore provide very efficient cooling. However, they lead to broadening of the distribution leading to the decay of the oscillatory signal. To show this explicitly, we study the time-dependent magnetisation of the system starting from the $\ket{\Uparrow}$ state, making use of a simplification due to translational invariance to access very large system sizes sup .

Fixing one of $\gamma t_{R}$ and $\gamma t_{C}$ and varying the other leads to the results shown in the two leftmost columns of Fig. 2: for either finite $\gamma t_{C}$, finite $\gamma t_{R}$, or both finite, the time-crystalline response of $m$ decays exponentially with $t$ while the lifetime of the subharmonic oscillations grows exponentially with $\gamma t_{C}$. On the other hand it decays at least polynomially with $\gamma t_{R}$, so overall stronger dissipation has opposite effects during each part of the driving. Nevertheless broadening of the distribution means the oscillations always decay via the mechanism shown in the lower panel of Fig. 1.

The magnetisation vanishes with time due to the state strobscopically being in a mixture of both ground state manifolds, with opposite magnetisations.

One then expects that an observable finite and equal in both the ground states will remain finite in a statistical mixture of the two, such as is obtained at long times in the present case. Such an observable is the correlator

[TABLE]

The persistent oscillations of $C(t)$ are shown in the rightmost column of Fig. 2. This is a fundamental difference between this dissipative Floquet phase and the $\pi$-spin glass, where persistent oscillations of $C$ imply those of an initially finite $m$. The amplitude of the oscillations decreases with decreasing $\gamma t_{C}$ because the signal is strongest in the two fixed point states $\ket{\Uparrow/\Downarrow}$ and $\gamma t_{C}$ controls how well the system is cooled into the two ground states. In this sense, this order is also induced by dissipation. Increasing $\gamma t_{R}$ causes $C$ to remain close to unity (its value in the ground states) and resist rotation, consistent with the earlier general arguments.

The direct jump operators demonstrate that the broadening of the distribution in magnetisation is fatal to the time-crystalline order, even when the dissipation cleanly separates the two sectors.

Domain-wall annihilating jump operators. We now introduce a set of jump operators which avoid broadening in a natural way. These operators cause dynamics that only move domain walls (DWs): a free standing DW can move but not disappear. However two DWs can move into each other and annihilate. Such dynamics are fundamentally different in the 1- and 2-dimensional cases.

Denoting the neighbours of a site $\ell$ by $\{r_{\ell}\}$, to each product state $\ket{\alpha}$ and site $\ell$ there corresponds a jump operator

[TABLE]

where $M_{\{r_{\ell}\},\alpha}$ is the net magnetisation of the spins in $\{r_{\ell}\}$ and $\mathcal{P}_{i,s_{i}^{\alpha}}$ is a projector onto the spin at site $i$ in spin-state $s_{i}^{\alpha}=\uparrow/\downarrow$. The corresponding Lindblad dynamics along the diagonal is governed by a Pauli master equation

[TABLE]

with the $\gamma$ determined according to the rules above while the off-diagonals decay exponentially.

In 1-$d$, the dynamics along the diagonal amounts to the DWs executing a random walk, i.e. diffusing. The probability distribution of magnetisation starting from a sharp value $m$ broadens at short times (Fig. 3 top left), and at long times becomes bimodal with two peaks at $\pm 1$ of height such that $-p(-1)+p(1)=m$ (Fig. 3 top right), as found by solving Eq. (7) using a classical kinetic Monte Carlo approach. The resulting destruction of the DTC is shown in the bottom two panels.

For $d\geq 2$, this type of domain wall dynamics eventually eliminates the minority phase by effectively causing a line (or surface) tension, tending to minimise the area of the interface between two non-conserved phases (see Fig. 4 for two examples of allowed transitions and Sup. Mat. sup for a demonstration of how the dynamics minimises the interface length) so that the dissipation cleanly separates the two sectors. In general, local dissipative processes lowering the energy of (ferromagnetic) Ising-type hamiltonians in $d\geq 2$ will behave in a qualitatively similar way. Less obviously, the dissipative dynamics does not broaden the magnetization distribution starting from a state sharp in magnetization, Fig. 5. It then follows that a DTC phase may be stable; this is supported by the lower panels of Fig. 5 where a rapid rotation is shown to result in persistent subharmonic oscillations while a slow rotation in a ferromagnetic phase in which the magnetization never changes sign.

In order to show this explicitly on finite-sized systems a numerical solution of the full Lindblad equation is desirable. However note that Eq.(6) dictates that the number of jump operators is exponentially large in $N$, making a numerical solution all but impossible. We surmount this by noting that one can have another set of jump operators which lead to identical (to those of Eq.(6)) dynamics for the diagonal elements of the density matrix but whose number grows polynomially with $N$.

Each ‘majority rule’ operator, for a given site, consists of a product of projectors onto the site’s neigbouring spins (the set is over all possible configurations of $\mathcal{P}^{\uparrow/\downarrow}$ for the neighbours) multiplied by the spin raising (lowering) operator for the site if the projector configuration has fewer $\mathcal{P}^{\downarrow(\uparrow)}$ compared to $\mathcal{P}^{\uparrow(\downarrow)}$, see e.g. Fig. 4. The resulting dynamics is displayed in Fig. 5

.

Conclusions and outlook. We have discussed general mechanisms leading to dissipative (de)stabilisation of DTCs in disorder-free dissipative Floquet systems. Our example of a dissipation-stabilised DTC is completely distinct, relying on a fundamentally different mechanism, from the $\pi$-spin glass introduced in Khemani et al. (2016), which in turn is unstable to dissipation Lazarides and Moessner (2017). We also uncover other, non-DTC but still dissipation-induced phases.

Phenomenologically, a crucial difference between the dissipative Floquet system and the $\pi$-SG Khemani et al. (2016) is that in the former the oscillations in the magnetisation can decay even though its correlation function, $C(t)$, can synchronise and oscillate persistently, while in the latter one implies the other.

While the mechanism of obtaining period-doubling from periodic switching between distinct sectors of Hilbert space is intuitively transparent, our analysis of its failure modes we believe also sheds light on recent work finding stable subharmonic oscillations Russomanno et al. (2017); Gong et al. (2018); Gambetta et al. (2019). In these works the system Hamiltonian is fully connected. This typically leads to a stochastic description in which the noise vanishes with diverging system size, so that the master equation for the density matrix results in no broadening over the time evolution Gelhausen and Buchhold (2018); Benatti et al. (2018).

We believe that generally, treatments for short-range models based on approximate mean-field and other analyses involving only a few effective degrees of freedom may erroneously find stable time-crystalline behaviour in the absence of such a noise suppression mechanism. Our proposal is that the role provided by long-range interactions can, however, be replaced by the effectively macroscopic rigidity of the ordered component of a symmetry-broken system such as the Ising magnet subjected to the dissipative processes discussed here.

Finally, we have only considered Markovian dissipation. An open and interesting problem is to understand whether the physics unveiled here is changed qualitatively in the non-Markovian case, and whether there are non-fine-tuned non-markovian environments that lead to interesting new examples of oscillatory dynamics in quantum systems.

I Direct jump operators

In this section of the supplementary material, we provide details of the direct jump operators [Eq. (4) of the main text] and describe how translation-invariance of the Hamiltonian as well as that of the jump operators can be exploited to simulate relatively large system sizes. Note that the action of the direct jump operators, of the form $L_{\alpha}=\sqrt{\gamma}\left[\Theta(m_{\alpha})\ket{\Uparrow}\bra{\alpha}+\Theta(-m_{\alpha})\ket{\Downarrow}\bra{\alpha}\right],$ on a basis state $\ket{\alpha}$ depends only on the magnetisation density of the state. Similarly, the action of the $\hat{H}_{\mathrm{TFIM}}$ (in the $g\to 0$ limit) as well as the global rotation operator $U_{R}(\theta)=\exp\left[-i\theta\sum_{\ell}\sigma^{x}_{\ell}\right]$ on $\ket{\alpha}$ also depends only on $m_{\alpha}$. This allows for the fact that the density matrix elements in the Ising configuration basis depend only the magnetisation densities of the basis states. Mathematically,

[TABLE]

Hence, instead of the $2^{N}$-dimensional density matrix $\rho$, it suffices for us to work with the $(2N+1)$-dimensional matrix $\varrho$. Normalisation of the density matrix in this notation is ensured via $\sum_{m}\varrho_{mm}\mathcal{N}_{m}=1$, where $\mathcal{N}_{m}=\binom{N}{N(1+m)/2}$ is the multiplicity of the product states with magnetisation density $m$. The probability distribution of the magnetisation density, $\phi_{m}$, is given by $\varrho_{mm}\mathcal{N}_{m}$ and the expectation values of the magnetisation $m$ and the correlation function $C$ are then given by

[TABLE]

respectively.

In order to solve for the dynamics of the system, we require the equation of the motion for $\varrho_{mm^{\prime}}$ under the action of the both, the unitary rotation operator as well as the jump operators. The former is described by the equation

[TABLE]

As that the jump operators transfer weight from all positive (negative) magnetisation states to the all-up (all-down) state directly with a rate $\gamma$, the evolution of $\varrho$ is described by the set of equations

[TABLE]

To see that the jump operators spread the state out in magnetisation, one can simply analyse the solutions of Eq. (S4) with the initial conditions such that the state is narrowly distributed in magnetisation around a value $m(0)$. The solutions yield that all the $\varrho_{m,m^{\prime}}$ decay exponentially with a rate $\gamma$ except for $\varrho_{1,1}$ or $\varrho_{-1,-1}$ depending on if $m(0)\gtrless 0$, which approaches unity exponentially. Specifically, for the initial condition $\varrho_{m(0)m(0)}(0)=1/\mathcal{N}_{m(0)}$, the solution to the time dependent distribution $\phi_{m}$ can be expressed as

[TABLE]

for $m(0)\gtrless 0$, which manifestly shows that the time-dependent state is not sharply distributed in $m$ and the distribution has a finite standard deviation at any finite time $t$, see Fig. S1 for results. Eqs. (S4) and (S3) together describe the evolution of the density matrix of the system in the rotation cycle, whereas the former suffices in the cooling cycle if we work in the $g\to 0$ limit

II Domain-wall annihilating jump operators

In this section, we discuss further details of the domain-wall annihilating operators. The operators in Eq. (6) (main text) locally shift domain walls such that two of them can annihilate each other upon meeting. However the set of operators are rather inconvenient from the point of view of a numerical solution of the corresponding Lindblad equation as there are exponentially (in $N$) many jump operators. We surmount this problem by considering a different set of jump operators which have the same effect on basis product states as the ones in Eq. (6) (main text) but which contains only polynomially in $N$ of them. In particular, this new set contains $2^{Z}\times N$ jump operators, ($2^{Z}$ for each spin) where $Z$ is the coordination number of the lattice and $N$ is the number of spins. As mentioned in the main text, each operator in the set consists of a spin-flip (or lowering or raising) term for a given site, say $\ell$, multiplied to projectors onto each configuration of its neighbours (hence $2^{Z}$ of them). If a neighbour-configuration has a net magnetisation which is positive (negative), the corresponding projector is multiplied to $\sigma^{+}_{\ell}$ ($\sigma^{-}_{\ell}$) and if the neighbour-configuration has equal number of up and down spins, the projector is multiplied to $\sigma^{x}_{\ell}$, the latter with much smaller rate.

For a square lattice in 2-$d$, for a given site $\ell$ at position $(x,y)$, there are three classes of jump operators:

• •

All four of the neighbours are aligned with each other but anti-aligned with the spin at $\ell$. There are two jump operators for this case, $\mathcal{P}_{x+1,y+1}^{\uparrow}\mathcal{P}_{x+1,y-1}^{\uparrow}\mathcal{P}_{x-1,y+1}^{\uparrow}\mathcal{P}_{x-1,y-1}^{\uparrow}\sigma^{+}_{x,y}$ and the same with $\uparrow\rightarrow\downarrow$ and $+\rightarrow-$. See Fig. 4 (right) [main text] for a visual example.

• •

Three of the neighbours are anti-aligned with the spin at $\ell$ and one is aligned. There are four configurations each for a the net magnetisation of the neighbours being positive or negative. So there are a total of eight operators in this class. As example, one of the operators in this class is $\mathcal{P}_{x+1,y+1}^{\uparrow}\mathcal{P}_{x+1,y-1}^{\uparrow}\mathcal{P}_{x-1,y+1}^{\uparrow}\mathcal{P}_{x-1,y-1}^{\downarrow}\sigma^{+}_{x,y}$, and the projector on the down spin could be on any of the four neighbours which generates the other three operators. Similarly $\uparrow\rightarrow\downarrow$ and $+\rightarrow-$ generates the other four operators in this class.

• •

Two of the neighbours are up and two of them are down. In this case, we always flip the spin at site $\ell$, but the operator has a much smaller rate. An example operator in this class is $\mathcal{P}_{x+1,y+1}^{\uparrow}\mathcal{P}_{x+1,y-1}^{\uparrow}\mathcal{P}_{x-1,y+1}^{\downarrow}\mathcal{P}_{x-1,y-1}^{\downarrow}\sigma^{x}_{x,y}$ There are six operators in this class (one for each of the two-up two-down configurations).

Note that the first two classes of the jump operators tries to reduce the length of the domain wall and favour the majority phase. However, they can get stuck if they encounter a straight domain wall. The third class of the jump operators serve to unfreeze such potentially frozen domain walls.

The effect of these jump operators on the diagonal elements of the density matrix can be understood from a classical Monte Carlo simulation using Glauber dynamics with the local energy cost function for a spin at site $\ell$ given by $\sigma^{z}_{\ell}\sum_{r\in\{r_{\ell}\}}\sigma^{z}_{r}$. While the Monte Carlo at zero temperature would try to generically take the system to the majority phase, there is a technical subtlety. The zero temperature classical Monte Carlo which tries to reduce the string tension of the domain wall can get stuck if it encounters a domain wall which straight; hence we run the Monte Carlo at a finite but small temperature. This essentially is the manifestation of the third class of the jump operators discussed above.

Fig. S2 shows a specific Monte Carlo trajectory with each panel showing the spin-configuration on the lattice at specific Monte Carlo times; green denotes regions of down-spin and yellow, up-spin. The evolution of the configuration clearly shows the shrinking of the domain-wall lengths and corresponds to a macroscopic coarse-grained version of the processes shown in Fig. 4 (main text).

As a final remark, we mention that the full solution of the Lindblad equation in the 2-$d$ case (results of Fig. 5 (main text)) was obtained using the Monte Carlo wave-function method, see [Mølmer et al., J. Opt. Soc. Am. B 10, 524 (1993)] for details.