Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Spatiotemporal patterns in the active cyclic Potts model

Hiroshi Noguchi noguchi@issp.u-tokyo.ac.jp Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan    Jean-Baptiste Fournier jean-baptiste.fournier@u-paris.fr Laboratoire Matière et Systèmes Complexes (MSC), Université Paris Cité & CNRS, 75013 Paris, France
Abstract

The nonequilibrium dynamics of a cycling three state Potts model is studied on a square lattice using Monte Carlo simulations and continuum theory. This model is relevant to chemical reactions on a catalytic surface and to molecular transport across a membrane. Several characteristic modes are formed depending on the flipping energies between successive states and the contact energies between neighboring sites. Under cyclic symmetry conditions, cycling homogeneous phases and spiral waves form at low and high flipping energies, respectively. In the intermediate flipping energy regime, these two modes coexist temporally in small systems and/or at low contact energies. Under asymmetric conditions, we observed small biphasic domains exhibiting amoeba-like locomotion and temporal coexistence of spiral waves and a dominant non-cyclic one-state phase. An increase in the flipping energy between two successive states, say state 0 and state 1, while keeping the other flipping energies constant, induces the formation of the third phase (state 2), owing to the suppression of the nucleation of state 0 domains. Under asymmetric conditions regarding the contact energies, two different modes can appear depending on the initial state, due to a hysteresis phenomenon.

I Introduction

Spatiotemporal patterns are widely observed in nonequilibrium conditions, like waves of target and spiral shapes in two-dimensional (2D) systems Murray (2003); Mikhailov and Showalter (2006); Mikhailov and Ertl (2009); Beta and Kruse (2017); Bailles et al. (2022); Okuzono and Ohta (2003); Qu (2006); Sugimura and Kori (2015). Cyclic dynamics of three or more states are one of the conditions to generate such waves. The cyclic states are observed in various systems, such as chemical reactions, gene expression, and ecological systems Mikhailov and Showalter (2006); Mikhailov and Ertl (2009); Kerr et al. (2002); Kelsic et al. (2015). These dynamics are often explained using deterministic continuum equations Mikhailov and Showalter (2006); Mikhailov and Ertl (2009); Szolnoki et al. (2014); Reichenbach et al. (2007, 2008) and lattice Lotka–Volterra models (also called rock–paper–scissors models) Kerr et al. (2002); Kelsic et al. (2015); Szolnoki et al. (2014); Reichenbach et al. (2007, 2008); Itoh and Tainaka (1994); Tainaka (1994); Szabó et al. (1999); Szabó and Szolnoki (2002); Johnson and Seinen (2002); Szczesny et al. (2013); Juul et al. (2013); Mir et al. (2022).

The spatiotemporal patterns under noise and fluctuations are not yet fully understood. Noise is typically used to trigger an excitable wave and to generate stochastic resonance Gammaitoni et al. (1998); McDonnell and Abbott (2009). The wave propagation of subexitable media of photosensitive Belousov–Zhabotinsky reaction is enhanced by the random variation of light intensity Mikhailov and Showalter (2006); Kádár et al. (1998); Wang et al. (1999); Alonso et al. (2001). In partial differential equations, random noises are added to include mesoscale fluctuations Hildebrand and Mikhailov (1996); Reichenbach et al. (2007, 2008). In the lattice Lotka–Volterra models Kerr et al. (2002); Kelsic et al. (2015); Szolnoki et al. (2014); Reichenbach et al. (2007, 2008); Itoh and Tainaka (1994); Tainaka (1994); Szabó et al. (1999); Szabó and Szolnoki (2002); Johnson and Seinen (2002); Szczesny et al. (2013); Juul et al. (2013); Mir et al. (2022) for the predator–prey systems, noise is included as a random selection of lattice sites. Since populations increase by self-reproduction, the nucleation of species domains is not considered. In contrast, in small molecular systems, thermal fluctuations have dominant effects, and nucleation and growth are the crucial kinetic processes. Spiral and target wave patterns have been observed in chemical reactions on noble metal surfaces Ertl (2008); Brär et al. (1994); Gorodetskii et al. (1994); Barroo et al. (2020); Zeininger and et al. (2022). In a previous paper Noguchi et al. , we studied a nonequilibrium three-state Potts model under cyclic symmetry (the three states being equivalent), and reported that nucleation and growth can alter spatiotemporal patterns. We found two modes: (i) a homogeneous cycling mode (HC), where each state dominantly covers the system while changing cyclically through nucleation and growth. (ii) a spiral wave (SW) mode where the three states coexist while spiral waves are formed from the contact points of three states. These two modes can temporally coexist in small systems but not in larger systems. In this study, we examine the dynamics of the three-state active Potts model in asymmetric cycling conditions. We will show several modes and hysteresis and amoeba-like motions of small biphasic domains. A non-cyclic one-state dominant phase and temporal coexistence with the SW mode newly appear.

The cyclic Potts model and method are described in Sec. II. It is a 2D lattice model for chemical reactions on a catalytic surface and molecular transport through a membrane. The dynamics under cyclic-symmetry conditions are briefly described in Sec. III. The dynamics in asymmetric flipping energy and asymmetric contact energies are described in Secs. IV and V, respectively. The theoretical analysis by a continuum theory is described in Sec. VI. Finally, a summary is presented in Sec. VII.

Refer to caption
Figure 1: (a) Three-state Potts model in a square lattice. (b) Nonequilibrium cyclic flipping of the three states driven by h01+h12+h20=hcyc0subscript01subscript12subscript20subscriptcyc0h_{01}+h_{12}+h_{20}=h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0. (c) Possible realizations with energy levels depicted: (1) Chemical reaction on a catalytic surface, (2) Molecular transport through a membrane, (3) Three-state process with 11+superscript1superscript11^{-}\to 1^{+}1 start_POSTSUPERSCRIPT - end_POSTSUPERSCRIPT → 1 start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT excitation.

II Active Cyclic Potts Model

We consider a 2D square lattice with sites i𝑖iitalic_i having three states si=αsubscript𝑠𝑖𝛼s_{i}=\alphaitalic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_α, with α{0,1,2}𝛼012\alpha\in\{0,1,2\}italic_α ∈ { 0 , 1 , 2 }, as shown in Fig. 1(a). Two nearest neighbor sites have a contact energy Jsisjsubscript𝐽subscript𝑠𝑖subscript𝑠𝑗-J_{s_{i}s_{j}}- italic_J start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT. In addition, the three states have different self-energies εαsubscript𝜀𝛼\varepsilon_{\alpha}italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT, so that the total Hamiltonian reads

H=Hint+iεsi,Hint=ijJsisj.formulae-sequence𝐻subscript𝐻intsubscript𝑖subscript𝜀subscript𝑠𝑖subscript𝐻intsubscriptdelimited-⟨⟩𝑖𝑗subscript𝐽subscript𝑠𝑖subscript𝑠𝑗H=H_{\mathrm{int}}+\sum_{i}\varepsilon_{s_{i}},\quad H_{\mathrm{int}}=-\sum_{% \langle ij\rangle}J_{s_{i}s_{j}}.italic_H = italic_H start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT + ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_ε start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT , italic_H start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT = - ∑ start_POSTSUBSCRIPT ⟨ italic_i italic_j ⟩ end_POSTSUBSCRIPT italic_J start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT . (1)

To model the attraction between like states, we mainly use symmetric contact energies Jαα=Jsubscript𝐽𝛼𝛼𝐽J_{\alpha\alpha}=Jitalic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = italic_J, and set Jαβ=0subscript𝐽𝛼𝛽0J_{\alpha\beta}=0italic_J start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = 0 for αβ𝛼𝛽\alpha\neq\betaitalic_α ≠ italic_β. This corresponds to the standard three-state Potts model with external fields Potts (1952); Binder (1981). We define the equilibrium flipping energies hαβsubscript𝛼𝛽h_{\alpha\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT as the variations hαβ=εαεβsubscript𝛼𝛽subscript𝜀𝛼subscript𝜀𝛽h_{\alpha\beta}=\varepsilon_{\alpha}-\varepsilon_{\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT. Hence, h01+h12+h20=0subscript01subscript12subscript200h_{01}+h_{12}+h_{20}=0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 0 by definition in thermal equilibrium conditions.

We extended this model to the nonequilibrium situation in which hαβεαεβsubscript𝛼𝛽subscript𝜀𝛼subscript𝜀𝛽h_{\alpha\beta}\neq\varepsilon_{\alpha}-\varepsilon_{\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ≠ italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT, but where still hβα=hαβsubscript𝛽𝛼subscript𝛼𝛽h_{\beta\alpha}=-h_{\alpha\beta}italic_h start_POSTSUBSCRIPT italic_β italic_α end_POSTSUBSCRIPT = - italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT. Then we have h01+h12+h20=hcyc0subscript01subscript12subscript20subscriptcyc0h_{01}+h_{12}+h_{20}=h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0 (see Fig. 1(b)). A possible way to do this is to add a nonequilibrium or active contribution to some or all of the flipping energies, in the form hαβ=εαεβ+hαβneqsubscript𝛼𝛽subscript𝜀𝛼subscript𝜀𝛽subscriptsuperscriptneq𝛼𝛽h_{\alpha\beta}=\varepsilon_{\alpha}-\varepsilon_{\beta}+h^{\mathrm{neq}}_{% \alpha\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT + italic_h start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT. Hence, the detailed-balance condition is locally satisfied for one flip between states s=α𝑠𝛼s=\alphaitalic_s = italic_α and β𝛽\betaitalic_β but not globally for cycles (s=0120𝑠0120s=0\to 1\to 2\to 0italic_s = 0 → 1 → 2 → 0).

Three possible designs are described below: chemical reaction, molecular transport, and excitation, as shown in Fig. 1(c). (1) Reaction on a catalytic surface: molecules bind and unbind to the surface, with the state s=0𝑠0s=0italic_s = 0 corresponding to an empty surface site, and the other two states to a site occupied by a molecule, either in the form s=1𝑠1s=1italic_s = 1 or s=2𝑠2s=2italic_s = 2. The surface catalyzes the reaction from s=1𝑠1s=1italic_s = 1 to s=2𝑠2s=2italic_s = 2, whereas this reaction is kinetically frozen in the bulk. (2) Molecular transport between two sides of a membrane: amphiphilic molecules bind to both sides of a bilayer membrane (states s=1𝑠1s=1italic_s = 1 and 2222) and switch between them (flip–flop). The molecules are transported by a difference in chemical potential between the solutions on the two sides of the membrane Miele et al. (2020); Holló et al. (2021); Noguchi (2023). (3) An excitation process: s=1𝑠1s=1italic_s = 1 has a low-energy state s=1𝑠superscript1s=1^{-}italic_s = 1 start_POSTSUPERSCRIPT - end_POSTSUPERSCRIPT and a high-energy state s=1+𝑠superscript1s=1^{+}italic_s = 1 start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT triggered by external means (e.g., photoexcitation or ATP hydrolysis). Then, if the backward transformation from s=1+𝑠superscript1s=1^{+}italic_s = 1 start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT to s=0𝑠0s=0italic_s = 0 is negligibly slow, we have effectively hcyc=ε1+ε1>0subscriptcycsubscript𝜀superscript1subscript𝜀superscript10h_{\mathrm{cyc}}=\varepsilon_{1^{+}}-\varepsilon_{1^{-}}>0italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT 1 start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT 1 start_POSTSUPERSCRIPT - end_POSTSUPERSCRIPT end_POSTSUBSCRIPT > 0. Experimentally, chemical waves have been observed at catalytic surfaces (H2 or CO oxidation and NO or NO2 reduction on noble metal surfaces such as palladium and ruthenium) Ertl (2008); Brär et al. (1994); Gorodetskii et al. (1994); Barroo et al. (2020); Zeininger and et al. (2022). The transfer of water molecules through a chiral liquid-crystalline monolayer induces a target wave Tabe and Yokoyama (2003). Although more than two reactions or complicated molecular interactions occur in these experimental systems, we constructed a simple minimum model to capture the essence.

For simplicity, no state exchange between neighboring sites (diffusion) is considered. The thermal energy and the lattice spacing are normalized to unity. We mainly use J00=J11=J22=J=2subscript𝐽00subscript𝐽11subscript𝐽22𝐽2J_{00}=J_{11}=J_{22}=J=2italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = italic_J = 2 to induce a phase separation between different states as a symmetric condition of the contact energy. To examine the effects of asymmetric contact energies, J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT is varied from 1.751.751.751.75 to 2.82.82.82.8 with keeping J11=J22=2subscript𝐽11subscript𝐽222J_{11}=J_{22}=2italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = 2, in Sec. V. Our lattice has N𝑁Nitalic_N sites with a side length of N𝑁\sqrt{N}square-root start_ARG italic_N end_ARG under periodic boundary conditions. The state of a site is flipped according to a Monte Carlo (MC) method. A site is chosen at random, and then it is flipped to either of the other two states with 1/2121/21 / 2 probability. The new state is accepted or rejected with Metropolis probability

psisi=min(1,eΔHsisi)subscript𝑝subscript𝑠𝑖subscriptsuperscript𝑠𝑖1superscript𝑒Δsubscript𝐻subscript𝑠𝑖subscriptsuperscript𝑠𝑖p_{s_{i}s^{\prime}_{i}}=\min\left(1,e^{-\Delta H_{s_{i}s^{\prime}_{i}}}\right)italic_p start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT = roman_min ( 1 , italic_e start_POSTSUPERSCRIPT - roman_Δ italic_H start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT end_POSTSUPERSCRIPT ) (2)

in which ΔHsisi=HintHinthsisiΔsubscript𝐻subscript𝑠𝑖subscriptsuperscript𝑠𝑖subscriptsuperscript𝐻intsubscript𝐻intsubscriptsubscript𝑠𝑖subscriptsuperscript𝑠𝑖\Delta H_{s_{i}s^{\prime}_{i}}=H^{\prime}_{\mathrm{int}}-H_{\mathrm{int}}-h_{s% _{i}s^{\prime}_{i}}roman_Δ italic_H start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT = italic_H start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT - italic_H start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_s start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT is the energy variation in the change from the old state to the new one. This procedure is performed N𝑁Nitalic_N times per MC step (time unit). Statistical errors are calculated from three or more independent runs.

Refer to caption
Figure 2: Cyclic flipping in symmetric condition (h01=h12=h20=hsubscript01subscript12subscript20h_{01}=h_{12}=h_{20}=hitalic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h and Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2). (a) System of size N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Probability pphase=pthreesubscript𝑝phasesubscript𝑝threep_{\mathrm{phase}}=p_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_phase end_POSTSUBSCRIPT = italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT of the three-phase state (Nα>0.02Nsubscript𝑁𝛼0.02𝑁N_{\alpha}>0.02Nitalic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT > 0.02 italic_N, N1>0.02Nsubscript𝑁10.02𝑁N_{1}>0.02Nitalic_N start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT > 0.02 italic_N, and N2>0.02Nsubscript𝑁20.02𝑁N_{2}>0.02Nitalic_N start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT > 0.02 italic_N). At h0.820.82h\leq 0.82italic_h ≤ 0.82, depending on the initial condition, both the SW and the HC modes (indicated by the vanishing of pthreesubscript𝑝threep_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT) can appear. At h0.830.83h\geq 0.83italic_h ≥ 0.83, only the SW mode can exist. (b) System of size N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. Probability pthreesubscript𝑝threep_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT of the three-phase state and probability ponesubscript𝑝onep_{\mathrm{one}}italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT of one-phase state (Ns>0.98Nsubscript𝑁𝑠0.98𝑁N_{s}>0.98Nitalic_N start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT > 0.98 italic_N with s{0,1,2}𝑠012s\in\{0,1,2\}italic_s ∈ { 0 , 1 , 2 }). The HC and SW modes temporally coexist for 1h1.1511.151\leq h\leq 1.151 ≤ italic_h ≤ 1.15. (c) Average difference flow qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT between forward and backward flips, showing that detailed balance is more violated in the SW mode than in the HC mode. The solid and dashed lines represent the data at N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT and 1282superscript1282128^{2}128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, respectively. (d) Snapshots at h=0.90.9h=0.9italic_h = 0.9 and 1.61.61.61.6 from bottm to top for N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. The light gray, blue (dark gray), and red (medium gray) represent s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively.
Refer to caption
Figure 3: Behavior of the system for low contact energy, Jαα=1.2subscript𝐽𝛼𝛼1.2J_{\alpha\alpha}=1.2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 1.2, in symmetric condition (h01=h12=h20=hsubscript01subscript12subscript20h_{01}=h_{12}=h_{20}=hitalic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h) at N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Sequential snapshots at h=0.250.25h=0.25italic_h = 0.25. (b) Probability pthreesubscript𝑝threep_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT of the three-phase state and probability ponesubscript𝑝onep_{\mathrm{one}}italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT of one-phase state (Ns>0.95Nsubscript𝑁𝑠0.95𝑁N_{s}>0.95Nitalic_N start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT > 0.95 italic_N with s{0,1,2}𝑠012s\in\{0,1,2\}italic_s ∈ { 0 , 1 , 2 }). The HC and SW modes temporally coexist for 0.17h0.20.170.20.17\leq h\leq 0.20.17 ≤ italic_h ≤ 0.2.

III Dynamics in Cyclic-Symmetry Conditions

First, we briefly describe the dynamics of the active cyclic Potts model in symmetric conditions, i.e., for h01=h12=h20=hsubscript01subscript12subscript20h_{01}=h_{12}=h_{20}=hitalic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h and J00=J11=J22=Jsubscript𝐽00subscript𝐽11subscript𝐽22𝐽J_{00}=J_{11}=J_{22}=Jitalic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = italic_J. Detailed results at J=2𝐽2J=2italic_J = 2 are given in Ref. Noguchi et al. . For large systems (N1922𝑁superscript1922N\geq 192^{2}italic_N ≥ 192 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT), at low hhitalic_h, either the HC or SW mode appears (see their definitions in Sec. I), depending on the initial state, while only the SW mode appears at high hhitalic_h. Transitions from the HC to SW modes were observed, but the reverse transition never occurred (within accessible simulation times). In the SW mode, spiral waves propagate from the contact points of three states, and the number of each state fluctuates around N/3𝑁3N/3italic_N / 3 (see the snapshots in Fig. 2(d)). For small systems (N1542𝑁superscript1542N\leq 154^{2}italic_N ≤ 154 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT), the HC mode appears at low hhitalic_h, the SW mode appears at high hhitalic_h, and the two modes temporally coexist at medium hhitalic_h. Representing temporal coexistence by a ‘+’, we denote this mixed state by SW+HC.

Refer to caption
Figure 4: Behavior of the system for weakly asymmetric flipping energies. Temporal coexistence of the SW and HC modes (SW+HC) for h01=1.1subscript011.1h_{01}=1.1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.1, h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2, and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Sequential snapshots. (b) Time evolution of the fraction of sites in each state (corresponding to snapshots in (a)). The light gray, blue (dark gray), and red (medium gray) colors represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively.

The ratios of these two modes can be quantified using the number densities Nα/Nsubscript𝑁𝛼𝑁N_{\alpha}/Nitalic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N of the three states. For Nα>0.98Nsubscript𝑁𝛼0.98𝑁N_{\alpha}>0.98Nitalic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT > 0.98 italic_N, we consider that the lattice is dominantly covered by a unique state s=α𝑠𝛼s=\alphaitalic_s = italic_α, and we call this global state ‘one’ phase (except for the condition of J=1.2𝐽1.2J=1.2italic_J = 1.2). The fraction of time spent in this state is denoted by pphase=ponesubscript𝑝phasesubscript𝑝onep_{\mathrm{phase}}=p_{\mathrm{one}}italic_p start_POSTSUBSCRIPT roman_phase end_POSTSUBSCRIPT = italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT (Fig. 2). In the HC mode, this ‘one’ phase global state cyclically changes from s=α𝑠𝛼s=\alphaitalic_s = italic_α to [α+1]delimited-[]𝛼1[\alpha+1][ italic_α + 1 ], where [α+1]=(α+1)mod3delimited-[]𝛼1modulo𝛼13[\alpha+1]=(\alpha+1)\bmod 3[ italic_α + 1 ] = ( italic_α + 1 ) roman_mod 3. Since the transient dynamics of nucleation and growth is rapid, we identify ponesubscript𝑝onep_{\mathrm{one}}italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT with the fraction of time spent in the HC mode in the cyclic-symmetry condition. Moreover, we consider that the system is in the ‘three’ phases global state when N0>0.02Nsubscript𝑁00.02𝑁N_{0}>0.02Nitalic_N start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT > 0.02 italic_N, N1>0.02Nsubscript𝑁10.02𝑁N_{1}>0.02Nitalic_N start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT > 0.02 italic_N, and N2>0.02Nsubscript𝑁20.02𝑁N_{2}>0.02Nitalic_N start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT > 0.02 italic_N. The fraction of time spent in this state is called pphase=pthreesubscript𝑝phasesubscript𝑝threep_{\mathrm{phase}}=p_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_phase end_POSTSUBSCRIPT = italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT (Fig. 2). The ratio pone/pthreesubscript𝑝onesubscript𝑝threep_{\mathrm{one}}/p_{\mathrm{three}}italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT / italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT corresponds to the ratio of the times spent in the HC and SW modes in the cyclic-symmetry condition (see Figs. 2(a) and (b)).

Refer to caption
Figure 5: Behavior of the SW mode in the case of a large asymmetry between the flipping energies, for h01=1.6subscript011.6h_{01}=1.6italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.6, h12=1.7subscript121.7h_{12}=1.7italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 1.7, h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2, and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Sequential snapshots. The interfaces between the s=1𝑠1s=1italic_s = 1 and s=2𝑠2s=2italic_s = 2 domains, which move towards the s=1𝑠1s=1italic_s = 1 domain, are faster than the other boundaries. (b) Time evolution of the fractions of sites in each state (corresponding to snapshots in (a)). The light gray, blue (dark gray), and red (medium gray) colors represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively.

Since the nucleation barrier decreases with increasing hhitalic_h (at fixed contact energy J𝐽Jitalic_J) and the nucleation more frequently occurs in larger systems, the mean lifetime of the ‘one phase’ state was found to be roughly proportional to exp(h)/N𝑁\exp(-h)/Nroman_exp ( - italic_h ) / italic_N Noguchi et al. . Furthermore, the domains were found to grow with a velocity of vwave0.07h+0.009similar-to-or-equalssubscript𝑣wave0.070.009v_{\mathrm{wave}}\simeq 0.07h+0.009italic_v start_POSTSUBSCRIPT roman_wave end_POSTSUBSCRIPT ≃ 0.07 italic_h + 0.009 at h0.50.5h\geq 0.5italic_h ≥ 0.5. Hence, at high hhitalic_h and/or large N𝑁Nitalic_N, during the growth of a s=[α+1]𝑠delimited-[]𝛼1s=[\alpha+1]italic_s = [ italic_α + 1 ] domain within a s=α𝑠𝛼s=\alphaitalic_s = italic_α domain, the nucleation of a smaller s=[α+2]𝑠delimited-[]𝛼2s=[\alpha+2]italic_s = [ italic_α + 2 ] domain often starts. Because domain growth is a stochastic process, three-state contacts are thus frequently produced. When a three-state contact appears, a spiral wave forms, as the three domain boundaries associated with this contact point exhibit a rotative motion (each s=α𝑠𝛼s=\alphaitalic_s = italic_α phase is invaded by the adjacent s=[α+1]𝑠delimited-[]𝛼1s=[\alpha+1]italic_s = [ italic_α + 1 ] phase). Conversely, the SW mode changes into the HC mode, when the three-state contact points stochastically disappear. This disappearance occurs in small systems but not in large systems (N1922𝑁superscript1922N\geq 192^{2}italic_N ≥ 192 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT). Similar size dependence of spiral waves was reported in excitable media Qu (2006); Sugimura and Kori (2015). Since the density of the three-state contacts linearly increases with increasing hhitalic_h (compare two snapshots in Fig. 2(d)), the SW mode more often changes into the HC mode at smaller hhitalic_h. More detail is given in Ref. Noguchi et al. .

During the cyclic flipping, the forward flip (s=α[α+1]𝑠𝛼delimited-[]𝛼1s=\alpha\to[\alpha+1]italic_s = italic_α → [ italic_α + 1 ]) more frequently occurs than the backward flip (s=[α+1]α𝑠delimited-[]𝛼1𝛼s=[\alpha+1]\to\alphaitalic_s = [ italic_α + 1 ] → italic_α). This is quantified by the flow qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT defined as the average difference between forward and backward flips per site, during one MC step. Figure 2(c) shows that qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT is much higher in the SW mode than in the HC mode and increases with increasing hhitalic_h in both modes.

As the contact energy J𝐽Jitalic_J decreases, the SW mode appears at lower hhitalic_h, and the cyclic change of dominant phases in the HC mode becomes faster. Thus, the temporal coexistence of these two modes can be observed in large systems. For N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, the coexistence is obtained at J=1.2𝐽1.2J=1.2italic_J = 1.2, as shown in Fig. 3. Note that the threshold Nα/N=0.95subscript𝑁𝛼𝑁0.95N_{\alpha}/N=0.95italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N = 0.95 is used for the one-phase state since the domains contain 3% of the other states at J=1.2𝐽1.2J=1.2italic_J = 1.2 (see Fig. 3(a)). The transition between the HC and SW modes occurs at higher hhitalic_h with increasing J𝐽Jitalic_J. For J=1.5𝐽1.5J=1.5italic_J = 1.5, the discontinuous transition is obtained at h3.5similar-to-or-equals3.5h\simeq 3.5italic_h ≃ 3.5 and the two modes maintain for t108similar-to𝑡superscript108t\sim 10^{8}italic_t ∼ 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT around the transition point owing to hysteresis. However, this transition likely becomes continuous for simulations of much longer periods.

Refer to caption
Figure 6: Dynamic phase diagrams in the case of asymmetric flipping energies for Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2 and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Diagram as a function of h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT and h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT for h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1. (b) Diagram as a function h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT and h12=h20subscript12subscript20h_{12}=h_{20}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT. The squares (light blue) and circles (gray) represent the HC and SW modes, respectively. The up-pointing and down-pointing triangles (green) represent the homogeneous phases s=0𝑠0s=0italic_s = 0 (denoted E0) and s=2𝑠2s=2italic_s = 2 (denoted E2), respectively. The pluses (blue), diamonds (magenta), and crosses (red) represent the temporal coexistence SW+HC, SW+E0, and SW+E2, respectively. The diagonal in (b) corresponds to the cyclic-symmetry condition (h01=h12=h20subscript01subscript12subscript20h_{01}=h_{12}=h_{20}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT).
Refer to caption
Figure 7: Probabilities to observe given fractions Ns/Nsubscript𝑁𝑠𝑁N_{s}/Nitalic_N start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT / italic_N of each state in the case of asymmetric flipping energies for Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2 and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) State E2 (homogeneous s=2𝑠2s=2italic_s = 2 phase) for h01=1.7subscript011.7h_{01}=1.7italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.7 and h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1. (b) State E2 (close to the HC mode) for h01=1.45subscript011.45h_{01}=1.45italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.45 and h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1. (c) Temporal coexistence of the SW and HC modes (SW+HC) for h01=1.1subscript011.1h_{01}=1.1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.1 and h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, corresponding to Fig. 4. (d) SW mode at h01=1.6subscript011.6h_{01}=1.6italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.6, h12=1.7subscript121.7h_{12}=1.7italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 1.7, h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, corresponding to Fig. 5.
Refer to caption
Figure 8: Density of states and mode lifetimes in the case of asymmetric flipping energies as a function of h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT for h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2, and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. The dynamic mode in which the system sits is indicated at the top of the figure. (a) Fraction of sites in each state. The fully symmetric case corresponds to the crossing of the three lines at h01=1subscript011h_{01}=1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1. For h01<1subscript011h_{01}<1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT < 1 the s=0𝑠0s=0italic_s = 0 state dominates even in the HC mode, since the system spends more time in this state. The gray up-pointing triangles, blue down-pointing triangles, and red circles represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively. (b) Probabilities of one-phase state and three-phase coexistence. (c) Mean lifetimes of the homogeneous phases (solid lines) and SW mode (dashed line).
Refer to caption
Figure 9: Dependence of the system’s properties on the asymmetry h01h12subscript01subscript12h_{01}-h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT when keeping h01+h12subscript01subscript12h_{01}+h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT constant, for Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2 and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a), (b) h01+h12=3subscript01subscript123h_{01}+h_{12}=3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 3 and h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1. (a) Fraction of sites in each state. The gray up-pointing triangles, blue down-pointing triangles, and red circles represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively. (b) Probabilities of the one-phase state and the three-phase coexistence. The corresponding dynamic modes are described at the top of the figure. (c) Average difference flow qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT between forward and backward flips. The solid lines represent the data for h01+h12=3subscript01subscript123h_{01}+h_{12}=3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 3, 3.63.63.63.6, and 4444 at h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1. The dashed lines represent the data for h01+h12=4subscript01subscript124h_{01}+h_{12}=4italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 4 and h20=1.5subscript201.5h_{20}=1.5italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1.5.
Refer to caption
Figure 10: Density of states and mode lifetimes in the case of asymmetric flipping energies as a function of h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT for h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2, and N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Fraction of sites in each state. The gray up-pointing triangles, blue down-pointing triangles, and red circles represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively. (b) Probabilities of one-phase state and three-phase coexistence. (c) Mean lifetimes of the homogeneous phases (solid lines) and SW mode (dashed line). The dynamic modes are described at the top of the figure. (d) Sequential snapshots at h01=1.3subscript011.3h_{01}=1.3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.3.

IV Dynamics for Asymmetric Flipping Energies

In this section, we describe the dynamics with asymmetric flipping energies, for symmetric contact energies Jαα=2subscript𝐽𝛼𝛼2J_{\alpha\alpha}=2italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 2. By asymmetric flipping, we mean that h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT, and h20subscript20h_{20}italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT are not all equal. We still assume, however, that h01+h12+h20=hcyc0subscript01subscript12subscript20subscriptcyc0h_{01}+h_{12}+h_{20}=h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0 and keep hβα=hαβsubscript𝛽𝛼subscript𝛼𝛽h_{\beta\alpha}=-h_{\alpha\beta}italic_h start_POSTSUBSCRIPT italic_β italic_α end_POSTSUBSCRIPT = - italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT. First, we consider a small system of N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT sites, in which the HC and SW modes can coexist in the symmetric condition. When the flipping energies slightly deviate from the symmetric condition, the dynamics exhibit only small changes. For example, the SW+HC mode is obtained at h01=1.1subscript011.1h_{01}=1.1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.1, h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, similar to that at h01=h12=h20=1subscript01subscript12subscript201h_{01}=h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, and the mean densities of s=0𝑠0s=0italic_s = 0 and s=2𝑠2s=2italic_s = 2 are slightly lower and higher 1/3131/31 / 3 in the SW mode, respectively (see Fig. 4 and Movies S1 and S2).

However, large asymmetry alters the dynamics qualitatively. In symmetric conditions, the three-phase contact points move isotropically, and the boundaries between domains move all at the same speed in the direction from s=[α+1]𝑠delimited-[]𝛼1s=[\alpha+1]italic_s = [ italic_α + 1 ] to s=α𝑠𝛼s=\alphaitalic_s = italic_α. In the case of asymmetric flipping energies, the speeds of the different boundaries are not the same. In the SW mode, when the flipping energies differ substantially, the three-phase contact points move ballistically rather than diffusively (see Fig. 5 and Movie S3). For h01=1.6subscript011.6h_{01}=1.6italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.6, h12=1.7subscript121.7h_{12}=1.7italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 1.7, h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, as in Fig. 5, small biphasic domains of s=1𝑠1s=1italic_s = 1 and 2222 move like amoeba locomotion and bacteria colony growth in the direction of the s=1𝑠1s=1italic_s = 1 region. Their large fluctuations often result in domain division and disappearance.

In the asymmetric case, the average fraction ρα=Nα/Nsubscript𝜌𝛼subscript𝑁𝛼𝑁\rho_{\alpha}=N_{\alpha}/Nitalic_ρ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT = italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N of sites in each state differs from 1/3131/31 / 3, as can be seen in Fig. 5. The contact energies turn out to be essential in determining those fractions. Indeed, for h01=1.6subscript011.6h_{01}=1.6italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 1.6, h12=1.7subscript121.7h_{12}=1.7italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 1.7, and h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, the mean-field approach developed in Sec. VI.1, which neglects both the contact energies and the spatial structures, predicts ρ0ρ10.32similar-to-or-equalssubscript𝜌0subscript𝜌1similar-to-or-equals0.32\rho_{0}\simeq\rho_{1}\simeq 0.32italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ≃ italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ≃ 0.32 and ρ20.37similar-to-or-equalssubscript𝜌20.37\rho_{2}\simeq 0.37italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ≃ 0.37, whereas for J=2𝐽2J=2italic_J = 2 the actual ratio ρ0/ρ15subscript𝜌0subscript𝜌15\rho_{0}/\rho_{1}\approx 5italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ≈ 5 (see Fig. 5). This is expected, since a single site flip within a domain involves the loss of four Jαα=Jsubscript𝐽𝛼𝛼𝐽J_{\alpha\alpha}=Jitalic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = italic_J contacts, whereas a boundary flip involves the loss of two contacts (2J=42𝐽42J=42 italic_J = 4). The nucleation of a domain and the motion of an interface are therefore cooperative events, and they strongly affect the state ratios.

Dynamic phase diagrams are shown in Fig. 6. In addition to the SW and HC modes, new one-phase modes Eα, in which the state s=α𝑠𝛼s=\alphaitalic_s = italic_α is predominant as in equilibrium, appear when either h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT or h20subscript20h_{20}italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT dominates. Once a one-phase mode Eα is installed, it does not transform into the SW mode nor into the next state s=[α+1]𝑠delimited-[]𝛼1s=[\alpha+1]italic_s = [ italic_α + 1 ] of the cycling process. Temporal coexistence SW+Eα are possible, as shown in Fig. 6. The various modes are distinguished as follows. When pthree>0.05subscript𝑝three0.05p_{\mathrm{three}}>0.05italic_p start_POSTSUBSCRIPT roman_three end_POSTSUBSCRIPT > 0.05 (see its definition in Sec. III), the system is in the SW mode. When pone>0.05subscript𝑝one0.05p_{\mathrm{one}}>0.05italic_p start_POSTSUBSCRIPT roman_one end_POSTSUBSCRIPT > 0.05, the system is either in the HC mode or in one of the Eα modes. In the case both of them exceed 0.050.050.050.05, the system is either in the coexistence mode SW+HC or in a coexistence mode SW+Eα. The HC and Eα modes are distinguished from the distribution of Nssubscript𝑁𝑠N_{s}italic_N start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT shown in Fig. 7. When a peak exists at Nα/N1subscript𝑁𝛼𝑁1N_{\alpha}/N\approx 1italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N ≈ 1 and it is ten times higher than their local minima close to Nα/N=1subscript𝑁𝛼𝑁1N_{\alpha}/N=1italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N = 1, we consider that the dominant phase of s=α𝑠𝛼s=\alphaitalic_s = italic_α exists in a recognizable period. When all three states take the dominant phase, it is in the HC mode. Otherwise, we consider it to be in the Eα mode, in which the s=α𝑠𝛼s=\alphaitalic_s = italic_α state has the maximum peak at Ns/N1subscript𝑁𝑠𝑁1N_{s}/N\approx 1italic_N start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT / italic_N ≈ 1. The distribution in Fig. 7(b) belongs to E2 but is close to the threshold of the HC mode.

As h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT increases while h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT and h20subscript20h_{20}italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT are fixed, the dynamics changes from E0, to SW+E0, SW, SW+E2, and E2 (or HC instead of E0), in order (see Figs. 6 and 8). Interestingly, high h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT generates the E2 mode, although h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT directly enhances the nucleation and growth of the s=1𝑠1s=1italic_s = 1 phase. This dynamics can be understood by seeing the lifetimes of one-state dominant phases shown in Fig. 8(c). The lifetime of s=α𝑠𝛼s=\alphaitalic_s = italic_α is calculated as the period during which the system stays in Nα>0.98Nsubscript𝑁𝛼0.98𝑁N_{\alpha}>0.98Nitalic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT > 0.98 italic_N. The mean lifetime of the s=0𝑠0s=0italic_s = 0 dominant phase exponentially decreases with increasing h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, by more frequent nucleation and growth of the s=1𝑠1s=1italic_s = 1 phase in the s=0𝑠0s=0italic_s = 0 phase, as expected. Conversely, the mean lifetime of the s=1𝑠1s=1italic_s = 1 phase is almost independent of h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, because the s=1𝑠1s=1italic_s = 1 phase is terminated by the nucleation and growth of the s=2𝑠2s=2italic_s = 2 phase, which is determined by h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT. However, the s=2𝑠2s=2italic_s = 2 phase remains for longer periods with increasing h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, owing to the suppression of the nucleation. To understand this, let us consider an isolated dimer of s=0𝑠0s=0italic_s = 0 sites in the s=2𝑠2s=2italic_s = 2 phase. One site in this dimer can go back to the s=2𝑠2s=2italic_s = 2 state via two pathways. One is the direct flipping to s=2𝑠2s=2italic_s = 2, whose probability is min(1,exp(2Jh20))12𝐽subscript20(1,\exp(2J-h_{20}))( 1 , roman_exp ( 2 italic_J - italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT ) ), since the number of contacts between s=0𝑠0s=0italic_s = 0 and s=2𝑠2s=2italic_s = 2 sites changes from six to four. The other is the two-step flipping via s=1𝑠1s=1italic_s = 1, in which the probabilities of s=01𝑠01s=0\to 1italic_s = 0 → 1 and 12121\to 21 → 2 are min(1,exp(J+h01))1𝐽subscript01(1,\exp(-J+h_{01}))( 1 , roman_exp ( - italic_J + italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ) ) and min(1,exp(3J+h12))13𝐽subscript12(1,\exp(3J+h_{12}))( 1 , roman_exp ( 3 italic_J + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ) ), respectively. The latter flipping increases for higher h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT, since the first step is rate-limiting. Thus, the s=2𝑠2s=2italic_s = 2 phase becomes temporally dominant. Conversely, the acceleration of the second step does not enhance the latter flipping, as seen in the constant lifetime of the s=1𝑠1s=1italic_s = 1 state.

A similar tendency in the flipping energies is also found in the SW mode. When hα[α+1]subscript𝛼delimited-[]𝛼1h_{\alpha[\alpha+1]}italic_h start_POSTSUBSCRIPT italic_α [ italic_α + 1 ] end_POSTSUBSCRIPT is the highest among the three flipping energies, the s=[α+2]𝑠delimited-[]𝛼2s=[\alpha+2]italic_s = [ italic_α + 2 ] state is the most occupied over the lattice (see Figs. 4 and 5). The lifetime of the SW mode is calculated as the lapse from the time entering 0.1N<Nα<0.75N0.1𝑁subscript𝑁𝛼0.75𝑁0.1N<N_{\alpha}<0.75N0.1 italic_N < italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT < 0.75 italic_N for all α𝛼\alphaitalic_α to that entering the one phase mode (Nα>0.98Nsubscript𝑁𝛼0.98𝑁N_{\alpha}>0.98Nitalic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT > 0.98 italic_N for one of the states). When the three forward flips are comparable (e.g., h011similar-to-or-equalssubscript011h_{01}\simeq 1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ≃ 1 while h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1 as in Fig. 8), the mean lifetime of the SW mode becomes longer and the system falls into this mode. As the flipping energies deviate more from the symmetric condition, one of the states (s=α𝑠𝛼s=\alphaitalic_s = italic_α) dominates the SW mode, and subsequently the SW mode changes into the s=α𝑠𝛼s=\alphaitalic_s = italic_α dominant phase, Eα, via the SW+ Eα mode or the HC mode (see Fig. 6).

The phase diagram in Fig. 6(a) is slightly asymmetric with respect to the symmetric axis h01=h12subscript01subscript12h_{01}=h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT. This is due to the fact that the transition from the s=α𝑠𝛼s=\alphaitalic_s = italic_α state to the s=[α+1]𝑠delimited-[]𝛼1s=[\alpha+1]italic_s = [ italic_α + 1 ] state is also dependent on the [α+1][α+2]delimited-[]𝛼1delimited-[]𝛼2[\alpha+1]\to[\alpha+2][ italic_α + 1 ] → [ italic_α + 2 ] flip, as mentioned above. This asymmetry can be more quantitatively captured by plotting the various physical quantities as a function of h01h12subscript01subscript12h_{01}-h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT along the diagonal line fixing h01+h12subscript01subscript12h_{01}+h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT, as shown in Fig. 9. At h01+h12=3subscript01subscript123h_{01}+h_{12}=3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 3 and h20=1subscript201h_{20}=1italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1, the SW mode more often occurs for h01h12<0subscript01subscript120h_{01}-h_{12}<0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT < 0 than for h01h12>0subscript01subscript120h_{01}-h_{12}>0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT > 0, and the center of the SW existence region is at h01h120.1similar-to-or-equalssubscript01subscript120.1h_{01}-h_{12}\simeq-0.1italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ≃ - 0.1 (see Figs. 9(a) and (b)). The flow qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT has a maximum around this center (compare the red line in Fig. 9(c) and the blue line in Fig. 9(b)). With increasing h01+h12subscript01subscript12h_{01}+h_{12}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT, the center position is shifted to the positive direction and the maximum of qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT becomes higher (see the three solid lines in Fig. 9(c)). At higher h20subscript20h_{20}italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT, qfsubscript𝑞fq_{\mathrm{f}}italic_q start_POSTSUBSCRIPT roman_f end_POSTSUBSCRIPT increases (compare the top two lines in Fig. 9(c)). Thus, the asymmetry direction and amplitude vary depending on the conditions.

Refer to caption
Figure 11: Effects of the asymmetric contact energies with the symmetric flipping energies h=h01=h12=h20subscript01subscript12subscript20h=h_{01}=h_{12}=h_{20}italic_h = italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT at J11=J22=2subscript𝐽11subscript𝐽222J_{11}=J_{22}=2italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = 2 and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a)–(c) Dependence on the contact energy J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT at h=h01=h12=h20=1.3subscript01subscript12subscript201.3h=h_{01}=h_{12}=h_{20}=1.3italic_h = italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1.3. The dynamic modes are described at the top of the figure. (a) Probabilities of three states. The gray up-pointing triangles or squares, blue down-pointing triangles or crosses, and red circles or diamonds represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively. (b) Probabilities of one-phase state and three-phase coexistence. (c) Lifetimes of the homogeneous phases (s=0𝑠0s=0italic_s = 0, 1111, and 2222) and SW mode. The dashed lines in (a–c) represent kinetically trapped modes. (d) Dynamic phase diagram as a function of J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT. The squares (light blue) and circles (gray) represent the HC and SW modes, respectively. The up-pointing and down-pointing triangles (green) represent the homogeneous phases of s=0𝑠0s=0italic_s = 0 and s=2𝑠2s=2italic_s = 2 (E0 and E2), respectively. The pluses (blue), open diamonds (magenta), closed diamonds (magenta), and crosses (red) represent the temporal coexistence: SW+HC, SW+E0, SW+E1, SW+E2, respectively. The asterisks indicate the hysteresis of two modes (E0 and SW+E1 or E1).

Last, we consider a large system of N=2562𝑁superscript2562N=256^{2}italic_N = 256 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT sites, in which the HC mode only exists in the hysteresis at low hhitalic_h in symmetric conditions. The effects of the asymmetric flipping energies are similar to those in the small system (N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT). As h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT increases, the number ratio N2/Nsubscript𝑁2𝑁N_{2}/Nitalic_N start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / italic_N of the s=2𝑠2s=2italic_s = 2 state increases in the SW mode, and subsequently, the E2 phase appears via the temporal coexistence with the SW mode (SW+E2), as shown in Fig. 10. Since the size of biphasic domains and the number density of three-contact points continuously decrease with increasing N2/Nsubscript𝑁2𝑁N_{2}/Nitalic_N start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT / italic_N, the transition between the SW and Eα modes likely occurs continuously via the mode coexistence at very large N𝑁Nitalic_N.

V Dynamics for Asymmetric Contact Energies

In this section, we describe the dynamics in the case of asymmetric contact energies J00J11=J22=Jsubscript𝐽00subscript𝐽11subscript𝐽22𝐽J_{00}\neq J_{11}=J_{22}=Jitalic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT ≠ italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = italic_J and symmetric flipping energies h01=h12=h20=hsubscript01subscript12subscript20h_{01}=h_{12}=h_{20}=hitalic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h, for N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. When J00>Jsubscript𝐽00𝐽J_{00}>Jitalic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT > italic_J, the s=0𝑠0s=0italic_s = 0 state is stabilized and becomes dominant at high J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT, so that the E0 mode appears (Fig. 11). Conversely, with decreasing J00<Jsubscript𝐽00𝐽J_{00}<Jitalic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT < italic_J, the E2 mode appears since the nucleation of a s=0𝑠0s=0italic_s = 0 domain is suppressed due to the low contact energy in the domain. This is captured by the lifetimes of one-state dominant phases, as in the case of the asymmetric flipping energies (see Fig. 11(c)).

Refer to caption
Figure 12: Coexistence of the SW and E0 modes at J00=2.5subscript𝐽002.5J_{00}=2.5italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT = 2.5, J11=J22=2subscript𝐽11subscript𝐽222J_{11}=J_{22}=2italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT = 2, h=h01=h12=h20=1.3subscript01subscript12subscript201.3h=h_{01}=h_{12}=h_{20}=1.3italic_h = italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1.3, and N=1282𝑁superscript1282N=128^{2}italic_N = 128 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT. (a) Sequential snapshots in the SW mode. A biphasic domain moves in the direction of the s=2𝑠2s=2italic_s = 2 region. (b) Time evolution of numbers of three states (corresponding to snapshots in (a)). The light gray, blue (dark gray), and red (medium gray) colors represent the data for s=0𝑠0s=0italic_s = 0, 1111, and 2222, respectively. (c) Time evolution of numbers of three states starting from a s=0𝑠0s=0italic_s = 0 dominant phase (E0 mode).

Interestingly, hysteresis occurs at high hhitalic_h and high J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT. Two different modes appear depending on the initial condition, as shown in Fig. 12 and Movie S4: the E0 mode and the SW+E1 mode. They do not change into each other in our simulation times 108similar-toabsentsuperscript108\sim 10^{8}∼ 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT, since their lifetimes become longer due to high energy barriers for the transitions. In the s=1𝑠1s=1italic_s = 1 phase, biphasic domains of s=2𝑠2s=2italic_s = 2 and s=0𝑠0s=0italic_s = 0 often appear but do not grow to cover the entire lattice (see Figs. 12(a) and (b)). At relatively low hhitalic_h, the SW+E1 branch discretely appears in the middle of the E0 region (see Fig. 11). Thus, the SW+E1 mode is likely a kinetically trapped metastable state but the E0 mode remains as the only solution in the long time limit. In contrast, at high hhitalic_h, the E0 branch is only connected with the E0-mode solution at high J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT, and the SW+E1 branch is continuously connected with the SW+E1-mode solution at low J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT (see the region of h1.51.5h\geq 1.5italic_h ≥ 1.5 in Fig. 11(d)). The SW+E1 mode becomes the E1 mode at high J00subscript𝐽00J_{00}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT (J002.5subscript𝐽002.5J_{00}\geq 2.5italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT ≥ 2.5 at h=1.51.5h=1.5italic_h = 1.5). Thus, the stable mode in the long time limit likely changes from the SW+E1 or E1 mode to E0 mode in the middle of the hysteresis region.

In the symmetric condition (J00=J11=J22subscript𝐽00subscript𝐽11subscript𝐽22J_{00}=J_{11}=J_{22}italic_J start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = italic_J start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT and h01=h12=h20=hsubscript01subscript12subscript20h_{01}=h_{12}=h_{20}=hitalic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h), three dominant phases (s=0𝑠0s=0italic_s = 0, 1111, and 2222) are trapped in the limit h1much-less-than1h\ll 1italic_h ≪ 1. In particular, they are degenerated thermal-equilibrium states at h=00h=0italic_h = 0. The hysteresis found here is a nonequilibrium extension of this degeneration.

Refer to caption
Figure 13: Density of the states s=0𝑠0s=0italic_s = 0, 1111, and 2222 as a function of h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT in an homogeneously mixed phase for h12=h20=1subscript12subscript201h_{12}=h_{20}=1italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 1 and Jαα=0subscript𝐽𝛼𝛼0J_{\alpha\alpha}=0italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 0. Note that h01+h12+h200subscript01subscript12subscript200h_{01}+h_{12}+h_{20}\neq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT ≠ 0, so the system is out of equilibrium. The solid and dashed lines represent the data for Metropolis rates and Glauber rates, respectively.

VI Theoretical Analysis

VI.1 Homogeneous Mixed State in the Absence of Interactions

Let us call ρα=Nα/Nsubscript𝜌𝛼subscript𝑁𝛼𝑁\rho_{\alpha}=N_{\alpha}/Nitalic_ρ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT = italic_N start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_N the number density of the state s=α𝑠𝛼s=\alphaitalic_s = italic_α. Since there are no empty sites, ρ0+ρ1+ρ2=1subscript𝜌0subscript𝜌1subscript𝜌21\rho_{0}+\rho_{1}+\rho_{2}=1italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 1. We work in units such that the thermal energy kBTsubscript𝑘B𝑇k_{\mathrm{B}}Titalic_k start_POSTSUBSCRIPT roman_B end_POSTSUBSCRIPT italic_T and the size of the lattice sites are unity. Let us first consider a simplified situation in which we disregard all spatial structures and where the interactions between adjacent sites are neglected, i.e., Jαα=0subscript𝐽𝛼𝛼0J_{\alpha\alpha}=0italic_J start_POSTSUBSCRIPT italic_α italic_α end_POSTSUBSCRIPT = 0. In agreement with Eq. (1), the energy density (per site) is given by

fself=ρ1ε1+ρ2ε2+(1ρ1ρ2)ε0,subscript𝑓selfsubscript𝜌1subscript𝜀1subscript𝜌2subscript𝜀21subscript𝜌1subscript𝜌2subscript𝜀0f_{\mathrm{self}}=\rho_{1}\varepsilon_{1}+\rho_{2}\varepsilon_{2}+(1-\rho_{1}-% \rho_{2})\varepsilon_{0},italic_f start_POSTSUBSCRIPT roman_self end_POSTSUBSCRIPT = italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + ( 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , (3)

and the flipping energies hαβsubscript𝛼𝛽h_{\alpha\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT are εαεβsubscript𝜀𝛼subscript𝜀𝛽\varepsilon_{\alpha}-\varepsilon_{\beta}italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT in equilibrium.

As mentioned earlier, we place the system in an arbitrary nonequilibrium state, by setting hαβ=εαεβ+hαβneqsubscript𝛼𝛽subscript𝜀𝛼subscript𝜀𝛽subscriptsuperscriptneq𝛼𝛽h_{\alpha\beta}=\varepsilon_{\alpha}-\varepsilon_{\beta}+h^{\mathrm{neq}}_{% \alpha\beta}italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT + italic_h start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT, with hβα=hαβsubscript𝛽𝛼subscript𝛼𝛽h_{\beta\alpha}=-h_{\alpha\beta}italic_h start_POSTSUBSCRIPT italic_β italic_α end_POSTSUBSCRIPT = - italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT, while assuming that the transition rates wαβsubscript𝑤𝛼𝛽w_{\alpha\beta}italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT obey the local detailed balance condition: wαβ/wβα=exp(hαβ)subscript𝑤𝛼𝛽subscript𝑤𝛽𝛼subscript𝛼𝛽w_{\alpha\beta}/w_{\beta\alpha}=\exp(h_{\alpha\beta})italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT / italic_w start_POSTSUBSCRIPT italic_β italic_α end_POSTSUBSCRIPT = roman_exp ( italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ). We therefore allow for h01+h12+h20=hcyc0subscript01subscript12subscript20subscriptcyc0h_{01}+h_{12}+h_{20}=h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0, with hcyc=h01neq+h12neq+h20neqsubscriptcycsubscriptsuperscriptneq01subscriptsuperscriptneq12subscriptsuperscriptneq20h_{\mathrm{cyc}}=h^{\mathrm{neq}}_{01}+h^{\mathrm{neq}}_{12}+h^{\mathrm{neq}}_% {20}italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = italic_h start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_h start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_h start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT, which we assume positive to favor the cycling s=0120𝑠0120s=0\to 1\to 2\to 0italic_s = 0 → 1 → 2 → 0. In practice, we take Metropolis rates wαβ=min[1,exp(hαβ)]subscript𝑤𝛼𝛽1subscript𝛼𝛽w_{\alpha\beta}=\min[1,\exp(h_{\alpha\beta})]italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = roman_min [ 1 , roman_exp ( italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ) ] or Glauber rates wαβ=[1+exp(hαβ)]1subscript𝑤𝛼𝛽superscriptdelimited-[]1subscript𝛼𝛽1w_{\alpha\beta}=[1+\exp(-h_{\alpha\beta})]^{-1}italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT = [ 1 + roman_exp ( - italic_h start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT.

The dynamical equations are then

ρ˙1subscript˙𝜌1\displaystyle\dot{\rho}_{1}over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT =w01ρ0(w10+w12)ρ1+w21ρ2,absentsubscript𝑤01subscript𝜌0subscript𝑤10subscript𝑤12subscript𝜌1subscript𝑤21subscript𝜌2\displaystyle=w_{01}\rho_{0}-\left(w_{10}+w_{12}\right)\rho_{1}+w_{21}\rho_{2},= italic_w start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - ( italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ) italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , (4)
ρ˙2subscript˙𝜌2\displaystyle\dot{\rho}_{2}over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT =w02ρ0+w12ρ1(w20+w21)ρ2,absentsubscript𝑤02subscript𝜌0subscript𝑤12subscript𝜌1subscript𝑤20subscript𝑤21subscript𝜌2\displaystyle=w_{02}\rho_{0}+w_{12}\rho_{1}-\left(w_{20}+w_{21}\right)\rho_{2},= italic_w start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - ( italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT ) italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , (5)
ρ˙0subscript˙𝜌0\displaystyle\dot{\rho}_{0}over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT =(w01+w02)ρ0+w10ρ1+w20ρ2.absentsubscript𝑤01subscript𝑤02subscript𝜌0subscript𝑤10subscript𝜌1subscript𝑤20subscript𝜌2\displaystyle=-\left(w_{01}+w_{02}\right)\rho_{0}+w_{10}\rho_{1}+w_{20}\rho_{2}.= - ( italic_w start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT ) italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT . (6)

The third equation follows from the first two due to the density conservation equation. In the stationary state, ρ˙α=0subscript˙𝜌𝛼0\dot{\rho}_{\alpha}=0over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT = 0, hence replacing ρ0subscript𝜌0\rho_{0}italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT by 1ρ1ρ21subscript𝜌1subscript𝜌21-\rho_{1}-\rho_{2}1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT in Eqs.(4)–(5) yields a linear system for the stationary densities ρ1subscript𝜌1\rho_{1}italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and ρ2subscript𝜌2\rho_{2}italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT, which gives

ρ1ρ0subscript𝜌1subscript𝜌0\displaystyle\frac{\rho_{1}}{\rho_{0}}divide start_ARG italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG =w02w21+w01(w20+w21)w12w20+w10(w20+w21),absentsubscript𝑤02subscript𝑤21subscript𝑤01subscript𝑤20subscript𝑤21subscript𝑤12subscript𝑤20subscript𝑤10subscript𝑤20subscript𝑤21\displaystyle=\frac{w_{02}w_{21}+w_{01}(w_{20}+w_{21})}{w_{12}w_{20}+w_{10}(w_% {20}+w_{21})},= divide start_ARG italic_w start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ( italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT ) end_ARG start_ARG italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ( italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT ) end_ARG , (7)
ρ2ρ0subscript𝜌2subscript𝜌0\displaystyle\frac{\rho_{2}}{\rho_{0}}divide start_ARG italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG =w01w12+w02(w10+w12)w12w20+w10(w20+w21).absentsubscript𝑤01subscript𝑤12subscript𝑤02subscript𝑤10subscript𝑤12subscript𝑤12subscript𝑤20subscript𝑤10subscript𝑤20subscript𝑤21\displaystyle=\frac{w_{01}w_{12}+w_{02}(w_{10}+w_{12})}{w_{12}w_{20}+w_{10}(w_% {20}+w_{21})}.= divide start_ARG italic_w start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT ( italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ) end_ARG start_ARG italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ( italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT ) end_ARG . (8)

For Metropolis rates, assuming h010subscript010h_{01}\geq 0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ≥ 0, h120subscript120h_{12}\geq 0italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ≥ 0, and h200subscript200h_{20}\geq 0italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT ≥ 0, we obtain

ρ1ρ0subscript𝜌1subscript𝜌0\displaystyle\frac{\rho_{1}}{\rho_{0}}divide start_ARG italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG =eh011+eh12+eh021+eh12+eh02ehcyc,absentsuperscript𝑒subscript011superscript𝑒subscript12superscript𝑒subscript021superscript𝑒subscript12superscript𝑒subscript02superscript𝑒subscriptcyc\displaystyle=e^{h_{01}}\frac{1+e^{h_{12}}+e^{h_{02}}}{1+e^{h_{12}}+e^{h_{02}}% e^{h_{\mathrm{cyc}}}},= italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT end_POSTSUPERSCRIPT end_ARG , (9)
ρ0ρ2subscript𝜌0subscript𝜌2\displaystyle\frac{\rho_{0}}{\rho_{2}}divide start_ARG italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG =eh201+eh01+eh211+eh01+eh21ehcyc.absentsuperscript𝑒subscript201superscript𝑒subscript01superscript𝑒subscript211superscript𝑒subscript01superscript𝑒subscript21superscript𝑒subscriptcyc\displaystyle=e^{h_{20}}\frac{1+e^{h_{01}}+e^{h_{21}}}{1+e^{h_{01}}+e^{h_{21}}% e^{h_{\mathrm{cyc}}}}.= italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT divide start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT + italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_e start_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT end_POSTSUPERSCRIPT end_ARG . (10)

Equilibrium is achieved for hcyc=0subscriptcyc0h_{\mathrm{cyc}}=0italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = 0. In this case, Eq. (9) give ρ1/ρ0=exp(h01)subscript𝜌1subscript𝜌0subscript01\rho_{1}/\rho_{0}=\exp(h_{01})italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT / italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = roman_exp ( italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ) and ρ0/ρ2=exp(h02)subscript𝜌0subscript𝜌2subscript02\rho_{0}/\rho_{2}=\exp(h_{02})italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT / italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = roman_exp ( italic_h start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT ), which produces the Boltzmann distribution ραexp(εα)proportional-tosubscript𝜌𝛼subscript𝜀𝛼\rho_{\alpha}\propto\exp(-\varepsilon_{\alpha})italic_ρ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ∝ roman_exp ( - italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ), as expected.

When hcyc0subscriptcyc0h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0, the stationary densities differ from the equilibrium ones. As h01subscript01h_{01}italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT increases at fixed h12subscript12h_{12}italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT and h20subscript20h_{20}italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT, we find that ρ1subscript𝜌1\rho_{1}italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT increases and ρ0subscript𝜌0\rho_{0}italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT decreases, as expected (see Fig. 13). When Glauber rates are used, the density changes become slightly larger (compare the dashed and solid lines in Fig. 13).

VI.1.1 Equilibrium State and Free Energy

At equilibrium, i.e., for hcyc=0subscriptcyc0h_{\mathrm{cyc}}=0italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = 0, the densities can be obtained from the free-energy density f=fself+fmix𝑓subscript𝑓selfsubscript𝑓mixf=f_{\mathrm{self}}+f_{\mathrm{mix}}italic_f = italic_f start_POSTSUBSCRIPT roman_self end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_mix end_POSTSUBSCRIPT, with

fmix=ρ1lnρ1+ρ2lnρ2+(1ρ1ρ2)ln(1ρ1ρ2),subscript𝑓mixsubscript𝜌1𝜌1subscript𝜌2subscript𝜌21subscript𝜌1subscript𝜌21subscript𝜌1subscript𝜌2f_{\mathrm{mix}}=\rho_{1}\ln\rho 1+\rho_{2}\ln\rho_{2}+(1-\rho_{1}-\rho_{2})% \ln(1-\rho_{1}-\rho_{2}),italic_f start_POSTSUBSCRIPT roman_mix end_POSTSUBSCRIPT = italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT roman_ln italic_ρ 1 + italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT roman_ln italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + ( 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) roman_ln ( 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) , (11)

arising from the entropy of mixing. Indeed, minimizing f𝑓fitalic_f with respect to ρ1subscript𝜌1\rho_{1}italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and ρ2subscript𝜌2\rho_{2}italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT yields ρα/ρ0=exp(ε0εα)subscript𝜌𝛼subscript𝜌0subscript𝜀0subscript𝜀𝛼\rho_{\alpha}/\rho_{0}=\exp(\varepsilon_{0}-\varepsilon_{\alpha})italic_ρ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT / italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = roman_exp ( italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ).

Note that while the equilibrium densities are obtained by minimizing the full free energy density f𝑓fitalic_f, the transition rates wαβsubscript𝑤𝛼𝛽w_{\alpha\beta}italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT result from the variations of fself=ffmixsubscript𝑓self𝑓subscript𝑓mixf_{\mathrm{self}}=f-f_{\mathrm{mix}}italic_f start_POSTSUBSCRIPT roman_self end_POSTSUBSCRIPT = italic_f - italic_f start_POSTSUBSCRIPT roman_mix end_POSTSUBSCRIPT.

Refer to caption
Figure 14: Continuous model with h01=h12=h20=0.3subscript01subscript12subscript200.3h_{01}=h_{12}=h_{20}=0.3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 0.3, obtained from ε0=0subscript𝜀00\varepsilon_{0}=0italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, ε1=0.3subscript𝜀10.3\varepsilon_{1}=-0.3italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - 0.3, ε2=0.6subscript𝜀20.6\varepsilon_{2}=-0.6italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 0.6, and h20neq=hcyc=0.9superscriptsubscript20neqsubscriptcyc0.9h_{20}^{\mathrm{neq}}=h_{\mathrm{cyc}}=0.9italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = 0.9. The interaction and gradient coefficients are u=4𝑢4u=4italic_u = 4 and k=1𝑘1k=1italic_k = 1, respectively. (a) Free energy landscape f𝑓fitalic_f as a function of ρ1subscript𝜌1\rho_{1}italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and ρ2subscript𝜌2\rho_{2}italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT. The bottom left, bottom right and top left minima correspond to the phases s=0𝑠0s=0italic_s = 0, s=1𝑠1s=1italic_s = 1, and s=2𝑠2s=2italic_s = 2, respectively. (b) Numerical integration of Eqs. (12)–(14) in one dimension, showing a biphasic s=2,1𝑠21s=2,1italic_s = 2 , 1 (red, blue) solitary wave in an s=0𝑠0s=0italic_s = 0 (gray) background. Time flows from top to bottom. The time lapse between two snapshots: 4×1054superscript1054\times 10^{5}4 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT time steps dt=1.25×104𝑑𝑡1.25superscript104dt=1.25\times 10^{-4}italic_d italic_t = 1.25 × 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT. The RGB color encodes the system state as ρ0𝒢+ρ1+ρ2subscript𝜌0𝒢subscript𝜌1subscript𝜌2\rho_{0}\mathcal{G}+\rho_{1}\mathcal{B}+\rho_{2}\mathcal{R}italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT caligraphic_G + italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT caligraphic_B + italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT caligraphic_R, with 𝒢𝒢\mathcal{G}caligraphic_G gray, \mathcal{B}caligraphic_B blue, and \mathcal{R}caligraphic_R red. (c) Two different-sized bands moving in opposite directions, colliding and annihilating each other. The same parameters as in the previous image.

VI.2 Continuum Theory with Spatial Gradients and Interactions

We now assume that the densities ρi(𝐱)subscript𝜌𝑖𝐱\rho_{i}(\mathbf{x})italic_ρ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x ) vary in space. Since we consider only local state flips (no state exchange between adjacent sites), the dynamical equations still have the same form locally:

ρ˙1(𝐱)subscript˙𝜌1𝐱\displaystyle\dot{\rho}_{1}(\mathbf{x})over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( bold_x ) =w01ρ0(w10+w12)ρ1+w21ρ2,absentsubscript𝑤01subscript𝜌0subscript𝑤10subscript𝑤12subscript𝜌1subscript𝑤21subscript𝜌2\displaystyle=w_{01}\rho_{0}-\left(w_{10}+w_{12}\right)\rho_{1}+w_{21}\rho_{2},= italic_w start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - ( italic_w start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT ) italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , (12)
ρ˙2(𝐱)subscript˙𝜌2𝐱\displaystyle\dot{\rho}_{2}(\mathbf{x})over˙ start_ARG italic_ρ end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ( bold_x ) =w02ρ0+w12ρ1(w20+w21)ρ2,absentsubscript𝑤02subscript𝜌0subscript𝑤12subscript𝜌1subscript𝑤20subscript𝑤21subscript𝜌2\displaystyle=w_{02}\rho_{0}+w_{12}\rho_{1}-\left(w_{20}+w_{21}\right)\rho_{2},= italic_w start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - ( italic_w start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT + italic_w start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT ) italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , (13)
ρ0(𝐱)subscript𝜌0𝐱\displaystyle\rho_{0}(\mathbf{x})italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( bold_x ) =1ρ1ρ2.absent1subscript𝜌1subscript𝜌2\displaystyle=1-\rho_{1}-\rho_{2}.= 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT . (14)

The transition rates wij(𝐱)subscript𝑤𝑖𝑗𝐱w_{ij}(\mathbf{x})italic_w start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ( bold_x ), however, now depend on the local free-energy variations, including the interactions between adjacent sites and the penalty due to density gradients. The latter are required in order to take into account that the nucleation of a new phase behaves differently within a domain and at the interface between two domains. Thus, the free-energy density is now f=fself+fmix+fint+fgrad𝑓subscript𝑓selfsubscript𝑓mixsubscript𝑓intsubscript𝑓gradf=f_{\mathrm{self}}+f_{\mathrm{mix}}+f_{\mathrm{int}}+f_{\mathrm{grad}}italic_f = italic_f start_POSTSUBSCRIPT roman_self end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_mix end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_grad end_POSTSUBSCRIPT, with

fintsubscript𝑓int\displaystyle f_{\mathrm{int}}italic_f start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT =u12ρ12u22ρ22u02(1ρ1ρ2)2,absentsubscript𝑢12superscriptsubscript𝜌12subscript𝑢22superscriptsubscript𝜌22subscript𝑢02superscript1subscript𝜌1subscript𝜌22\displaystyle=-\frac{u_{1}}{2}\rho_{1}^{2}-\frac{u_{2}}{2}\rho_{2}^{2}-\frac{u% _{0}}{2}\left(1-\rho_{1}-\rho_{2}\right)^{2},= - divide start_ARG italic_u start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - divide start_ARG italic_u start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - divide start_ARG italic_u start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ( 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (15)
fgradsubscript𝑓grad\displaystyle f_{\mathrm{grad}}italic_f start_POSTSUBSCRIPT roman_grad end_POSTSUBSCRIPT =k2(ρ1)2+k2(ρ2)2+k2[(1ρ1ρ2)]2,absent𝑘2superscriptsubscript𝜌12𝑘2superscriptsubscript𝜌22𝑘2superscriptdelimited-[]1subscript𝜌1subscript𝜌22\displaystyle=\frac{k}{2}\left(\nabla\rho_{1}\right)^{2}+\frac{k}{2}\left(% \nabla\rho_{2}\right)^{2}+\frac{k}{2}\left[\nabla\left(1-\rho_{1}-\rho_{2}% \right)\right]^{2},\hskip 11.38092pt= divide start_ARG italic_k end_ARG start_ARG 2 end_ARG ( ∇ italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG italic_k end_ARG start_ARG 2 end_ARG ( ∇ italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG italic_k end_ARG start_ARG 2 end_ARG [ ∇ ( 1 - italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (16)

where uα>0subscript𝑢𝛼0u_{\alpha}>0italic_u start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT > 0 quantifies the attractive couplings between like states, and k>0𝑘0k>0italic_k > 0 the penalty associated with density gradients (common to all states for simplicity).

As discussed in Sec. VI.1.1, the transition rates must be calculated from the free energy minus the contribution arising from the entropy of mixing. The latter is given by the functional

F[ρ1,ρ2]=d2x(fself+fint+fgrad).𝐹subscript𝜌1subscript𝜌2superscript𝑑2𝑥subscript𝑓selfsubscript𝑓intsubscript𝑓gradF[\rho_{1},\rho_{2}]=\int\!d^{2}x\left(f_{\mathrm{self}}+f_{\mathrm{int}}+f_{% \mathrm{grad}}\right).italic_F [ italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] = ∫ italic_d start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_x ( italic_f start_POSTSUBSCRIPT roman_self end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_int end_POSTSUBSCRIPT + italic_f start_POSTSUBSCRIPT roman_grad end_POSTSUBSCRIPT ) . (17)

The energy change associated with the flip of one site from the state s=01𝑠01s=0\to 1italic_s = 0 → 1, at point 𝐱𝐱\mathbf{x}bold_x, supplemented by the nonequilibrium contribution h10neqsuperscriptsubscript10neqh_{10}^{\mathrm{neq}}italic_h start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT, is given by

ΔH01(𝐱)Δsubscript𝐻01𝐱\displaystyle\Delta H_{01}(\mathbf{x})roman_Δ italic_H start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT ( bold_x ) =δFδρ1(𝐱)|ρ2=ε1ε0u1ρ1+u0ρ0absentevaluated-at𝛿𝐹𝛿subscript𝜌1𝐱subscript𝜌2subscript𝜀1subscript𝜀0subscript𝑢1subscript𝜌1subscript𝑢0subscript𝜌0\displaystyle=\left.\frac{\delta F}{\delta\rho_{1}(\mathbf{x})}\right|_{\rho_{% 2}}=\varepsilon_{1}-\varepsilon_{0}-u_{1}\rho_{1}+u_{0}\rho_{0}= divide start_ARG italic_δ italic_F end_ARG start_ARG italic_δ italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( bold_x ) end_ARG | start_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_u start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_u start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT (18)
k2(ρ1ρ0)h01neq=ΔH10(𝐱).𝑘superscript2subscript𝜌1subscript𝜌0superscriptsubscript01neqΔsubscript𝐻10𝐱\displaystyle-k\nabla^{2}\left(\rho_{1}-\rho_{0}\right)-h_{01}^{\mathrm{neq}}=% -\Delta H_{10}(\mathbf{x}).- italic_k ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) - italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = - roman_Δ italic_H start_POSTSUBSCRIPT 10 end_POSTSUBSCRIPT ( bold_x ) .

Likewise, ΔH02=ε2ε0u2ρ2+u0ρ0k2(ρ2ρ0)h02neq=ΔH20Δsubscript𝐻02subscript𝜀2subscript𝜀0subscript𝑢2subscript𝜌2subscript𝑢0subscript𝜌0𝑘superscript2subscript𝜌2subscript𝜌0superscriptsubscript02neqΔsubscript𝐻20\Delta H_{02}=\varepsilon_{2}-\varepsilon_{0}-u_{2}\rho_{2}+u_{0}\rho_{0}-k% \nabla^{2}(\rho_{2}-\rho_{0})-h_{02}^{\mathrm{neq}}=-\Delta H_{20}roman_Δ italic_H start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_u start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_u start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - italic_k ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) - italic_h start_POSTSUBSCRIPT 02 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = - roman_Δ italic_H start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT. From the energy change in the sequence s=201𝑠201s=2\to 0\to 1italic_s = 2 → 0 → 1, we infer also ΔH21=ε1ε2u1ρ1+u2ρ2k2(ρ1ρ2)h21neq=ΔH12Δsubscript𝐻21subscript𝜀1subscript𝜀2subscript𝑢1subscript𝜌1subscript𝑢2subscript𝜌2𝑘superscript2subscript𝜌1subscript𝜌2superscriptsubscript21neqΔsubscript𝐻12\Delta H_{21}=\varepsilon_{1}-\varepsilon_{2}-u_{1}\rho_{1}+u_{2}\rho_{2}-k% \nabla^{2}(\rho_{1}-\rho_{2})-h_{21}^{\mathrm{neq}}=-\Delta H_{12}roman_Δ italic_H start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_u start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_u start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_k ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) - italic_h start_POSTSUBSCRIPT 21 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = - roman_Δ italic_H start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT.

In the following, we shall take u0=u1=u2=usubscript𝑢0subscript𝑢1subscript𝑢2𝑢u_{0}=u_{1}=u_{2}=uitalic_u start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_u start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_u start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = italic_u and break detailed balance by setting h01neq=0superscriptsubscript01neq0h_{01}^{\mathrm{neq}}=0italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = 0, h12neq=0superscriptsubscript12neq0h_{12}^{\mathrm{neq}}=0italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = 0, and h20neq=hcyc>0superscriptsubscript20neqsubscriptcyc0h_{20}^{\mathrm{neq}}=h_{\mathrm{cyc}}>0italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT > 0. For convenience, we assume Glauber rates:

wαβ(𝐱)=11+eΔHαβ(𝐱).subscript𝑤𝛼𝛽𝐱11superscript𝑒Δsubscript𝐻𝛼𝛽𝐱w_{\alpha\beta}(\mathbf{x})=\frac{1}{1+e^{\Delta H_{\alpha\beta}(\mathbf{x})}}.italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ( bold_x ) = divide start_ARG 1 end_ARG start_ARG 1 + italic_e start_POSTSUPERSCRIPT roman_Δ italic_H start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ( bold_x ) end_POSTSUPERSCRIPT end_ARG . (19)

Note that they obey the local detailed balance condition wαβ/wβα=exp(ΔHαβ)subscript𝑤𝛼𝛽subscript𝑤𝛽𝛼Δsubscript𝐻𝛼𝛽w_{\alpha\beta}/w_{\beta\alpha}=\exp(-\Delta H_{\alpha\beta})italic_w start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT / italic_w start_POSTSUBSCRIPT italic_β italic_α end_POSTSUBSCRIPT = roman_exp ( - roman_Δ italic_H start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ), even in the nonequilibrium case hcyc0subscriptcyc0h_{\mathrm{cyc}}\neq 0italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT ≠ 0.

Refer to caption
Figure 15: Evolution of a biphasic s=1,2𝑠12s=1,2italic_s = 1 , 2 band in a s=0𝑠0s=0italic_s = 0 backgroud, for asymmetric flipping energies. (a) Weak asymmetry: h01=0.34subscript010.34h_{01}=0.34italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 0.34, h12=0.3subscript120.3h_{12}=0.3italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 0.3, and h20=0.3subscript200.3h_{20}=0.3italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 0.3, obtained from ε0=0subscript𝜀00\varepsilon_{0}=0italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, ε1=0.34subscript𝜀10.34\varepsilon_{1}=-0.34italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - 0.34, ε2=0.64subscript𝜀20.64\varepsilon_{2}=-0.64italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 0.64, and h20neq=hcyc=0.94superscriptsubscript20neqsubscriptcyc0.94h_{20}^{\mathrm{neq}}=h_{\mathrm{cyc}}=0.94italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = 0.94. (b) Large asymmetry: h01=0.4subscript010.4h_{01}=0.4italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = 0.4, h12=0.45subscript120.45h_{12}=0.45italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = 0.45, and h20=0.3subscript200.3h_{20}=0.3italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 0.3, obtained from ε0=0subscript𝜀00\varepsilon_{0}=0italic_ε start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, ε1=0.4subscript𝜀10.4\varepsilon_{1}=-0.4italic_ε start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = - 0.4, ε2=0.85subscript𝜀20.85\varepsilon_{2}=-0.85italic_ε start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = - 0.85, and h20neq=hcyc=1.15superscriptsubscript20neqsubscriptcyc1.15h_{20}^{\mathrm{neq}}=h_{\mathrm{cyc}}=1.15italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_neq end_POSTSUPERSCRIPT = italic_h start_POSTSUBSCRIPT roman_cyc end_POSTSUBSCRIPT = 1.15.

VI.2.1 Numerical Analysis

Equations (12)–(14) with the rates of Eq. (19), obtained from the energy variations ΔHαβ(𝐱)Δsubscript𝐻𝛼𝛽𝐱\Delta H_{\alpha\beta}(\mathbf{x})roman_Δ italic_H start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ( bold_x ) defined in Eq. (18) and below, allow to calculate the spatially resolved evolution of the system. We have integrated the discretized version of these equations in one dimension with periodic boundary conditions in a lattice of L=120𝐿120L=120italic_L = 120 sites using the classical Runge-Kutta method. The first and second derivatives were discretized using central differences with fourth-order accuracy.

We first investigated the dynamics in cyclic-symmetry conditions, with h01=h12=h20=0.3subscript01subscript12subscript200.3h_{01}=h_{12}=h_{20}=0.3italic_h start_POSTSUBSCRIPT 01 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT 20 end_POSTSUBSCRIPT = 0.3, in which case the three states are equivalent and the dynamics promotes the sequence s=0120𝑠0120s=0\to 1\to 2\to 0italic_s = 0 → 1 → 2 → 0. We show in Fig. 14(b) that a biphasic s=1,2𝑠12s=1,2italic_s = 1 , 2 band in a s=0𝑠0s=0italic_s = 0 background forms a travelling wave, which behaves as a perfect soliton. The biphasic band is initialized with (ρ1,ρ2)subscript𝜌1subscript𝜌2(\rho_{1},\rho_{2})( italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) equals (0,0)00(0,0)( 0 , 0 ) for the s=0𝑠0s=0italic_s = 0 (gray) phase, (0.8,0)0.80(0.8,0)( 0.8 , 0 ) for the s=1𝑠1s=1italic_s = 1 (blue) phase and (0,0.8)00.8(0,0.8)( 0 , 0.8 ) for the s=2𝑠2s=2italic_s = 2 (red) phase. The stationary densities (ρ1,ρ2)subscript𝜌1subscript𝜌2(\rho_{1},\rho_{2})( italic_ρ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_ρ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) of the three phases are indicated as the green points in Fig. 14(a), matching the minima of the free energy f𝑓fitalic_f. Such travelling waves stochastically appear also in two dimensions, as can be seen in Fig. 3(a) and Movie S1. They also produce the spiralling dynamics around the contact point of three different domains (see Figs. 3(a) and 10(d)).

Biphasic bands travelling in opposite directions annhilate each other, as shown in Fig. 14(c), even when their sizes are different, as a natural consequence of the s=0120𝑠0120s=0\to 1\to 2\to 0italic_s = 0 → 1 → 2 → 0 transformation. Such events can also be seen in two dimensions (see Fig. 10(d) and Movie S1).

In the case of asymmetric flipping energies, we find that three interfaces s=01𝑠01s=0\to 1italic_s = 0 → 1, s=12𝑠12s=1\to 2italic_s = 1 → 2 and s=20𝑠20s=2\to 0italic_s = 2 → 0 have different velocities. Depending on the relative asymmetries, this can lead to the widening of a biphasic band, as shown in Fig. 15(a), or its disappearance (see Fig. 15(b)).

VII Summary

We have studied the dynamics of the active cyclic Potts model under asymmetric conditions. We found that biphasic domains and non-cyclic one-state dominant phases are induced by either asymmetric flipping energies or asymmetric contact energies. The spiral wave mode can temporally coexist with either the one-state dominant phases or the homogeneous cycling mode. Biphasic domains move like amoeba and exhibit division and disappearance. When the flipping energy from the s=0𝑠0s=0italic_s = 0 to the s=1𝑠1s=1italic_s = 1 state is increased, or the contact energy between s=0𝑠0s=0italic_s = 0 sites is decreased while other flipping energies are fixed, the s=2𝑠2s=2italic_s = 2 dominant phase appears owing to the suppression of the nucleation of s=0𝑠0s=0italic_s = 0 domains in the s=2𝑠2s=2italic_s = 2 phase. Two separate modes can be formed depending on the initial state due to hysteresis: one-state dominant phase and coexistence of spiral waves and one-state dominant phase, or two types of one-state dominant phases.

In reaction-diffusion systems, the length scales of spatiotemporal patterns are usually controlled by the diffusion coefficients of the reactants Mikhailov and Ertl (2009). However, in the present system, the diffusion of the various states is not taken into account. The observed wavelengths fluctuate largely and are determined, on average, by the competition between domain nucleation and growth. We mentioned three experimental situations relevant to the present cyclic model in the Introduction, including reaction on a catalytic surface. Although oxidation or reduction on catalytic surfaces has been intensively studied in experiments Ertl (2008); Brär et al. (1994); Gorodetskii et al. (1994); Barroo et al. (2020); Zeininger and et al. (2022), the effects of nucleation and growth have not been well understood. We believe that the rich dynamics of the present system, including the amoeba-like motion of biphasic domains, can be experimentally observed on catalytic surfaces. Here, we simulated small systems. Such small systems can be realized using metal nanoparticles Tang et al. (2020). Chemical waves are known to induce shape deformations of metal particles Tang et al. (2020); Ghosh et al. (2022), lipid membranes Wu et al. (2018); Tamemoto and Noguchi (2021); Noguch (2023), and macroscopic gel sheets Ionov (2014); Maeda et al. (2008); Levin et al. (2020). In this study, we have only considered non-deformable flat surfaces. Surface deformation and the dynamics on curved surfaces are likely to modify the dynamics of the active cyclic Potts model and yield new phenomena.

Here, we used the three-state Potts model in the square lattice. It is applicable systematically from thermal equilibrium to far-from-equilibrium. The dynamics may be changed when other types of lattices are used. When four or more states are considered, more complex dynamics are expected. Hence, the cyclic Potts model is an excellent model system to study nonequilibrium dynamics coupled with nucleation and growth. It is simple and easy to extend to several directions for further studies.

Acknowledgements.
We thank Frédéric van Wijland for stimulating discussion. This work was supported by JSPS KAKENHI Grant Number JP21K03481 and JP24K06973.

References