Chi-Jen Wang, Da-Jiang Liu, Jim W. Evans

​

AMC2025, 3 Aug 2025

Triplet approximation in Schloegl’s second model for autocatalysis on the square lattice

DISCONTINUOUS PHASE-TRANSITIONS IN

NON-EQUILIBRUM SYSTEMS

Background: Discontinuous Transitions in Equilibrium Systems

P

pressure

n

non-ideal gas

Liquid = L

Gas = G

density

Van der Waals equation of state

(with van der waals loop)

P

pressure

n

non-ideal gas

Liquid = L

Gas = G

density

Metastable Gas = MG

ML

Maxwell

construction

Actual

behavior

P_e

Discontinuous phase transition

at “equistability” pressure P_e

gas

liquid

Coexistence at p=pe

Stationary planar interface

CONDENSATION P>P_e

L

L

MG

Super-

critical

liquid

droplet

grows

Sub-critical

shrinks

SCHLOEGL’S FIRST MODEL

SCHLOEGL’S SECOND MODEL

Phase Diagrams

SCHLOEGL’s models for autocatalytic chemical reactions

SCHLOEGL’S SECOND MODEL FOR AUTOCATALYSIS

QUADRATIC CONTACT PROCESS FOR SPATIAL EPIDEMICS

Concentrations of X� C

p�annihilation rate

Population Density of H�� H

p

recovery rate

MODELS FOR CATALYSIS AND EPIDEMICS

NUCLEATION & GROWTH OF DROPLETS

exchange rate h=0

​

P =0.098 t =455, 2275, 4090

​

P_eq = 0.094 P_s = 0.101

P

p_eq p_s

S

H

h=1

​

P=0.208

t= 408, 960,

1440,…

​

P_eq =0.197

​

P_s = 0.21

P

p_eq p_s

Guo, De Decker, Evans, Phys. Rev. E 82 (2010) 021121

NUCLEATION & GROWTH OF DROPLETS

S

vacuum

droplet

PROPRAGATION VELOCITY FOR VARIOUS SLOPE OF INTERFACES h>0

h=0

h=0.001

h=0.1

wave solution

Interface slope=1, size 256X256

droplet solution

1. h=0, 0.001, and 0.1; p<peq
2. h=0, 0.001, and 0.1; p>peq

�Size ~10⁴ sites

critical droplets

KMC SIMULATIONS for small h=0.001

h>0 removes the extreme V=0 behavior

• small p, infected state can displace the healthy state
• large p, healthy state can displace the infected state

t=0 8000 24000 40000

p= 0.092

​

​

p= 0.095

​

​

p= 0.095

​

​

p= 0.0097

< 2PC region <

healthy droplets

@ p = 0.214

R > R_c grow

R < R_c shrink

R = R_c critical

grow

shrink

grow

shrink

infected droplets

@ p = 0.206

R > R_c⁺ grow

R < R_c^- shrink

Generic two-phase

Coexistence

(2PC)

​

​

​

​

​

Both droplets

shrink

p

R

Droplet

radius

Droplet area

A ≡ π R²

stationary

DROPLET DYNAMICS

Durrett �h=0.01

Dimension dependence

MF approximation

PRE 2012

Set P[o_i] = Cm , P[∙_i] = 1 – Cm, where

​

EXTENSIONS OF THE QUADRATIC CONTACT PROCESS

Possible Refinements

Liu et al. Phys. Rev. Lett. 98 (2007) 050601; Guo et al. PRE 75 (2007); Physica A 387 (2008)

Liu, J. Stat. Phys. 135 (2009); Guo, Liu and Evans, J. Chem. Phys. 130 (2009) 074106

Guo et al. PRE 82 (2010); Physica A 391 (2012)

Wang et al., J. Stat. Phys. 144 (2011); Phys. Rev. E (2012); J. Chem. Phys. (2015); Phys. Rev. E (2020);

Cubic lattice realization

Spontaneous annihilation

at rate P=1/τ_filled

creation by diagonal neighbor pairs

at rates k/12 where k = # neighboring diagonal ● pairs

p

EXTENSIONS OF THE QUADRATIC CONTACT PROCESS

Possible Rate Change

0 ⅙ ½ 1

p

Phys. Rev. E (2023)

EXTENSIONS OF THE QUADRATIC CONTACT PROCESS

Bethe lattice Realization

spontaneous �annihilation

at rate p=1/τ_filled

creation by neighbor pairs

at rates k/3 where �k = # neighboring ● pairs

p

1/3

1

Labeling the nodes j on Bethe lattice

Labeling the rings k

�

i i i

Radius symmetry cases�set C_k = P∙_k for the “concentration” of particles in ring k.

1. annihilation rate p at filled site i

2. particle creation at empty site i

P∙_k + Po_k, = 1

Liu Wang Evans �Phys. Rev. E. 104, 014135 (2021)

Creation rate: r_n=_nC₂/_zC₂ =n(n-1)/z(z-1) , �for n occupied neighbors out of z neighbors

Bethe lattice

an infinite connected cycle-free graph where each node is connected to z neighbors

z=3

z=4

finite graph�z=3

infinite graph

Related Works – finite graph

Kouvaris, Sebek, Mikhailov, Kiss 2016

u_i is the chemical concentration, f(u) = −u(u − h)(u − a) with 0<h<a,

K is the strength of diffusive coupling,�Aij=1 if there is a edge connects node i and node j, Aij=0 otherwise.

​

d/dt C_k = -p C_k + (1-C_k) C_k+1 [2 C_k-1 + (z-2) C_k+1]/z, for k ≥1,

MF version lattice differential equations (discrete reaction-diffusion equations) (dRDE)

p

1/3

1

Mean field approximation

d/dt C₀ = -pC₀ + (1-C₀)(C₁)² for k = 0.

C_±(MF) = ½ ± ½ (1-4p)^1/2 for 0 ≤ p ≤ p_s(MF)= ¼, and C_vac = 0.

Here C₊ and C_vac =0 are stable, and C_- is an unstable steady state. �C_± are a populated states, and C_vac is the trivial absorbing vacuum state. C₊ and C_- disappear beyond the spinodal point p_s, a sn-bifurcation.

MF approximation : � P∙_k-1–o_k–∙_k+1 ≈ P∙_k-1 Po_k P∙_k+1⁼ C_k-1(1-C_k)C_k+1

Visualizaton for z=3

Simple truncation

Pair approximation

Pair approximation : P∙_k-1–o_k–∙_k+1

≈ (P∙_k-1–o_k) (Po_k–∙_k+1 ) / Po_k

​

p

p

having correlations �with neighbors

k k k+1

k-2 k-1 k-1

k

k

k-1 k-1

k

(a) d/dt P∙_k = - p P∙_k + z^-1[2 P∙_k-1o_k Po_k∙_k+1 + (z-2)(Po_k∙_k+1)²]/Po_k 

�(b) d/dt Po_k-1o_k = +p P∙_k-1o_k + p Po_k-1∙_k

- (z-2)z^-1(z-1)^-1[2P∙_k-2o_k-1 Po_k-1∙_k + (z-3) (Po_k-1∙_k)²]� Po_k-1o_k /(Po_k-1)²  

- (z-2)z^-1 (Po_k∙_k+1)² Po_k-1o_k /(Po_k)²,

Spatially homogeneous case - pair

Set K = P∙_k-1–o_k / Po_k = Po_k–∙_k+1 / Po_k = P∙o/Po� for z=3, d/dt Poo = 2 P∙o [p – (z-2)K(1-K)/z] becomes

d/dt (C-K(1-C)) = 2 K(1-C) [p - K(1-K)/3] .   

​

​

​

​

​

​

stable populated steady state with 

K_±(pair) = ½ ± ½ [1 - 4zp/(z-2)]^1/2

C_±(pair) = (K_±)²/[p + (K_±)²] = 1/{1+1/[K_±/p - z/(z-2)]}  

= 1 - ½ {-1 + 4p/(z-2) ± [1 - 4zp/(z-2)]^1/2}/[1+4p/(z-2)²],

for 0 ≤ p ≤ p_s(pair) = (z-2)/(4z), the spinodal point. ��C₊ and C_- correspond to stable and unstable populated steady states, �as do K₊ and K_-. Vacuum steady state C_vac(pair) = K_vac(pair) = 0.

pair-approximation for homogeneous states yields

  d/dt C = -pC + (1-C)K²

d/dt (C-K(1-C)) = 2 K(1-C) [p – (z-2)K(1-K)/z]

Triplet-like approximation

Triplet-like approximation : Po_k-2∙_k-1o_k∙_k+1o_k+2 ≈ Po_k-2∙_k-1o_k Po_k∙_k+1o_k+2 / Po_k

P∙_k-1o_k∙_k+1o_k+2 ≈ P∙_k-1o_k∙_k+1 Po_k∙_k+1o_k+2 / Po_k∙_k+1

Higher-order truncation

Rewrite

(b) d/dt Po_k-1o_k = +p P∙_k-1o_k + p Po_k-1∙_k

 - (z-2)z^-1(z-1)^-1 [2P∙_k-2o_k-1o_k Po_k-1o_k∙_k + (z-3)(Po_k-1o_k∙_k)²] /Po_k-1o_k

- (z-2)z^-1 (Po_k-1o_k∙_k+1)² /Po_k-1o_k

(c) d/dt Po_k-1o_ko_k+1 = +p (P∙_k-1o_ko_k+1 + Po_k-1∙_ko_k+1 + Po_k-1o_k∙_k+1)

- (z-2)z^-1(z-1)^-1 (2P∙_k-2o_k-1∙_k + (z-3)Po_k-1∙_k∙_k) Po_k-1o_ko_k+1 /Po_k-1

- (z-2)(z-3)z^-1(z-1)^-1 Po_k∙_k+1∙_k+1 Po_k-1o_ko_k+1 /Po_k

- (z-2)z^-1 Po_k+1∙_k+2∙_k+2 Po_k-1o_ko_k+1 /Po_k+1, for k ≥ 2,

(d) d/dt Po_ko_k+1o_k+1 = +p (P∙_ko_k+1o_k+1 + 2Po_ko_k+1∙_k+1)

- (z-3)z^-1(z-1)^-1 (2P∙_k-1o_k∙_k+1 + (z-4)Po_k∙_k+1∙_k+1) Po_ko_k+1o_k+1 /Po_k

- 2(z-2)z^-1 Po_k+1∙_k+2∙_k+2 Po_ko_k+1o_k+1 /Po_k+1, for k ≥ 1,

(e) d/dt Po_k-1∙_ko_k+1 = +p (P∙_k-1∙_ko_k+1 - Po_k-1∙_ko_k+1 + Po_k-1∙_k∙_k+1)

- 2z^-1(z-1)^-1 (P∙_k-2o_k-1∙_k + (z-2)Po_k-1∙_k∙_k) Po_k-1∙_ko_k+1 /Po_k-1∙_k

+ (z-2)(z-3)z^-1(z-1)^-1 Po_k∙_k+1∙_k+1 Po_k-1o_ko_k+1 /Po_k

- 2z^-1 P∙_ko_k+1∙_k+2 Po_k-1∙_ko_k+1 /P∙_ko_k+1, for k ≥ 2,

(f) d/dt P∙_ko_k+1o_k+1 = +p (- P∙_ko_k+1o_k+1 + 2P∙_k∙_k+1o_k+1)

+ (z-3)z^-1(z-1)^-1 (2P∙_k-1o_k∙_k+1 + (z-4)Po_k∙_k+1∙_k+1)Po_ko_k+1o_k+1 /Po_k

- 4z^-1 P∙_ko_k+1∙_k+2 P∙_ko_k+1o_k+1 /P∙_ko_k+1, for k ≥ 1.

Triplet approximation

Higher-order truncation

KMC simulation results

global steady-state concentration C versus p for Schloegl’s second model on a finite Bethe lattice with z = 3 and BC at R₁₇. Each data point is obtained as an average over 10⁴ Monte Carlo Steps (MCS).

​

dP(S)/dt = - pP(S) + P(H)P(S)² = - pP(S) + [1- P(S)] P(S)² = R(S)

S→H H+2S→3S � loss gain

P(S) = fraction of sick individuals

P(H) = fraction healthy = 1 – P(S)

KINETICS AND STEADY-STATES

BISTABILITY

stable infected

stable all healthy

P_S =1/4

sn bifurcation

or “spinodal”

P(S)

​

unstable

steady state

P(S)=[1+√(1-4p)] /2

Steady-States

HYDRODYNAMIC LIMIT h>>1 rapid exchange/ strong mixing

recovery rate

A.S. Mikhailov, Foundations of Synergetics I (1991)

Equistability p_eq=2/9=0.222 �(stationary interface)

1

​

​

​

​

​

​

​

​

S

​

​

​

​

​

​

​

​

​

0

x

all healthy

infected

V>0 for p<p_eq

wave of epidemic spread

V<0 for p>p_eq

TRIGGER WAVES & EQUISTABILITY

∂S/∂t = R(S) + h ∇²S

−

HYDRODYNAMIC LIMIT h>>1 rapid exchange/ strong mixing

2D KINETIC MONTE CARLO (KMC) SIMULATION h=0

Liu, Guo, Evans, PRL 98 (2007) 050601

Extreme case: h=0 (no exchange mobility)

all healthy state

Infected state

2D KINETIC MONTE CARLO (KMC) SIMULATION h=0

Liu, Guo, Evans, PRL 98 (2007) 050601

Extreme case: h=0 (no exchange mobility)

stable infected state exists up to p_e

p_eq for stationary interface now

depends on orientation

EXACT MASTER EQUATION for HOMOGENEOUS STATES: square lattice

S

d/dt P[S] = – p P[S] + P[ H S ] �

S

To deal with the term P[ H S ] :

P[S] : probability that individual at the site is sick.

​

P[H] = 1 – P[S]: probability that individual at the site is healthy.

​

P[H S] : probability that healthy and sick individuals are neighbors.

master equation

…example for h=0

S

mean field approximation

d/dt P[S] = – p P[S] + P[H] P[S] P[S]

d/dt P[S] = – p P[S] + P[H S ] P[H S] /P[H]

d/dt P[HS] = much more complicated

pair approx.

Triplet Approximation

dRDE vs KMC PREDICTION: square lattice h=0

Corresp to �KMC 0.09440

Corresp to �KMC 0.0871

V

p

0.21138

diagonal 11

Near-vertical

0.20505

horiz./vert. 10

0.20711

0.10831

0.10602

Near-vertical

0.10560

p_e = 0.09440

p_f= 0.0871

V

p

V

p

mean-

field

(site)

pair

KMC

0

0

0

H

H

H

H

S

S

S

S

H

H

H

H

H

S

S

S

S

H

S

S

S

S

H

H

PRE (2012)

Near-vertical

Higher order approximation: Triplet

Kirkwood superposition approximation method

Higher order approximation: Triplet

• Kirkwood superposition approximation method, P(x1,x2,x3)

Higher order approximation: Triplet

Another example

Higher order approximation: Triplet

… more coupled equations

Homogeneous Hierarchical Master Equations

Site independent P[o]= P[o_i] for all i

Homogeneous Hierarchical Master Equations

Simplification Setting

Homogeneous Hierarchical Master Equations

The steday stats of homogeneous triplet approximation equation occurs

when the time derivative is zero. Therefore, the steady states satisfy

a) pC=((1–Q) –(1–L)Q)(1–C)

b) p2(1–Q)³ = ((1–Q) –(1–L)Q)(1–L)Q(1–K)

c) p(2(1–K)Q+ R(1–Q))(1–Q)² = KQ(1–L)(1–K)((1–Q) –(1–L)Q),

d) 4p(2(1–L)Q+S(1–Q)) (1–Q)⁴ ((1–L)Q+(1–S)(1–Q)) (LQ+S(1–Q)) =

QL(1–K)((1–Q) –(1–L)Q)

[ 2CS(1–Q) (1–L)² + ( ((1–L)Q+(1–S)(1–Q)) (LQ+S(1–Q)) ) (2(1–L)(1–Q)² + (1–K)³ ) ]

e) p(2–3R)(1–Q)⁴ = R((1–Q) –(1–L)Q)[ ((1–Q) –(1–K)Q)((1–Q) –(1–L)Q) + (1–Q)³ ],

(C–((1–L)Q+(1–S)(1–Q) )(1–C)) ((LQ+S(1–Q))) (1–Q)³

f) [ 4p(2–3S) (1–Q)³ + LQ (1–L)² ((1–Q) –(1–L)Q) KQ (1–K)Q (1–C)⁴ ]

=

2 S² C² ((1–Q) –(1–K)Q) ((1–Q) –(1–L)Q) [ (1–Q)³ + ((1–Q) –(1–L)Q)² ]

+2(C–((1–L)Q+(1–S)(1–Q) )(1–C)) ((LQ+S(1–Q))) (1–Q)² ((1–Q) –(1–L)Q) S

[ (1–Q)³ + C [((1–Q) –(1–K)Q)] [((1–Q) –(1–L)Q)]² ]

Working

for h≥0

Perturbed Durrett models: bistability of steady-states

p

infected steady state

all healthy state

bistability

bistability

unstable steady state

MF approx. with exchanging

S

SQUARE LATTICE: MF AND PAIR STEADY-STATES

unstable steady state

all healthy state

0.0 0.05 0.1 0.15

p

Pair approx. with exchanging

infected steady state

EXACT MASTER EQUATION for INHOMOGENEOUS STATES: square lattice

S_i,j+1 S_i,j+1

d/dt P[S_i,j] = – p P[S_i,j] + (1/4){P[ H_i,j S_i+1,j] + P[S_i-1,j H_i,j ] + P[H_i,j S_i+1,j ] + P[S_i-1,j H_i,j] } S_i,j-1 S_i,j-1

S_i,j+1

Need to deal with the term P[ H_i,j S_i+1,j ]

Considering P[S_i,j] , P[H_i,j], P[H_i,j S_i+1,j]

master equation

…example for h=0

(i, j)

mean field approximation

pair approx.

Wang, Liu, Evans, Phys. Rev. E 85 (2012) 041109

S_i,j+1

P[ H_i,j S_i+1,j ] ≈ P[H_i,j] P[S_i+1,j] P[S_i,j+1]

S_i,j+1

P[ H_i,j S_i+1,j ] ≈ P[H_i,j S_i+1,j] P[H_i,j S_i,j+1] / P[H_i,j]

EXTENSIONS OF THE QUADRATIC CONTACT PROCESS

Nonlocal

Phys. Rev. L (2018)

• Why different shapes?

KMC SIMULATIONS for small h=0.01

• small p, infected droplet, diagonal interface moves faster than other interfaces.
• large p, healthy droplet, vertical interface moves faster than other interfaces.

diagonal

vertical

V

p

0

healthy droplets

@ p = 0.214

R > R_c grow

R < R_c shrink

R = R_c critical

grow

shrink

grow

shrink

infected droplets

@ p = 0.206

R > R_c⁺ grow

R < R_c^- shrink

Generic two-phase

Coexistence

(2PC)

​

​

​

​

​

Both droplets

shrink

p

R

Droplet

radius

Droplet area

A ≡ π R²

stationary

DROPLET DYNAMICS

Durrett �h=0.01

R

p

Generic

Two-phase

coexistence

grow

grow

shrink

shrink

healthy

& infected

droplets

shrink

@ p = 0.211

Durrett with h=0.01 (slow hopping)

Droplet area

A ≡ π R²

R

t

Families of distinct

stationary droplets

with 4-fold symmetry

Durrett with h=0.01 (slow hopping)

Another families of distinct

stationary droplets with

a rectangle like shape

p

R

Phys. Rev. E (2020)

Families of distinct

stationary droplets

with 4-fold symmetry

Durrett with h=0.01 (slow hopping)

R

p

R is too large, hard to distinct �the droplet families

R

t

p=0.21274, h=0.01

0.3% difference in radius of these 4 droplets

Durrett with h=0.01 (slow hopping)

1/R

p

linear vanishing behavior of curvature κ_c = 1/R_c 1/R_c ~ δp=|p-p_eq| for p not too close to p_eq

p(~diag)&p₊(diag) are similiar

|p-p_eq|/p_eq )^0.5

|p-p_eq|/p_eq )^0.22

nonlinear

nonlinear

|p-p_eq|/p_eq )^0.26

Durrett with ε=0.001 (weak spontaneous creation)

Generic two-phase

Coexistence

(2PC)

​

​

​

​

​

Both droplets

shrink

p

R

Droplet

radius

Droplet area

A ≡ π R²

stationary

Relations between wave and droplet solutions

Durrett �h=0.01

Durrett with h=0.01 (slow hopping)

Critical Curvature

|κ_c(S=∞) | ≈ 0.08023 + 21.6|V_p(S=∞)| �in the narrow range of V_p(S=∞) analyzed here, �

|κ_c(S=1)| ≈ 0.05558 + 38.85|V_p(S=1)| �at least provided that |V_p(S=1)| is not too small.

κ_c(S=1) ≈ f(V_p(S=1)) vanishes �as V_p(S=1)→0

Expect

Families of distinct

stationary droplets

with 4-fold symmetry

Liu, Guo, Evans, Phys. Rev. Lett. 98, 050601 (2007)

Guo, Liu, Evans, J. Chem. Phys. 130, 074106 (2009)

CJ Wang, Liu, Evans, Phys. Rev. E 85, 041109 (2012)

CJ Wang, Liu, Evans, J. Chem. Phys., 142, 164105 (2015)

Liu , CJ Wang, Evans, Phys. Rev. Lett. 121, 120603 (2018)�CJ Wang , Liu, Evans, Phys. Rev. E (2020)�Liu, CJ Wang , Evans, Phys. Rev. E (2021)�Shen , Liu, Evans, Phys. Rev. E (2023)

Thanks for �your attention