1 of 32

Chi-Jen Wang�Da-Jiang Liu, James W. Evans

NSTC, 26 Nov 2025

Triplet approximation in Schloegl’s second model for non-equilibrium phase transitions in lattice bistable systems

Schloegl第二模型於晶格雙穩態系統中非平衡相變的三階近似

p

1/3

1

2 of 32

NSTC Grand 111-2628-M-194-001-MY3Schloegl第二模型於空間流行病模型中的高階近似

…共15組方程式

TODAY�…共6組方程式

  • 研究主題 1

二維網格的三階近似

  • 研究主題 2

Bethe 網格的三階近似

  • 研究主題3

二維網格的近似滴液解

3 of 32

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

Pe

Discontinuous phase transition

at “equistability” pressure Pe

gas

liquid

Coexistence at p=pe

Stationary planar interface

P

L

G

MG

ML

Pe

n

CONDENSATION P>Pe

L

L

MG

Super-

critical

liquid

droplet

grows

Sub-critical

shrinks

4 of 32

Bethe lattice

Bethe lattice is an infinite connected cycle-free graph where

each node is connected to z neighbors

z=3

z=4

finite graph�z=3

infinite graph

5 of 32

Related Works – finite graph

Kouvaris, Sebek, Mikhailov, Kiss 2016

ui 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.

6 of 32

Related Works

Kouvaris, Sebek, Mikhailov, Kiss 2016

7 of 32

Related Works – 1st order phase trans

Chatterjeea and Durrett, 2013

  • If θ ≥ 2, r ≥ θ + 2, ϵ1 is small and p ≥ p1(ϵ1), then starting from all vertices occupied the fraction of occupied vertices is ≥1 − 2ϵ1 up to time exp(γ1(r)n) with high probability.
  • The process on the r-tree has a first-order phase transition.

Consider the discrete time threshold-θ contact process on a random r-regular graph and on trees in which all vertices have degree r.

Sites having at least θ many occupied neighbors at time t become occupied at time t + 1 with probability p.

8 of 32

Schloegl’s second model (quadratic contact process) on a lattice involves �1. spontaneous particle annihilation at rate p, and

2. autocatalytic particle creation at empty

sites with two or more occupied neighbors.

  • Various boundary conditions and unconventional simulation ensembles on the Bethe lattice to extract behavior for infinite size.
  • A discontinuous transition to the vacuum state when p exceeds a critical value.

Models of the work

p

1/3

1

9 of 32

p

1/3

1

Models of the work on Bethe lattice(z=3)

Labeling the nodes j on Bethe lattice

Labeling the rings k

Quadratic contact process

i i i

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

1. annihilation rate p at filled site i

2. particle creation rate rn at empty site i

P∙k + Pok = 1

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

Creation rate: rn=nC2/zC2 =n(n-1)/z(z-1) , �for n occupied neighbors out of z neighbors

10 of 32

p

1/3

1/3

Mean Field (MF) approximation

Visualizaton for z=3

Radius symmetry assumption

Master Equation:

P∙k-1–ok–∙k+1 ≈ P∙k-1 Pok P∙k+1 = Ck-1(1-Ck)Ck+1

k+1 � Pok–∙k+1 ≈ Pok P∙k+1 P∙k+1 = (1-Ck)Ck+1Ck+1

Ck = P∙k for the “concentration” of particles in ring k.

Mean Field

Approximation

11 of 32

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

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

p

Mean Field (MF) approximation

d/dt C0 = -p C0 + (1-C0)(C1)2 for ring k = 0.

C±(MF) = ½ ± ½ (1-4p)1/2 for 0 ≤ p ≤ ps(MF)= ¼, and Cvac = 0.

Here C+ and Cvac =0 are stable, and C- is an unstable steady state. �C+ and C- disappear beyond the spinodal point ps, a sn-bifurcation.

P∙k-1–ok–∙k+1 ≈ P∙k-1 Pok P∙k+1 = Ck-1(1-Ck)Ck+1

k+1 � Pok–∙k+1 ≈ Pok P∙k+1 P∙k+1 = (1-Ck)Ck+1Ck+1

Radius symmetry assumption

Ck = P∙k for the “concentration” of particles in ring k.

1/3

12 of 32

dC/dt = - pC + (1-C)C2� loss gain

Homogeneous assumption Ck = C for all sites.

KINETICS AND STEADY-STATES

BISTABILITY

populated stable

vacuum stable

PS =1/4

sn bifurcation

or “spinodal”

C

unstable

steady state

active stable =C+(MF)

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

Steady-States

Spatial homogeneous case - MF

annihilation rate

ubstable = C-(MF) � =[1-√(1-4p)] /2

13 of 32

Radius Symmetric (Non-homogeneous)

- MF estimate of the propagation velocity

(Left) the populated steady state embedded in the vacuum state;

(Right) the vacuum state embedded in the populated steady state.

V>0

V<0

V>0

V<0

14 of 32

MF estimate of the propagation velocity

V>0

V<0

15 of 32

Pair approximation

Pair approximation : P∙k-1–ok–∙k+1

≈ (Pk-1–ok) (Pokk+1 ) / Pok

p

p

with neighboring correlations

k k

k-1 k-1

k k

k-2 k-1 k-1 k+2

k

k+2

not ODE yet

16 of 32

Spatially homogeneous – pair approximation

 

   

pair-approximation for homogeneous states

  d/dt P∙ = -pP∙ + (P∙o)2/ Po

d/dt Poo = 2 P∙o [p – (z-2)P∙oPoo / z(Po)2 ]

Assumption Pk = Ck = C for all sites.

Pair approximation : P∙–o–∙ ≈ (P∙–o) (Po–∙) / Po

17 of 32

Simplification

Set K = P∙k-1–ok / Pok = Pok–∙k+1 / Pok = P∙o/Po

(a) d/dt P∙ = -pP∙ + (P∙o)2/ Po becomes � d/dt C = -pC + (1-C)K2

� (b) d/dt Poo = 2 P∙o [p – (z-2)P∙oPoo / z(Po)2 ]

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

   

homogeneous pair-approximation for z=3

  (a) d/dt C = -pC + (1-C)K2

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

pair-approximation for homogeneous states

  (a) d/dt P∙ = -pP∙ + (P∙o)2/ Po

(b) d/dt Poo = 2 P∙o [p – (z-2)P∙oPoo / z(Po)2 ]

18 of 32

Spatially homogeneous case - pair

   

stable populated steady state with 

K±(pair) = ½ ± ½ [1 - 12p]1/2

C±(pair) = (K±)2/[p + (K±)2] = 1/{1+1/[K±/p - 3]}  

= 1 - ½ {-1 + 4p ± [1 - 12p]1/2}/[1+4p],

for 0 ≤ p ≤ ps(pair) = 1/12≈0.08333, the spinodal point. ��C+ and C- correspond to� stable and unstable steady states.

Vacuum steady state

Cvac(pair) = Kvac(pair) = 0.

homogeneous pair-approximation for z=3

  (a) d/dt C = -pC + (1-C)K2

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

populated stable

unstable

steady state

vacuum stable

PS =1/12

sn bifurcation

or “spinodal”

BISTABILITY

Steady-States

C

annihilation rate p

19 of 32

Pair estimate of the propagation velocity

20 of 32

p

1/3

1/3

Triplet approximation

for z=3

Radius symmetry assumption

Master Equation:

p

p

k k

k-1 k-1

k

k-1 k+2

k+2

21 of 32

Conti – Master Equation

for z=3

k+1

k-1 k+1

k

k+2

k

k+3

k+1

k+2

k+3

change rate of

P(ok-1–ok–ok+1)

change rate of

ok+1

P(ok–ok+1)

22 of 32

Conti – Master Equation

for z=3

P(ok-1–∙k–ok+1) can’t be expressed as a linear combination of 1, P(oi), P(oi oj), P(oi oj ol).

Unable to enclose (a),(b),(c),(d).

Need the changing rate description of the probability density P(ok-1–∙k–ok+1).

6 ODEs.

23 of 32

Kirkwood approximation

Kirkwood approximation [Kirkwood 1935 JCP]

 

 

 

24 of 32

Triplet approximation

 

 

Apply the triplet approximation to the master equations

(a)-(f).

25 of 32

Simplifications

 

Triplet approximation in Bethe lattice

…3 more eqs

26 of 32

Triplet approximation in Bethe lattice

 

27 of 32

Homogeneous (site independent)

Spatial homogeneous case - Triplet

3 more eqs for dF/dt, dG/dt, dH/dt.

28 of 32

Steady Staes - Homegeneous

Difficult to find algebraic form of steady states

Numerical� Results

29 of 32

Compare with KMC simulation

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

30 of 32

Boundary conditions

How to characterize behavior on an infinite Bethe lattice?

annihilations and creations on the first k* rings 1≤k≤k*.

boundary conditions (BCs) imposed on ring k*+1

conventional �simulation �ensemble

period BCs in infinite Bethe lattice?

period BCs�

31 of 32

choices (K1, K2, K3) of BCs

KMC simulation results

(K1) sets K = K1 = <C>/(3-2<C>),

(K2) sets K = K2 = [p<C>/(1-<C>)]1/2

motivated by the pair approximation

(K3) sets K = K3, exactly from simulated configurations.

K = P∙o/Po

32 of 32

  • a discontinuous transition for Schloegl’s second model on an infinite Bethe lattice with z = 3 occurring at around p ≈ 0.06.

  • choice of BCs (Kj) for the finite lattice which accounts for strong spatial correlations in the model can best mimic behavior for infinite lattice.

  • expect a discontinuous transition also for the cases z>3.

  • behavior for bigger z reflects more closely the mean-field prediction of bistability, but the presence of noise in the stochastic model ensures the presence of a discontinuous transition.

SUMMARY