統計物理在化學主方程式、艾根模型及表面機械研磨處理的應用

(1)

國立臺灣大學理學院物理學系博士論文

Department of Physics College of Science

National Taiwan University Doctoral Thesis

統計物理在化學主方程式、艾根模型及表面機械研磨處理的應用

Application of Statistical Physics in Chemical Master Equation, Eigen Model, and Surface Mechanical Attrition Treatment

黃冠榮

Guan-Rong Huang

指導教授：胡進錕博士 Advisor: Chin-Kun Hu , Ph.D.

中華民國 104 年 9 月

September, 2015

(2)

謝

在台大五年多來，我能夠完成這篇論文，首先必須感謝我所有的指導教授：David Saakian，黃志青，余瑞琳和胡進錕博士。在論文研究問題的討論和修改，總是能給予許多的方向。讓我在每次討論中總是學到新的知識，觀念和方法。

感謝每次解析計算的討論中，David 教授總是一再的講解生物演化和化學主方程式的解析計算方法直到我可以完全理解。以及不遺餘力對於論文每個方程式的檢查與校正。感謝黃志青老師在材料科學計算中給予理論的方向，論文的修改和實驗數據的講解，其中，我也感謝在中山大學共同研究的時期，黃幫實驗室給予的照顧：生活起居和實驗數據的討論。特別是蔡濰猷學弟的幫忙，讓我心無旁鶩專心做研究。

感謝余瑞琳老師在數值方法上的指導，總是耐心地解說每一個步驟每一個原理，讓我在論文的數值計算部分，可以順利完成。感謝胡進錕博士提供我足夠的兼任助理薪水，讓我能夠專心在台大學習研究，以及在每篇論文的修改與審視和論文投稿的協助。

另外，我也感謝在博士班各位同學們互相支援幫助，感謝我的家人在博士班求學過程中全力支持。最後，感謝神讓我有這機會完成這篇論文。

(3)

Abstract

In macro, the method and concept of statistical physics can be a powerful tool in analytical and numerical computation and applied to many fields such as chemical reaction, bio-evolution, and material science. In our work, the methods of statistical physics: chemical master equation (CME), Hamilton- Jacobi equation (HJE), and canonical ensemble are used to calculate various physical quantities. Our work is organized in three topics: the CME with the Gaussian and compound Poisson noise, bio-evolution of Eigen model, and energy conversion in the surface mechanical attrition treatment (SMAT) experiment.

In the CME part, the chemical reaction among DNA, mRNA, and protein can be regarded as a stochastic process. We consider the CME with compound Poisson and Gaussian noises and obtain the exact solution of steady- state probability density function (PDF) verified by the algorithm of forward finite difference in large-scale time. Without Gaussian white noise, the solu- tion of CME (set diffusion coefficient ϵ = 0) can be returned to that of CME derived by Long Cai, et al.

In the bio-evolution part, we use the method of expansion in O(_N¹) to obtain the HJE for probability distribution in Hamming class which is applied to calculate the correction of O(_N¹) accuracy for the steady-state probability distribution in Hamming class and mean fitness in Eigen model. The steady- state distributions of O(_N¹) correction are well-consistent with the Runge- Kutta simulation with relative errors less than 1 %, while the mean fitness

(4)

of O(_N¹) is the same one derived by Michael Deem, et al. in quantum field theory.

In the SMAT part, we consider the collisions among the 304-steel balls, motor top, and chamber bottom, where the chamber or motor can be treated as a hot reservoir. Since we assume that all the collisions among them are elastic except the ball-sample collisions, the balls with negligible potential among them can be regarded as the canonical ensemble. By this concept, we construct the link for energy conversion among the motor top, sample bottom, and balls, where the kinetic energy, heat energy, and internal energy are included in the energy conversion. We also introduce the one-dimensional heat equation with uniform-distributed heat source to obtain the temperature distribution of sample, and we use this temperature distribution of sample to connect the Zenner-Hollmann parameter and the heat energy and surface hardness of sample.

Key words: chemical master equation, Gaussian white noise, compound Poisson noise, bio-evolution, Eigen model, SMAT, collision, energy conversion

(5)

List of Figures

1.1 (a) The desiged-dimension of chamber in the SMAT experiment. (b) The schematic drawing showing that the sample material gains the heat and strain energy from the kinetic energy loss of sample and flying balls. . . . 117 2.1 The Schematic of sample configuration. . . 117 3.1 (a) The average speed of flying balls (in Eq. (3.59)) versus the SMAT

amplitude for the parameters, H = 20 mm, D = 3 mm, ω = 40π krad/s, and _m^m^b

m = 10⁻⁶. (b) The average speed of flying balls (in Eq. (3.59)) versus the SMAT angular frequency for the parameters, H = 20 mm, D = 3 mm, A = 60 µm, and _m^m^b

m = 10⁻⁶. . . 118 3.2 (a) The variation trend of ∆E_k,loss,b predicted by Eq. (3.60) as a function

of ^m_m^b

s for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s,

mb

mm = 10⁻⁶, and e = 0.25, (b) the variation trend of ∆E_k,loss,s predicted by Eq. (3.61) as a function of ^m_m^b

s for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s, _m^m^b

m = 10⁻⁶, and e = 0.25, (c) the vari- ation trend of total energy loss, i.e., the sum of ∆Ek,loss,b and ∆Ek,loss,s

as a function of ^m_m^b

s for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s, _m^m^b

m = 10⁻⁶, and e = 0.25. . . 118 3.3 (a) The averaged time period of flying balls predicted by Eq. (3.64) versus

the SMAT amplitude for the parameters H = 20 mm, ω = 40π krad/s,

mb

mm = 10⁻⁶, and e = 0.25. (b) The averaged time period of flying balls predicted by Eq. (3.64) versus the SMAT angular frequency for the pa- rameters H = 20 mm, A = 60 µm, _m^m^b

m = 10⁻⁶, and e = 0.25. . . 119

(9)

3.4 (a) The variation trend of P_loss,b predicted by Eq. (3.65) as a function of

mb

ms for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s, _m^m^b

m =

10⁻⁶, and e = 0.25. (b) The variation trend of P_loss,s predicted by Eq.

(3.66) as a function of ^m_m^b

s for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s,_m^m^b

m = 10⁻⁶, and e = 0.25. (c) The variation trend of the total power loss, i.e., the sum of P_loss,band P_loss,s as a function of ^m_m^b

s for the parameters H = 20 mm, A = 60 µm, ω = 40π krad/s, _m^m^b

m = 10⁻⁶, and e = 0.25. . . 119 4.1 The mechanism for DNA-mRNA-protein process. . . 120 4.2 The transition PDF for mRNA-protein process. . . 120 4.3 The simulation for the dynamical state of PDF with parameters: a = 0.5

and b = 5 from t = 14 s∼ 4200 s. . . 121 4.4 The simulation for the dynamical state of PDF with parameters: a = 0.5

and b = 5 from t = 5600 s∼ 9800 s. . . 121 4.5 The simulation at t = 25200 s for the dynamical state of PDF and analyt-

ical solution with parameters: a = 0.5 and b = 5. . . 122 4.6 The simulation from t = 8 s ∼ 1000 s for the dynamical state of PDF

with parameters: a = 5 and b = 5. . . 122 4.7 The simulation from t = 1200 s∼ 2000 s for the dynamical state of PDF

and analytical solution with parameters: a = 5 and b = 5. . . 123 4.8 The simulation from t = 2400 s ∼ 12000 s for the dynamical state of

PDF and analytical solution with parameters: a = 5 and b = 5. . . 123 4.9 The simulation from t = 40 s ∼ 1000 s for the dynamical state of PDF

with parameters: a = 8 and b = 8. . . 124 4.10 The simulation from t = 1600 s∼ 9600 s for the dynamical state of PDF

and analytical solution with parameters: a = 8 and b = 8. . . 124 4.11 The steady-state of PDF for Eq. (4.11) with parameters: a = 2 , k = 1,

and γ₂ = 1. . . 125

(10)

4.12 The steady-state of PDF for Eq. (4.11) with parameters: a = 2, ϵ = 0.1, and γ₂ = 1. . . 125 4.13 The simulation from t = 0 ∼ 4480 for the dynamical state of PDF and

analytical solution with parameters: a = 0.5, ϵ = 2×10⁻⁶, γ2 = 2×10⁻³, and k = 1. . . 126 4.14 The simulation from t = 0 ∼ 1200 for the dynamical state of PDF and

analytical solution with parameters: a = 0.5, ϵ = 2×10⁻⁴, γ₂ = 2×10⁻³, and k = 1. . . 126 4.15 The simulation from t = 0 ∼ 4320 for the dynamical state of PDF and

analytical solution with parameters: a = 2, ϵ = 2× 10⁻⁶, γ₂ = 2× 10⁻³, and k = 1. . . 127 4.16 The simulation from t = 0 ∼ 4000 for the dynamical state of PDF and

analytical solution with parameters: a = 2, ϵ = 2× 10⁻⁴, γ₂ = 2× 10⁻³, and k = 1. . . 127 4.17 The simulation from t = 0 ∼ 1800 for the dynamical state of PDF and

analytical solution with parameters: a = 2, ϵ = 2× 10⁻⁴, γ₂ = 2× 10⁻³, and k = 10. . . 128 4.18 The simulation from t = 0 ∼ 1000 for the dynamical state of PDF and

analytical solution with parameters: a = 2, ϵ = 0.02, γ2 = 2× 10⁻³, and k = 1. . . 128 4.19 The probability distributions predicted by Eq. (4.65) and numerical results

with the fitness function and parameters: N = 100, f (m) = e^m², and γ = 1.129 4.20 The probability distributions predicted by Eq. (4.65) and numerical results

with the fitness function and parameters: N = 100, f (m) = e^2m², and γ = 2. . . 129 4.21 The probability distributions predicted by Eq. (4.65) and numerical results

with the fitness function and parameters: N = 100, f (m) = e^2m², and γ = 1. . . 130

(11)

4.22 The temperature distributions for pure Cu predicted by Eq. (4.69) for narrow region near the surface in (a) and for wider region in (b) with the parameters k₀ = 401 W/m· K, q = 0.772 × 10³ ∼ 2.28 × 10³W/mm³, As= 800 mm², L = 1 mm, l = 5 µm, Tb = 398 K, and Tt = 358 K. The temperature distributions for 304 stainless steel predicted by Eq. (4.69) for narrow region near the surface in (c) and for wider region in (d) with the parameters k = 14.9 W/m· K, q = 0.773 × 10³ ∼ 2.28 × 10³W/mm³, A_s = 800 mm², L = 1 mm, l = 5 µm, T_b = 378 K, and T_t = 318 K.

The different colored lines correspond to various percentages of kinetic energy loss which is converted into the heat energy of sample. . . 131 4.23 The temperature distributions for pure Cu predicted by Eq. (4.71) for nar-

row region near the surface in (a) and for wider region in (b) with the parameters k₀ = 401 W/m· K, q = 0.772 × 10³ ∼ 2.28 × 10³W/mm³, A_s= 800 mm², L = 1 mm, l = 5 µm, T_b = 398 K, and T_t = 358 K. The temperature distributions for 304 stainless steel predicted by Eq. (4.71) for narrow region near the surface in (c) and for wider region in (d) with the parameters k = 14.9 W/m· K, q = 0.773 × 10³ ∼ 2.28 × 10³W/mm³, A_s = 800 mm², L = 1 mm, l = 5 µm, T_b = 378 K, and T_t = 318 K.

The different colored lines corresponds to various percentages of kinetic energy loss which is converted into the heat energy of sample. . . 132 5.1 (a) The cross-sectional SEM micrograph taken from the sample subject to

SMAT with the 2 mm flying balls and 40 µm SMAT amplitude. (b) The gradient variation trend of hardness of selected SMAT 304 SS samples. . 132 5.2 The relationship between the resulting grain size and Zener-Holloman pa-

rameter with the sample processed by different SMAT conditions. . . 133

(12)

List of Tables

1 The coefficient table for the forward finite difference of f^′(x). . . 110 2 The coefficient table for the forward finite difference of f^′′(x). . . 110 4.1 The comparison of our results among P (m), P₁(m), and numerics for the

fitness function and parameters: f (m) = e^m², γ = 1, N = 100. . . 130 4.2 The comparison of our results among P (m), P₁(m), and numerics for the

fitness function and parameters: f (m) = e^2m², γ = 2, N = 100. . . 130 4.3 The comparison of our results among P (m), P₁(m), and numerics for the

fitness function and parameters: f (m) = e^2m², γ = 1, N = 100. . . 130

(13)

Chapter 1 Introduction

This chapter is organized as the three following topics: the chemical master equation, bio-evolution, and SMAT experiment in a brief and general introduction. It makes the connection and correspondence among the techniques and concepts of statistical physics.

Each section states the topic background and physical meaning of each equation.

1.1 Chemical Master Equation

The coupled chemical differential equations which are equivalent to a number of chemical reactions are common model to describe the chemical reactions among molecules, where such equations are described by variables: time-dependent concentration of each molecule and constants of temperature-dependent reaction rate. The changes of concen- trations with time can be modelled by differential equations with the large number of each interacting molecules. In the same circumstance, two or more reactions can take place simultaneously. The meaning for the collection of coupled ordinary differential equations is that these reactions occur concurrently in the solution.

For the small number of interacting molecules, however, the simple deterministic process breaks down. Molecules are collided by stochastic process (drift or diffusion) described well by Ito lemma, so chemical reactions can’t well described by some simulta- neous processes. With the introduction of probability density function (PDF), P (x, t), in terms of concentration for each molecule at a given time, the time-evolution of PDF

(14)

represents that reactions take place randomly among any possible reactions.

In recent years, biophysicists have paid more attention to stochastic dynamics in cell biology [1–6]. Friedman, Cai and Xie (FCX) obtained an partial differential equation (PDE) to describe the steady-state PDF of protein concentration for living cells in gene expression problem [1]. Losick and Desplan found that noises can induce cells to switch between different gene states in their experiments [2]. Thus, the PDF of stochastic process is a tool to describe stochastic reactions in cell biology.

From Itõ’s lemma [7], the total differential of concentration

dx = b(x)dt +√

2ϵdB, (1.1)

, where B is the Brownian motion, we can derive the corresponding Kolmogorov forward equation (KFE), i.e. the chemical master equation (CME),

∂P (x, t)

∂t = ϵ∂²P (x, t)

∂x² − ∂

∂x[b(x)P (x, t)], (1.2)

obeyed by PDF related to large deviation function or WKB expansion in quantum me- chanics, P (x, t) = e⁻¹^ϵ ^u^ϵ^(x,t). The equation for large deviation function is:

∂u_ϵ(x, t)

∂t = −[(∂u_ϵ(x, t)

∂x )²+ b(x)∂u_ϵ(x, t)

∂x ] + ϵ[∂²u_ϵ(x, t)

∂x² +db(x)

dx ], (1.3)

which reduces to Hamilton-Jacobi equation (HJE) with ϵ = 0, where x is the time- dependent concentration, and ϵ is a small perturbed constant. In classical mechanics (CM), the right-handed term is negative Hamiltonian, and u_ϵ(x, t) is the generating function for corresponding canonical transformation (CT).

For correspondence, we can consider a chemical reaction whose concentration x de- scribed by ordinary differential equation (ODE),

dx

dt = b(x), (1.4)

(15)

where the ODE can be re-written into Itõ’s lemma by the introduction of B white noise after diffusive perturbation. Equation (1.4) is equivalent to Eq. (1.1) with ϵ = 0, and the u(x, t) = lim_ϵ_→0u_ϵ(x, t) is called the principal function furnishing the entire family of or- bits corresponding to Hamiltonian system in phase space. The corresponding Hamiltonian has the form,

H(q, p) = p²+ b(q)p, (1.5)

which follows the equation of motion,

˙ q = ∂H

∂p = 2p + b(q)

˙

p =−∂H

∂q =−pdb(q)

dq , (1.6)

and has the following Lagrangian,

L(q, ˙q) = [p ˙q− H(q, p)]p=¹₂( ˙q−b(q)) = 1

4( ˙q− b(q))², (1.7)

which corresponds to the action functional,

S₀[q(t); (0, q(0)) → (t, q(t))] =^∫ ^t

0

dτ1 4[dq(τ )

dτ − b(q(τ))]². (1.8)

This is equivalent to path integral in quantum mechanics by making the integration path along imaginary axis in complex plane, and the probability of the system is proportional to e^−S⁰^/ϵ⁰ which is exactly the probability for a path in stochastic dynamics in Eq. (1.4) with ϵ = 0. For finite ϵ, the action functional generalized by Onsager and Machlup is S_ϵ = S₀+₂^ϵb^′(q) [8–10].

On the other hand, the corresponding CME can be obtained by the chemical system with n molecules,

∂P (n, t)

∂t = N [R+(n− 1)P (n − 1, t)

+ R₋(n + 1)P (n + 1, t)− [R+(n) + R₋(n)]P (n, t)], (1.9)

(16)

where N is a large integer and a maximal allowed number of molecules, R₊is the growth rate, and R₋ is the degradation rate. And the equation should be modified at the border.

Let us use the variable, x = _Nⁿ, to re-write Eq. (1.9) in terms of x as

1 N

∂P (x, t)

∂t = R₊(x− 1

N)P (x− 1 N, t) + R₋(x + 1

N)P (x + 1

N, t)− [R+(x) + R₋(x)]P (x, t). (1.10)

With the largeness of N , x becomes continuous from discrete. By using the ansatz (Or WKB method),

P (x, t) = exp[N u(x, t)], (1.11) to construct a equation for P (x, t). By using HJE with N → ∞ [11–15], the partial differential equations can be solved exactly as the following steps:

∂u

∂t + H(x, u^′) = 0,

H(x, u^′) = [R+(x)− R−(x)]u^′,

and the corresponding equation of motion for x:

˙x(t) = ∂H

∂u^′

= R₊(x)− R₋(x) = b(x), (1.12)

where u^′ = ^∂u_∂x. The distribution variance has been derived in [15]:

b²(x)

∫ _x

x0

c(y)

b³(y), (1.13)

where c(y) = R+(y) + R₋(y). With these results, we can formulate new CME with different process or noise.

The rest of CME topic is organized as follows: Firstly, we use HJE as a tool to investigate the physical of chemical system and solve the CME. Secondly, we introduce the drift-diffusion process to generalize PDF for chemical reactions [16, 17], which de-

(17)

scribes how to make the correspondence probability of chemical reaction by path integral and formulate CME with white noise [18]. Finally, we introduce the hybrid model of Gaussian and Poisson noise solved exactly and verified by the forward finite difference simulation [19], and we make conclusion about this hybrid model.

1.2 Bio-evolution

It’s widely acknowledged that DNA can carry hereditary code to determine the life performance, and it has the important influence on the reproduction and survival for each specie. As we all know that the positive self-regulation or mutation of gene can help species adapt to their surroundings for survival, and such processes called bio-evolution.

Exact bio-evolution results, the mean fitness, steady state distribution , and dynamics for genes are obtained by the tool of statistical physics and mathematics [20], HJE method, partial differential equation, and numerical simulation, which makes possibility to realize the mechanism of virus or cancer evolution.

In past decades, the bio-evolution process of virus is described very well by the Eigen and Crow-Kimura models for large population size or genome size [20–27]. In the two models, the fractions of population for different types p_i are described by a set of deterministic partial differential equations. The corresponding HJE for the two models can be obtained by expanding p_i in the order of O(_N¹), p_i = exp (N u_i), where N is the genome length. Since N is possibly 40∼ 100 [28], it’s important to investigate the finite popula- tion effect of _N¹ in a small N genome.

For the sake of simplicity, we can assume there are two different genotypes (denoted as ±1) for each nucleotide, and the N-nucleotide genome has 2^N types. We have the following system of equations in Eigen model for each probability p_i,

∂p_i

∂t =

2^N

∑

j=1

(Qij)rjpj − pi(

2^N

∑

j=1

rjpj), (1.14)

where the p_i satisfies^∑_ip_i = 1, the mean fitness r_iis the mean number of offspring per unit cycle for type i sequence, and mutation matrix with the mean nucleotide incorporation

(18)

fidelity q is expressed as:

Q_ij = q^N^−d^ij(1− q)^d^ij, (1.15) where d_ij is the Hamming distance (HD) between two sequences, S_i and S_j, defined as:

d_ij = (N −^∑^N

l=1

s^l_is^l_j)/2, (1.16)

where s^l_iis the spin with possible values±1 at l-th site in Si. To simplify the HD between sequences, we can choose a reference sequence S₀ with all spins being +1. Without loss of generality, sequences with the same number of−1 spin are assumed to have the same probability, namely, it’s symmetrical distribution. Thus, 2^N types are divided into N + 1 Hamming class, and the HD can be written as

d_l0 = (N −^∑^N

i=1

sⁱ_l)/2 = l, (1.17)

where 0 ≤ l ≤ N and l is the number of −1 spin for a sequence. With symmetrical distribution, we can also assume the mean fitness r_iin terms of m is symmetrical,

r_i = f (m), (1.18)

where m = 1 − ^2l_N between ±1 called magnetization and f(m) is called the fitness function. Therefore, the original Eigen model with 2^N equations is transformed into Eigen model with N + 1 equations with which is easy to be treated. This model with p_i = exp[N u(m, t) + u₁] expansion is lead to HJE equation as we did in CME and other high order of_N¹ equations on which the work of finite correction is based.

The rest of bio-evolution topics is organized as follows: Firstly, we introduce Crow- Kimura model and derive its some properties. Secondly, HJE is obtained by WKB expan- sion of p_i to investigate the characteristic of Crow-Kimura model, and then we develop HJE method and derive some useful formula from Crow-Kimura model which can be a useful tool for the finite correction of Eigen model. Finally, we solve Eigen model with N + 1 class in O(¹) accuracy, and it’s verified by the numerical simulation, Runge-Kutta

(19)

method.

1.3 SMAT modelling

In traditional engineering treatments, shot peening by using steel balls to bombard onto metal surfaces has been adopted to leave compressive residual strains within the affected region in promoting the fatigue properties [29, 30]. The balls have typical diameters of 0.1 ∼ 2 mm and gain their speed by compressed air. Normally, these balls bombard the metal surfaces in the frequency range of 20∼ 100 Hz and speed range of 50 ∼ 100 m/s.

A new physical treatment named ultrasonic surface mechanical attrition treatment (SMAT) was firstly introduced in 1999 [31, 32]. The SMAT balls are accelerated and bombarded by the ultrasonic motor on the chamber bottom, as shown in Fig. 1.1 (a).

The diameter and speed of flying balls and the bombarding frequency are in the range, 1 ∼ 10 mm, 1 ∼ 20 m/s, and 10 ∼ 100 kHz [33], respectively. The most important feature is that the incident direction onto the metal surface can be designed to vary with time lapse in making smaller grain size of metallic materials. This will lead to promising properties of metals such as grain refinement and gradient structure. And researchers have made extensive uses of this SMAT treatment on various metallic materials including pure iron [34], stainless steel [35], and pure copper [36] in making gradient and nano-crystalline structure [37–40].

Although SMAT has gradually developed into a matured engineering surface treatment, never before has SMAT been treated with rigorous analytical modelling. Therefore, a systematic SMAT model is actually needed. Here, we consider the interaction between flying balls and chamber imagined as a canonical ensemble, where chamber is reservoir giving balls the kinetic, internal, and heat energy. The chamber volume is much greater than balls volume, so collision frequency between balls is small compared to that between balls and chamber and balls interaction can be neglected. The motion of motor top is characterized by longitudinal harmonic motion,

v = 2Aπνsin(2πνt), (1.19)

(20)

where v_m is the velocity of motor top, A is the amplitude, and ν is the angular frequency.

To construct the relation of energy conversion between motor top and balls, the ball-motor collision can be counted as elastic collision. Thus, the induced velocity of ball, v_b, is described by

v_b = 2m_mv_m

m_b + m_m ≈ 2vm, m_m >> m_b, (1.20) where m_b and m_m are the mass of each ball and motor. On the other hand, the ball- sample collision is assumed to be inelastic collision with restitution constant e obeying conservation of momentum,

e = v_s^′ − vb^′

v_b− vs

, m_sv_s+ m_bv_b = m_sv^′_s+ m_bv_b^′, (1.21)

where m_sand v_sare sample mass and velocity, respectively.

The kinetic energy of flying balls is not conserved due to the inelastic ball-sample collision. The kinetic energy loss for flying balls and sample in the SMAT chamber can be mainly converted into three parts as indicated in Fig. 1.1 (b). Firstly, it is the strain energy of sample due to the formation of dislocations and vacancies [41–43]. Secondly, it is the heat energy of sample, where the heat flow and its temperature distribution are both important factors for the resulting metal micro-structure [44]. Recrystallization might be taken place in the sample while the experienced temperature reaches some critical values.

Finally, it would be the sonic energy and heat energy in the chamber originated from the inelastic collisions between flying balls.

In SMAT topics, we have established connections among the parameters of flying balls, the ball size, flying speed, the bombarding frequency and amplitude of motor motion, the height of chamber, and the energy and power of sample. During the SMAT processing time, we can find the input energy and power of sample through these connections. The condition for the frequency of flying ball reaching a stead speed can also be obtained in this approach. For the heat energy of sample, we have introduced the one- dimensional heat equation with the uniformly-distributed heat source to estimate the heat flow and temperature distribution of sample which are hard to be measured in the SMAT

(21)

experiment. With the temperature distribution of sample, we can make connection among the strain rate, hardness, and grain size of sample. With these connections and modelling, one can find an optimized approach to the mechanical performance of metal surface via the SMAT experiment.

(22)

Chapter 2 Basic theory for CME, Bio-evolution, and SMAT

2.1 Hamilton-Jacobi Equation

In Classical Mechanics (CM) [45], theoretical physicists use independent variable of position x and momentum p for each particle in constructing Hamiltonian to characterize the particle dynamics, where the correspond equation of motion constructed by Hamilto- nian for x and p is the Hamilton equation. To simplify the equation of motion, we prefer the physical system in which all of the generalized coordinates or all of the canonical mo- mentums are cyclic, namely, Hamiltonian with respect to x and p is a constant. To achieve the goal, a canonical transformation (CT) under which the equation of motion is invariant should be found, and any CT corresponds to a generating function consisted of half new and half old generalized coordinates. Here comes the Hamilton-Jacobi equation (HJE) which generating function satisfies.

In stochastic dynamics, a fictitious Hamiltonian which is similar to that of CM as a tool to formulate the Lagrangian and the action functional by using the WKB expansion, P (x_t, t) = e⁻¹^ϵ^u^ϵ^(x^t^,t)with small ϵ, where ϵ results from diffusion. And such WKB expan- sion we use in bio-evolution is P (x, t) = e^{N u(x,t)}with large N , where the N is genome length or population size. It’s obviously that both expansions are mathematically equiva-

(23)

lent for small ϵ or large N . In our work for CME or bio-evolution, we put such expansion into the equation of motion (KFE or evolution model) with condition ϵ→ 0⁺or N >> 1 to get the HJE for u. As ϵ→ 0⁺or N >> 1, lim_ϵ_→0⁺u_ϵ(x_t, t) is called the large deviation rate function in probability theory [46] or the principal function in CM [45]. The large deviation rate function asymptotes the behaviour of P (xt, t) as ϵ → 0⁺ or N → ∞, and the principal function furnishes the entire family of orbits corresponding to a Hamiltonian system in phase space.

In our research, HJE in chemical reaction derives the path probability of each reaction while HJE in bio-evolution derives a series of equation for finite correction of O(_N¹). Here we give some derivation and introduction to HJE in CM.

2.1.1 Canonical Transformation

In CM, the form of Hamilton’s equations are invariant under canonical transformation (CT) and Hamilton’s principle states that the most possible track of classical system makes the corresponding action functional minimized. Thus, we have the variation of action functional is 0 over the most possible track:

δ

∫ _t₂

t1

L(q, ˙q, t)dt = 0, δ

∫ _t₂

t1

L^′(q, ˙q, t)dt = 0, δ

∫ _t₂

t1

[p_iq˙_i − H(q, p, t)]dt = 0, δ

∫ _t₂

t1

[P_iQ˙_i− K(q, p, t)]dt = 0,

which implies:

λ[p_iq˙_i− H(q, p, t)] = PiQ˙_i− K(q, p, t) +dF dt ,

where q and p are old generalized coordinates, Q and P are new generalized coordinates, L, L^′, H, and K are the corresponding Lagrangian respectively, λ is a scale constant, and F is corresponding generating function in terms of half new and half old coordinates. For

(24)

the λ = 1 case, it’s the case called CT in CM. With λ = 1, the equation becomes:

piq˙i− H(q, p, t) = PiQ˙i− K(q, p, t) +dF

dt , (2.1)

where F called the generating function has four basic forms by [45],











F = F₁(q, Q, t)

F = F₂(q, P, t)− PiQ_i F = F₃(p, Q, t) + p_iq_i

F = F₄(p, P, t) + p_iq_i− PiQ_i

. (2.2)

2.1.2 Hamilton-Jacobi Equation

We start the derivation of Hamilton-Jacobi equation (HJE) from Eq. (2.1) and put the generating function F₂ in Eq. (2.2) to arrive:

p_iq˙_i− H(q, p, t) = PiQ˙_i− K(Q, P, t) + dF dt

= P_iQ˙_i− K(Q, P, t) + ∂F₂

∂t +∂F₂

∂qi

˙

q_i+∂F₂

∂pi

˙

p_i− ˙PiQ_i− PiQ˙_i

= −K(Q, P, t) − ˙PiQi+∂F₂

∂t +∂F₂

∂q_i q˙i+∂F₂

∂p_i p˙i.

Then we rearrange it and have the equation:

[H(q, p, t) +∂F2

∂t − K(Q, P, t)] + (∂F2

∂q_i − pi) ˙qi+ (∂F2

∂P_i − Qi) ˙Pi = 0, (2.3)

where ˙q_i and ˙P_i are separated independently. Since the three terms in Eq. (2.3) are independent of each other, the three terms must be 0 to hold the equality. Therefore,

∂F₂

∂q_i = pi,∂F₂

∂P_i = Qi, H(q, p, t) +∂F₂

∂t = K(Q, P, t). (2.4)

(25)

To make all generalized coordinates cyclic, we can set K(Q, P, t) = 0. And the corre- sponding Hamilton’s equations are:

Q˙_i = ∂K

∂P_i = 0, P˙_i =−∂K

∂Q_i = 0,

which means that all generalized coordinates are constants of motion. On the other hand, making K(Q, P, t) = 0 gives the corresponding equation for F₂:

H(q, p, t) +∂F₂

∂t = 0,

→ H(⃗q,∂F₂

∂⃗q , t) + ∂F₂

∂t = 0, (2.5)

where Eq. (2.5) is called Hamilton-Jacobi equation and F₂ is the Hamilton’s principal function in CM which is the counterpart of u(x, t) in stochastic dynamics. To investigate the physical meaning of F₂, we can take its total differential with respect to t:

dF₂

dt = ∂F₂

∂q_i q˙i+∂F₂

∂P_i

P˙i+∂F₂

∂t

= p_iq˙_i− H(q, p, t) = L(q, ˙q, t),

and then we integrate ^dF_dt² back with respect to t:

F₂ =

∫ _t

t0

L(q, ˙q, t^′)dt^′,

which states that the Hamilton’s principal function F2 is equivalent to action functional.

Thus, in physics, solving the HJE is equivalent to solve the variation equation of action functional, Euler-Lagrange equation.

2.1.3 HJE Application in CME

The HJE method in CM has been well developed for hundreds years, and it is a power and analytical tool to investigate the characteristic of physical system in macro. Thus, we

(26)

want to develop the HJE method in CME to help us realize the mechanism of chemical sys- tem. As stated in previous sections, the principal function u(x, t) in stochastic dynamics is similar to the role of generating function in CM. Here we introduce a fictitious Hamil- tonian corresponding to HJE in CME. In the introduction of CME, we have the following CME for general chemical system with n molecules:

∂P (x, t)

N ∂t = R₊(x− 1

N)P (x− 1

N, t) + R₋(x + 1

N)P (x + 1 N, t)

− [R+(x) + R₋(x)]P (x, t), (2.6)

where N is the maximally-allowed number of molecule, x = _Nⁿ, and R₊and R₋ are the rate of generation and degradation. By Taylor expansion with largeness of N at x and P (x, t) = exp [N u(x, t)], we have:

R₊(x− 1

N)≈ R+(x)− 1 N

∂R₊(x)

∂x , R₋(x + 1

N)≈ R₋(x) + 1 N

∂R₋(x)

∂x , P (x∓ 1

N, t)≈ P (x, t) ∓ 1 N

∂P (x, t)

∂x

= P (x, t)[1∓ ∂u(x, t)

∂x ],

∂P (x, t)

∂t = N∂u(x, t)

∂t P (x, t).

Put these expansions above into Eq. (2.6) to get the equation of zero order in _N¹:

∂u

∂tP (x, t) = [R₊(x)(1− ∂u

∂x) + R₋(x)(1 +∂u

∂x)− (R+(x) + R₋(x))]P (x, t), (2.7)

and divide both sides of Eq. (2.7) by P (x, t) and let u^′ = ^∂u_∂x to obtain:

∂u

∂t + [R₊(x)− R₋(x)]u^′ = ∂u

∂t + H(x, u^′) = 0. (2.8)

Equation (2.8) has the exact form of HJE, where u corresponds to generating function in CT and H corresponds to Hamilton in CM. Thus, the fictitious Hamiltonian and equation

(27)

of motion are shown as:

H(x, u^′) = [R+(x)− R−(x)]u^′, x(t) =˙ ∂H(x, u^′)

∂u^′ = R₊(x)− R₋(x) = b(x),

where b(x) = R₊(x)− R−(x). It is reasonable for chemical reaction of zero order whose concentration obeys the ordinary differential equation. And the variance of P (x, t) has been derived in [15]:

b²(x)

∫ _x

x0

c(y)

b³(y)dy, (2.9)

where c(y) = R₊(y) + R₋(y) and x₀ is the reference point. Therefore, each chemical system actually corresponds to a fictitious Hamiltonian system derived from CME of each chemical system.

2.2 Diffusion Process

Stochastic or random process is always related to diffusion process which can be traced back to Brownian motion. The random motion of particles suspended in water and re- ported by Robert Brown in 1827 is known as Brownian motion or Wiener process, which is the most important case in stochastic process. Some scientists in earlier periods con- sidered that Brownian motion is caused by living cells, and Poincaré thought this motion violates the second law of thermodynamics. Now scientists consider that such molecule motion is induced by the continuous collision of molecules around them. For general sit- uation, a specific molecule undergoes 10²⁰collisions per second. In 1905, Albert Einstein described Brownian motion in term of PDF equation by the kinetic molecular theory. He proved that the PDF motion satisfied the following partial differential equation (PDE),

∂P

∂t = D∂²P

∂x², (2.10)

(28)

where the positive D is a diffusion coefficient and x is the particle position. By changing of variables, y = ^√^x

2D, the PDE equation can be transformed into the form of heat equation:

∂P

∂t = 1 2

∂²P

∂x², (2.11)

whose well-known solution is

P (x, t) = 1

√4πDtexp[−(x− x0)²

4Dt ], (2.12)

where x₀ is the mean position of particle, t is the time, and it is Gaussian distribution for any given time.

Kiyoshi Itõ, a Japanese mathematician, came up a good idea to describe the drift- diffusion process by Itõ’s lemma:

dX_t= µ(X_t, t)dt + σ(X_t, t)dB_t, (2.13)

which states that the random variable dB_tmakes an impact on other deterministic variables in a small time interval ∆t; the expectation value E[dXt] is unchanged for an entire cycle.

The derivatives of Itõ’s lemma is the stochastic differential equation which is widely used in financial and biological physics. From Itõ’s lemma, we also derive the equation of motion for PDF, Kolmogorov forward equation (KFE), which is substantial to many fields related to stochastic process.

2.2.1 Brownian Motion

Brownian motion is the process of foundation for various stochastic processes, and here we derive it from central limit theorem. Considering one particle, it undergoes a collision to step ∆x displacement after a time-interval ∆t which is independent of the particle position. Also we consider the probabilities of stepping ∆x and−∆x are p and 1−

p respectively. With largeness of container volume, particles are far away from container border. The particle motion can be treated as the independent 1D-random walk. The

(29)

particle location at a specific time t, X(t), is expressed as:

X(t) = ∆x(I1+ I2+ ... + I_[ ^t

∆t]), (2.14)

where the I_i is 1 or−1 depending on i-th displacement which is +∆x or −∆x, [ ] is the Gauss function, and the corresponding probabilities for different displacement (±∆x) are:

P (I_i = 1) = p, P (I_i) = 1− p. (2.15) To simplify the derivation of Brownian motion, we can set the following variables:

∆x = σ√

∆t, p = 1 2+

√∆t 2σ µ, n = lim

∆t→0[ t

∆t] = lim

∆t→0

t

∆t. (2.16)

We can prove Eq. (2.16) by sandwich theorem:

t

∆t ≥ [ t

∆t]≥ t

∆t + 1,

∵ lim

∆t→0

t

∆t + 1 = lim

∆t→0

t + ∆t

∆t = lim

∆t→0

t

∆t,

∴ lim

∆t→0[ t

∆t] = lim

∆t→0

t

∆t. (2.17)

We have the expectation value of I_i∆x,

E[I_i∆x] = p× ∆x + (1 − p) × (−∆x)

= (1 2 +

√∆t 2σ µ−1

2 +

√∆t 2σ µ)σ√

∆t

= µ∆t, (2.18)

and the variance of I_i∆x,

V [I_i∆x] = E[(I_i∆x)²]− E²[I_i∆x]

= E[∆x²]− µ²(∆t)² = (σ²− µ²∆t)∆t,

(30)

where µ and σ are both time-independent constants. For p = ¹₂ case, the E[I_i∆x] = 0 and V [I_i∆x] = σ²∆t. Therefore, the expectation value and variance of X(t) for p = ¹₂ case are:

E[X(t)] = nE[I_i∆x] = 0, V [X(t)] = nσ²∆t = σ²t,

and the position distribution of particle can be derived from central limit theorem as ∆t→ 0:

z = X(t)− nµ∆t σ√

n∆t = X(t)− X0

σ√

t → 1

√2πe^−z2² ,

→ 1

√2πexp [−(X(t)− X0)²

2σ²t ]→ N(X0, σt), (2.19)

where N(X₀, σt) is the Gaussian distribution with expectation value X₀ = µt and variance σ²t for a given time t. The motion with µ = 0 and σ² = 1 is called standard Brownian motion (SBM). Since any Brownian motion can be transformed into SBM, only SBM have to be taken into account. Thus, such SBM lead to the important conclusion that each particle position is normally distributed with expectation value 0 and variance 1.

2.2.2 Itõ’s Lemma

In previous sections, the deterministic term µ which is so-called drift term is interfered by non-deterministic term σ which is so-called diffusion term, and any drift-diffusion process is well described by Itõ’s lemma. In physics, we also have the corresponding counterpart that the specific track of each quantum particle is unpredictable due to the wave-particle property, so do we have deterministic term (mean particle position) and non-deterministic term (track). For the total differential of any function f (x, t) in terms of x and t with the influence of B_t, the useful consequence,

df (x, t) = (µ∂f

∂x + σ² 2

∂²f

∂x² +∂f

∂t)dt + σ∂f

∂xdB_t, (2.20)

(31)

can be derived from Itõ’s lemma. Before proving the formula above, we need to prove the equality:

dB_t² = dt. (2.21)

To prove this, we at first set S with t_n= t as:

S = lim

n→∞

∑n k=1

(B_t_k − Btk−1)², (2.22)

where B_t_k − Btk−1 is Brownian motion in the time interval, ∆t = t_k− tk−1, and B_t_k − Bt_k−1 = dBt_k and ∆t = tk − tk−1 = dt in the case of n → ∞. Then we take the expectation of S:

E(S) = lim

n→∞

∑n k=1

(B_t_k− Btk−1)²

= lim

n→∞

∑n k=1

E[(B_t_k− Btk−1)²]]

= lim

n→∞

∑n k=1

(t_k− tk−1) = t, (2.23)

and the variance of S:

V (S) = lim

n→∞

∑n k=1

V [(B_t_k− Bt_k−1)²]

= lim

n→∞

∑n k=1

{E[(Bt_k− Bt_k−1)⁴]− E²[(Bt_k − Bt_k−1)²}

= lim

n→∞

∑n k=1

[3(t_k− tk−1)² − (tk− tk−1)²]

= lim

n→∞

∑n k=1

2(t_k− tk−1)²,

(32)

where here we have used the result,^∫_−∞^∞ x⁴N(0, ∆t)dx = 3∆t² = 3(t_k− tk−1)². We can derive V (S) = 0 as follows:

V (S) = lim

n→∞

∑n k=1

2(t_k− tk−1)²

≤ lim_n_→∞2 max

k (t_k− tk−1)

∑n k=1

(t_k− tk−1)

= 2t lim

n→∞max

k (tk− tk−1) = 0. (2.24)

As n→ ∞, the corresponding integral to E(S) is:

E(S) =

∫ _t

0

dB_t²′ = t, (2.25)

which states dB_t² = dt, and V (S) = 0 implies that dB_t² is a measurable variable. Thus, we can immediately deduce that dB_t² = dt. On the other hand, we have the corresponding counterpart in quantum mechanics, Heisenberg Uncertainty Principal, which states: If two physical observables are measured simultaneously, their commutator is 0. Here comes the simple proof from this principle,

[dB²_t, dt] = 0→ dB²t = cdt, dt = E[dB_t²] = cE[dt] = cdt,

→ c = 1, ∴ dBt² = dt, # (2.26)

where c is a phase constant. Now we turn to derive Eq. (2.20) by this equality. We have the following Taylor expansion near (x, t) point,

f (x + ∆x, t + ∆t) ≈ f(x, t) +∂f

∂x∆x + ∂f

∂t∆t

+ 1

2!

∂²f

∂x²(∆x)² +∂f

∂x

∂f

∂t∆t∆x + 1 2!

∂²f

∂t²(∆t)².

From Eq. (2.13), we also have:

∆x = µ∆t + σ∆B .

統計物理在化學主方程式、艾根模型及表面機械研磨處理的應用

國立臺灣大學理學院物理學系 博士論文

Department of Physics College of Science

National Taiwan University Doctoral Thesis

統計物理在化學主方程式、艾根模型及表面機械研磨處理的 應用

Application of Statistical Physics in Chemical Master Equation, Eigen Model, and Surface Mechanical Attrition Treatment

黃冠榮

Guan-Rong Huang

指導教授：胡進錕博士 Advisor: Chin-Kun Hu , Ph.D.

中華民國 104 年 9 月

September, 2015

謝

Abstract

Contents

List of Figures

List of Tables

Chapter 1 Introduction

1.1 Chemical Master Equation

1.2 Bio-evolution

1.3 SMAT modelling

Chapter 2

Basic theory for CME, Bio-evolution, and SMAT

2.1 Hamilton-Jacobi Equation

2.1.1 Canonical Transformation

2.1.2 Hamilton-Jacobi Equation

2.1.3 HJE Application in CME

2.2 Diffusion Process

2.2.1 Brownian Motion

2.2.2 Itõ’s Lemma

國立臺灣大學理學院物理學系博士論文

統計物理在化學主方程式、艾根模型及表面機械研磨處理的應用