Chi-Jen Wang�National Chung Cheng University

wiht Da-Jiang Liu, Jim W. Evans�Ames Lab, Iowa State University

​

WUSTL, July 23, 2024

Discontinuous non-equilibrium phase transition in Schloegl’s second model for autocatalysis

0 0 ¼ ½ 1

p

Stochastic cooperative processes �

dynamics and pattern formation in systems involving

chemical reactions, �population dynamics,

spread of epidemics or information

continuous/discontinuous phase transitions

in Hamiltonian systems

Metastability, interface propagation and nucleation

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=p_e

Stationary planar interface

Condensation P>P_e

L

L

MG

Super-

critical

liquid

droplet

grows

Sub-critical

shrinks

SCHLOEGL’s models for autocatalytic chemical reactions

aka Quadratic Contact Process

SCHLOEGL’S SECOND MODEL

Phase Diagrams

SCHLOEGL’S FIRST MODEL

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

​

dP(X)/dt = - pP(X) + P(∅)P(X)² = - pP(X) + [1- P(X)] P(X)² = R(X)

X→∅ ∅+2X→3X � loss gain

P(X) = concentration of X

P(∅) = concentration of ∅ = 1 – P(X)

KINETICS AND STEADY-STATES

BISTABILITY

stable actived

stable vacuum

P_S =1/4

sn bifurcation

or “spinodal”

C=P(X)

​

unstable

steady state

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

Steady-States

HYDRODYNAMIC LIMIT h>>1 rapid exchange/ strong mixing

annihilation rate

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

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

1

​

​

​

​

​

​

​

​

C

​

​

​

​

​

​

​

​

​

0

x

vacuum

actived

V>0 for p<p_eq

V<0 for p>p_eq

TRIGGER WAVES & EQUISTABILITY

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

−

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)

Vacuum state

active state

Vacuum

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

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

Extreme case: h=0 (no exchange mobility)

stable active state exists up to p_e

p_eq for stationary interface now

depends on orientation

vacuum

vacuum

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

C

1-C

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

C

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

t=0 8000 24000 40000

p= 0.092

​

​

p= 0.095

​

​

p= 0.095

​

​

p= 0.0097

< 2PC region <

h>0 removes the extreme V=0 behavior

• small p, active state can displace the vacuum state
• large p, vacuum state can displace the active state

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

Durrett �h=0.01

DROPLET DYNAMICS

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

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.

​

Kouvaris, Sebek, Mikhailov, Kiss (2016)

EXACT MASTER EQUATION for INHOMOGENEOUS STATES: square lattice

​

●_i,j+1

Need to deal with the term P[ ○_i,j ●_i+1,j ]

master equation

…example for h=0

Considering P[●_i,j] , P[○_i,j] , P[○_i,j ●_i+1,j]

​

​

(i, j)

●_i,j+1 ●_i,j+1

d/dt P[○_i,j] = p P[●_i,j] – (1/4){P[ ○_i,j ●_i+1,j ] + P[ ● _i-1,j ○_i,j ] + P[ ○_i,j ●_i+1,j ] + P[ ● _i-1,j ○_i,j ] }

●_i,j-1 ●_i,j-1

mean field approximation

pair approx.

●_i,j+1

P[ ○_i,j ●_i+1,j ] ≈ ��P[○_i,j] P[●_i+1,j] P[●_i,j+1]

●_i,j+1

P[ ○_i,j ●_i+1,j ] ≈

�P[○_i,j ●_i+1,j] P[○_i,j ●_i,j+1] / P[○_i,j]

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

dRDE for higher d is analogous

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

Current Work

… more coupled equations

for h≥0

Perturbed Durrett models: bistability of steady-states

p

activesteady state

Vacuum state

bistability

bistability

unstable steady state

MF approx. with exchanging

S

SQUARE LATTICE: MF AND PAIR STEADY-STATES

unstable steady state

Vacuum state

0.0 0.05 0.1 0.15

p

Pair approx. with exchanging

Active steady state

EXTENSIONS OF THE QUADRATIC CONTACT PROCESS

Nonlocal

Phys. Rev. L (2018)

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

Durrett �h=0.01

DROPLET DYNAMICS

R

p

Generic

Two-phase

coexistence

grow

grow

shrink

shrink

active & vacuum

droplets

all 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

• Why different shapes?

KMC SIMULATIONS for small h=0.01

• small p, actived droplet, diagonal interface moves faster than other interfaces.
• large p, vacuum droplet, vertical interface moves faster than other interfaces.

diagonal

vertical

V

p

0

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

Dimension dependence

MF approximation

PRE 2012

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

​

*Mean-field discrete reaction-diffusion equns (lattice diffeq) for the quadratic contact model with the bistability of homogeneous steady states

​

*Orientation-dependence Peq of propagation of planar interfaces

⇒ generic 2PC

​

*Propagation failure for vertical (or horizontal) and diagonal interfaces

​

*Extraordinarily rich droplet dynamics: vacuum droplets exhibit a classical critical size (corresponding to stationary droplet solution); active droplets can exhibit an entire family of stationary solutions

​

*Families of critical droplets exist near generic 2PC.

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

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

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

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]

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

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

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

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

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

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