二維磁偶極系統的隨機共振

(1)

國立交通大學

物理研究所

碩士論文

二維磁偶極系統的隨機共振

Stochastic resonance in two-dimensional Magnetic dipoles

system

研究生：李宜芳

指導教授：張正宏教授

(2)

二維磁偶極系統的隨機共振

Stochastic resonance in two-dimensional Magnetic dipoles

system

研究生：李宜芳 Student: Yi-Fang Li

指導教授：張正宏 Advisor: Cheng-Hung Chang

國立交通大學

物理研究所

碩士論文

A Thesis

Submitted to Institute of Physics National Chiao Tung University In partial Fulfillment of the Requirements

for the Degree of Master

in

Physics

August 2011

Hsinchu, Taiwan, Republic of China

中華民國一百年八月

(3)

摘要

隨機共振為噪音強化弱週期訊號的非平衡統計力學問題。過去對於隨機共振現象的理論大多是從數學模型出發，例如考慮雙穩態位能井，討論訊號與信噪比的結果。此論文從一個實際例子出發，發現兩個耦合的磁偶極系統也有隨機共振現象，利用現有的理論可以大致解釋觀察到的結果，而對模擬與理論的偏差我們提出解釋，例如質量或轉動慣量與第二路徑的影響。最後我們研究多個磁偶極的系統並與兩個磁偶極系統得到之結果做比較。

(4)

誌謝

物理所的氛圍一直都是溫馨和樂且適合做學術研究的，這兩年中謝謝所上的老師們給予我的幫助與指導，讓我的研究生活過得充實、踏實，而在行政業務上要特別的感謝親切的楦敏和俐君小姐，因為有妳們的提醒讓我在各方面都能夠順利的進行。也感謝我的同學們，和你們一起討論作業的時光令人難忘且值得回味再三，因為有大家互相惕勵著，讓我們都能夠完成口詴、而後畢業。在這裡我要特別的感謝我的指導教授張正宏老師，老師花了非常多時間耐心的與我們討論及分享研究及生活中的事物，讓我成長很多且深深的感受到研究之路的寬廣與壯美；另外還要感謝德明學長與宏慶學長給予我的協助，讓我能夠在有限的時間裡習得充足的資訊並且進一步的做應用。還要感謝我的家人們，因為有你們的支持讓我今天能夠站在這裡；總覺得感謝的話怎麼樣也說不完，容我在這裡對所有的人說一聲，謝謝大家。

(5)

LIST OF FIGURES

Fig. 1-1 雙向環型雷射輸出的光訊號強度對噪音強度的關係圖(信噪比) (圖片擷

取自 McNamara et al. Theory of stochastic resonance,physical

reviewA,volume39) ... 1

Fig.1-2 (a)無外場情形時的雙穩態位能圖；當系統中加入一弱外場A₀sin(t)，(b) 為當t不同的情況下位能之改變。其中位能井中的粒子若是在某個位能井中有較低的機率時，則粒子用虛線呈現。 ... 3 Fig. 2-1 逃逸問題示意圖。 ... 8 Fig. 2-2 14 Fig. 3-2 兩個磁偶極系統的交互作用能量關係圖。 ... 21 Fig. 3-1 兩磁偶極以頭尾相接的方式沿著 z 軸排列。... 21 Fig. 3-3 系統的能量面在隨時間變的外場 B=3sin(wt)下之變化。 ... 22 Fig. 3-4 系統兩個狀態的劃分示意圖。 ... 24 Fig. 3-5 達平衡時動能隨時間的變化。 ... 27 Fig. 3-6 D=1 時 x 之圖形。 ... 29 Fig. 3-7 粒子經過最可能路徑到達另一個位能井之躍遷速率隨著時間與噪音強度 D 之圖形。 ... 29 Fig. 3-8 粒子經過最可能路徑達到另一個位能井之平均首通時間。 ... 30 Fig. 3-9 D=1 時 x 之圖形。 ... 31 Fig. 3-10 粒子經過另一路徑之平均首通時間。 ... 31

(8)

Fig. 3-12 Welch 對資訊點之切割示意圖。 ... 33 Fig. 3-13 由 Matlab 所產生的高斯白噪音與其功率譜。左圖為一段時間序列之下噪音的振幅，右圖是利用 Welch 方法所得到的功率譜密度，注意 x 軸是 Log 的尺度。 ... 34 Fig. 3-14 系統處於第一個狀態的個數隨著噪音 D 改變之變化。 ... 35 Fig. 3-15 隨著噪音強度 D 的增加，處於第一個狀態的系統數目之頻譜密度；圖中的峰值處對應到外場的頻率。 ... 36 Fig. 3-16 在不同噪音強度之下，對應到外場頻率處的輸出訊號功率譜密度之強度。 ... 37 Fig. 3-17 系統由一個穩定態經過最可能路徑到達另一個穩定態之躍遷個數。 .. 38 Fig. 3-18 系統輸出之信噪比與噪音強度之關係圖。 ... 39 Fig. 3-19 理論計算之訊號、噪音及信噪比隨噪音強度增加之值。 ... 40 Fig. 3-22 質量相異的粒子在相對應於外場頻率處的信號功率譜之峰值。 ... 41 Fig. 3-23 黏滯系數改變造成系統隨著噪音變化之訊號峰值位置，其所對應的溫度。 ... 42 Fig. 3-24 多個磁偶極以頭尾相接的方式排列之輸出訊號。 ... 43 Fig. 3-25 多個磁偶極以環狀之頭尾相接的方式排列之輸出訊號。 ... 44 Fig. 3-26 三個磁偶極呈環狀排列之示意圖。右圖為達平衡時三個磁偶極的角度。 ... 44 Fig. 3-27 四個磁偶極呈環狀排列之示意圖。右圖為達平衡時四個磁偶極的角度。 ... 45 Fig. 3-28 三個磁偶極以環狀排列時之穩態示意圖(順時針與逆時針排列)。 ... 45

(9)

Fig. 1-1 雙向環型雷射輸出的光訊號強度對噪音強度的關係圖(信噪比) (圖片擷 取自 McNamara et al. Theory of stochastic

resonance,physical reviewA,volume39)

Chapter 1 Introduction

1.1 隨機共振現象

Benzi 等人在研究古冰川問題時提出了隨機共振的概念。他們發現過去幾十萬年間，地球冰川期和暖期的交替週期恰與地球繞太陽轉動的離心率幾乎相同，意味著太陽對地球施加了週期變化的訊號。但是這個週期訊號很小，本身並不足以造成如此大幅度的變化。在他們的氣候模型中，指出地球處於非線性的條件下可能使地球可以處於冰川期與暖期兩個狀態，地球離心率的週期變化是能夠讓地球在這兩個狀態之間變化的外力，而地球所受的隨機力則大大的提高了小週期訊號對系統的調控機制。這種非線性系統在微弱的週期性訊號之下，由於隨機力的作用增加了系統隨訊號變化的週期現象，稱之為隨機共振。在其他的物理實驗文章裡也有提出看到了隨機共振的現象。第一篇是由 Fauve 等人在斯密特觸發器電路系統發現的，斯密特觸發器可以假設為有兩個穩態的輸出，在 Fauve 的論文中還不清楚斯密特觸發器系統為何會產生隨機共振此一結果；而在 McNamara[1]等人的文章中對此一系統提出了理論的解釋。另一個是由 McNamara 等人利用雙向環型雷射所做的光學實驗中完成的，它有順時針與逆時針兩個模態形成雙穩的性質，其輸出的光強度顯示出隨機共振的現象，如 Fig. 1-1 所示；McNamara 同樣也對此系統做理論解釋，但是實驗上有太多自由的參數而無法完全的相互吻合。

(10)

當系統中施加的外場很微弱時，我們會直觀的認為系統並不太會受到此外力之影響而有大動作，且當系統中有噪音的干擾時，若輸入之噪音強度越大，此時輸出訊號之噪音亦會越大。上面的想法在線性系統中是成立的，但在非線性系統中則會出現有趣的現象；在系統、訊號與噪音之間存在某種關係時，如果增加噪音反而會增加輸出訊號的強度，輸入噪音高於或低於這個強度時，輸出的訊號都會降下來。這一個現象類似於力學中所謂的共振，因此上述的現象被稱之為”隨機共振”。這種效應讓我們在解釋一些類似的物理圖像時，有另一種詮釋的空間，例如鳥類與昆蟲對飛行方向之感測是否和這種非平衡統計問題相關仍是個未知的題目。

1.2 隨機共振的物理圖像

上節所討論的隨機共振現象包含三個不可缺少的因素：具有雙穩或多穩態的非線性系統、輸入訊號及噪音。我們用一維雙穩態來說明這個物理圖像。考慮一個具有雙穩態的位能系統，如 Fig.1-2(a)所示，一個慣性很小的粒子落於其中一個位能井中。此時加入一個微弱的低頻週期訊號A0sin(t)，則系統的位能會隨著訊號而改變如 Fig.1-2 (b)；當t  /2時，輸入訊號同時抬升與降低了左右兩邊的位能高度，因而降低了左邊位能井的深度、增加了右邊位能井的深度，故位於不穩定點左邊的粒子有更高的機率能夠經過不穩定點到達另一個位能井；當 2 / 3 t 時則剛好相反，粒子會有從右邊的位能井越過不穩定點抵達左邊的趨勢；而在t0,時，因為此訊號同時抬高或者同時降低了位能井的高度，所以待在不穩定點兩邊的粒子並沒有因為外場的影響而對某一邊的位能井有所偏好。必頇注意的是這個週期性訊號的強度A₀較微弱，在未達到最終的穩定態時並不足以令粒子越過位能障壁到另一個穩態，而只能讓粒子在原本的位能井中做震動。當加入噪音時，會影響位能井之間躍遷的速率；一開始噪音的強度增加時，位能井之間的躍遷速率會增加，使得系統更容易隨著訊號而有週期性的改變。而當噪音大到某個程度之後，則因為位能井之間的躍遷變得太過容易而導致粒子在兩穩態

(11)

Fig.1-2 (a)無外場情形時的雙穩態位能圖；當系統中加入一弱外場A₀sin(t)，(b)

為當t不同的情況下位能之改變。其中位能井中的粒子若是在某個位能井中有較

(12)

Chapter 2 Theory of stochastic resonance

2.1 隨機共振理論

當輸入信號和噪音強度很小時，系統的吸引域可視為仍由系統的位能決定；在此吸引域指的是位能低點附近的區域，在這個區域中較容易有粒子停駐而會有較高 的粒子密度。若輸入信號的頻率 w 很小，又可以進一步的認為系統在各個吸引域 達到局部平衡所需要的時間遠小於兩吸引域之間機率達整體平衡的時間，亦即訊號變化的週期長到足以令系統在局部區域達到波茲曼分布。與信號變化和兩個穩定態吸引域間概率交換所需的時間相比，在各個吸引域內達到概率平衡可以視為是瞬間完成的，這一個近似圖像即是絕熱近似。在絕熱近似下，考慮一個粒子在位能為雙穩態的一維系統中運動，設當xx'時此系統的位能在一次微分為零的不穩定點，而xx_時擁有最低能量。將xx'的區域視為系統的正狀態、x x'的區域視為系統的負狀態，則粒子處於正或負狀態的機率分別為



       ' ' ) ( ) ( x x dx x n dx x n   2-1 其中(x)是系統位於 x 的機率密度。n 、_ n 隨時間改變的機率交換為                   n t R n t R n n t R n t R n ) ( ) ( ) ( ) ( 2-2 ) (t R 是離開狀態的速率。又由於n n 1，故上式可改寫為         n R t  R t R t n n ( ) [ ( ) ( )] 2-3 為了方便討論，我們設x'=0, x_ c。式 2.3 可解出為

(13)



( ') ( ')



'] exp[ ) ( ' ) ' ( ) ' ( ) ( ) ( ) ( 0 0 0 1



                  t t t t dt t R t R t g dt t g t R t n t g t n 2-4 在上式中，隨意的R(t)並不一定能夠做積分。我們假設R(t)有以下形式： ) cos ( ) (t f ₀ t R_   _s 2-5 對 0cosst做展開：





... cos ... cos cos 2 1 2 2 0 2 0 2 2 0 2 0 1 0           t R R t t R s s s             2-6 其中 _n n n n d f d n f     ! ) 1 ( 2 1 ), ( 2 1 0    。將上式展開取到的一次項，代入 2-4 式：











t t



t t t n t t dt e t t n t t t n t t t g s s s s s s s s t t t t s                         sin cos 1 2 1 2 1 sin cos 1 2 1 2 1 ) ( )] ( exp[ ' ' cos 2 1 ) ( )] ( exp[ ) ( )] ( exp[ ) ( 0 0 1 2 2 0 0 0 1 2 2 0 0 0 ) ' ( 0 1 0 0 0 0 0 0 0 0 0                                   



2-7 設 2 0 2 0 0 2 0 2 0 0 sin , cos             ，則





) cos( sin sin cos cos sin cos 1 0 2 2 0                   t ts t t t s s s s s s s 故                            _  ₂ ₂ 0 0 1 2 2 0 0 1 0 0 0 ) cos( 1 ) cos( 1 ) ( 2 )] ( exp[ 2 1 ) ( s s s st t t n t t t n              2-8 利用nn 1即可得到n (t)。 

(14)

當t₀ 時，漸進解為              _    ₂ ₂ 0 1cos( ) 1 2 1 ) ( lim ) ( 0 s s t s t t n t n      ，與t₀時刻的分布無關。得到了 2-8 式，我們可以計算系統的自相關系數為 ]} 1 ) , | ( 2 [ )} , | ( ] 1 ) , | ( 2 1 ) , | ( 2 {[ ) , | ( ) , | ( ) , | ( ) , | ( ) , | ( ) , | ( ) , | ( ) , | ( ) , ( ) , | , ( | ) ( ) ( 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 , 0                                        

 

t c t n t x t n t c t n t c t n c t x t n t c t n c t x t n t c t n c t x t n t c t n c t x t n t c t n c dxdy t x t x t y n xy t x t x t x           2-9 取t₀ ， ) ( 2 ]} 2 ) 2 ( cos[ {cos ) ( cos 1 | ) ( ) ( lim ) ( ) ( 2 2 0 2 0 2 1 2 2 2 0 2 2 0 2 1 | | 2 0 , 0 0 0 s s s s s t t c t e c t x t x t x t x t x                                         2-10 此一相關函數不僅與時間間隔 有關，也與起始時間 t 有關。真正在實驗中測量相關函數時，是從不同的起始時間開始而進行系綜平均，所以從統計意義上的相關函數應該是對時間做平均，即 ) ( 2 cos ) ( 2 1 ) ( ) ( 2 ) ( ) ( 2 2 0 2 0 2 1 2 2 2 0 2 0 2 1 | | 2 2 0 0 s s s s average c e c dt t x t x t x t x s                                 



2-11 系統輸出變量的功率譜密度S()為相關函數的傅立葉變換，即 ) ( ) ( ) ( 2 cos ) ( 2 1 ) ( ) ( ) ( 2 1 2 2 0 2 0 2 1 2 2 2 0 2 0 2 1 | | 2 0                          S S d e c d e e c d e t x t x S i average s s i s i average               



             〉〈 2-12 其中

(15)

)] ( ) ( [ ) ( 2 ) ( 2 cos ) ( ) ( 2 1 2 1 ) ( 2 1 1 ) ( 2 1 ) ( 2 1 ) ( 2 1 ) ( 2 1 ) ( 2 2 0 2 0 2 1 2 2 2 0 2 0 2 1 2 2 2 2 0 2 0 2 1 2 2 0 0 2 0 2 2 0 2 0 2 1 2 0 2 2 0 2 0 2 1 2 0 2 2 0 2 0 2 1 | | 2 0 2 2 0 2 0 2 1 | | 2 2 2 0 2 0 2 1 | | 2 1 0 0 0 s s s i s s s s s i s i s i s c d e c S c i c i c d e e c d e e c d e e c S                                                                                                                                    



                2-13 由上式可看出S₂()是來自於輸出信號，它是一個和輸入訊號同頻率呈現 分布的函數，而S₁()則是來自輸出的噪音，它是具有羅倫茲形式的連續分布頻譜。在 ) ( 2  S 中包括兩個對稱於=0 的 delta 函數，我們只取正的頻譜來討論，定義 one-sided 頻譜密度為兩倍的S₂()在的範圍。而輸出的信噪比定義為輸出總信號功率譜與_s處的平均功率譜之比值，即 1 2 2 0 2 0 2 1 0 2 0 2 1 2 2 0 2 0 2 1 2 2 0 0 2 2 2 0 2 0 2 1 2 1 0 2 ) ( 2 1 4 ) ( 2 1 4 ) ( ) ( ) (                        



s s s s s c c S d S R                        2-14 上面的討論中僅只是在絕熱近似中對躍遷速率做展開而得到相關函數的資訊，所以我們只要能夠得到系統的躍遷速率，並將躍遷速率寫成 2-5 的形式，就能夠得到功率譜密度與信噪比。

(16)

Fig. 2-1 逃逸問題示意圖。

2.2 逃逸問題

上節中我們討論到對於隨機共振所展現的功率譜可以由穩態之間的躍遷速率所決定，因此如何求得躍遷速率是接下來幾節的重點；而系統在平均某個時間點第一次經過穩態間的不穩定點，即平均首通時間，與躍遷速率以及描述系統方程的運算子(Fokker-Planck Operator)之本徵值存在著某種關聯，因此在這章我們將探討首通時間、躍遷速率及 Fokker-Planck Equation 之本徵值的問題。考慮在一維空間位能場U(x)中運動的粒子，其運動方程為 ) ( ' x U x x m      _2-15 上式中 m 為此粒子的質量，為阻尼係數， dx U d x U   '( ) 為此位能場對粒子的作用力。在過阻尼的作用下，慣性項mx的作用可忽略，故方程式變為 ) ( ) ( ' ) ( ' x C x U x U x      2-16 即粒子的運動決定於位能 U(x)的結構。Fig.2-1 表示一種常見的位能圖，該位能在 min x 處有極小值，在x_max處有極大值，分別對應到系統的穩定態與不穩定態。當 x 趨近時，U(x)分別趨近 ，系統在 xx_max 是被束縛的，稱為束縛區；而在xxmax是不受束縛的，稱為逃逸區。考慮噪音的影響，方程式 2-16 可寫為 ) ( ) (x t C x   2-17 我們假設隨機力為白噪音，即 ) ' ( 2 ) (t' (t) , 0 (t)     D tt  〉〈〉  〈 2-18 逃逸問題是指系統從穩定態x_min處出發，在除了 由位能 U(x)所產生的力之外的其他力(例如噪 音)作用下，逃離位能井的演化問題。接下來我

(17)

們探討的是在微弱的外場作用下，系統由於噪音的作用而產生逃逸的現象。

2.2.1 克萊默斯逃逸速率

我們將計算 Fig.2-1 的位能圖中，首先位於xx_min附近的粒子越過xx_max 的逃逸

速率。對應於方程式 2-17 的 Fokker-Planck Equation 為 ) , ( )] , ( ) ( [ ) , ( 2 2 t x x D t x x C x t t x _ _           2-19 在絕熱近似下，系統達到局部平衡所需的時間遠小於整體平衡所需的時間，此局部平衡的狀態稱為準穩態。定義 Fokker-Planck 算子為 2 2 ) ( x D x C x L_FP        2-20 在準穩態時， ( , )0   t t x  ，即 0 ) , (x t  L_FP 對 x 積分得 S t x e x De D x U D x U     ) , ( ) ( ) (  2-21 其中 S 為通過x_max處由左往右的機率流。將上式兩邊同乘上 D x U e ) ( ，並對 x 從x_min積 分到圖中的點 A(Ax_max )，得到 dx e S t A e t x e D A x D x U D A U D x U



  min min ( ) ( ) min ) ( )] , ( ) , ( [   2-22 因為粒子處於 A 點的時間很短，會很快的跑向的區域且沒有來自束縛區的粒子補充，所以在 A 點的機率密度(A,t)可以視為零，故上式變為



 A x D x U D x U e t x De S min min ( ) min ) ( / ) , (  2-23 假設能量差U很大，則在 x=x 附近的機率密度可近似成

(18)

D x U x U e t x t x ) ( ) ( min min ) , ( ) , (     則在x_min左右鄰近的區間{x₁,x₂}內找到粒子的機率 p 為 dx e t x dx t x p x x D x U x U x x



    2 1 min 2 1 ) ( ) ( min, ) ( ) , (   2-24 由於機率 p 乘上逃逸速率 r 會等於機率流 S，故從式 2-23 和 2-24 可得到



     A x D x U x x D x U A x D x U D x U e e D p S r e t x De S rp min 2 1 min min ) ( ) ( ) ( min ) ( / ) , (  2-25 上式分母中的第一個積分主要貢獻來自於x_min附近的區域，第二個積分主要來自於 max x 附近處，所以可以將 U(x)近似於這兩處做展開得 2 max max max 2 min min min ) ( | ) ( '' | 2 1 ) ( ) ( ) )( ( '' 2 1 ) ( ) ( x x x U x U x U x x x U x U x U       2-26 將 2-26 帶回 2-25，可以得到 r D x U x U e x U x U r ) ( ) ( max min min max | ) ( '' | ) ( '' 2 1     2-27 r 表示機率密度由穩定區流入不穩定區的速率，被稱為克萊默斯逃逸速率，由系統 的位能差與位能在微分為零點處的曲率所決定。

2.2.2 首通時間

首通時間定義為在一維位能場中，粒子由x₁ x' x ₂出發，在 T 時間首次穿過 x、₁ x₂ 所需的時間。 設 P(x,t|x’,0)為初始處於x'的粒子在

(19)

t 時刻到達x₁ xx₂的機率密度分布，它遵循 Fokker-Planck Equation： ) ' ( ) 0 , ' | 0 , ( ) ( x x x x x L t FP         

2-28

並考慮吸收邊界條件為(x1,t|x',0)(x2,t|x',0)0，即當粒子一到達x、1 x2則立 即除去不考慮，類似於在逃逸過程中機率一旦進入不穩定區後便不再返回。W(x’,t) 為 t 時刻處於區間{x₁,x₂}內的機率總量



 2 1 ) 0 , ' | , ( ) , ' ( x x dx x t x t x W  2-29 在 T→T+dt 時間內粒子離開{x₁,x₂}的機率(即粒子具有首通時間 T→T+dT 的機率) 為  



 2 1 ) 0 , ' | , ( ) , ' ( x x dxdT x T x T x dW  2-30 則首通時間 T 的分布函數為w(x',T)



     2 1 ) 0 , ' | , ( ) , ' ( ) , ' ( x x dx x T x dT T x dW T x w  2-31 利用上式可算出首通時間分布的各階矩為



   0 2 1 ) ' , ( ) , ' ( ) ' ( x x n n n x T w x T dT P x x dx T 2-32



 _   0 ) 0 , ' | , ( ) ' , (x x T x T x dT Pn n 2-33 對 2-33 式做部分積分，可得到 1 , ) 0 , ' | , ( ) ' , ( 0 1  

_

   n dT x T x T n x x P_n n  2-34 其中的 ( , ') ( , | ',0) ( ,0| ',0) ( ') 0 0 x x x T x dT x x x x P 



        ，因為在t時所有的粒子都會離開邊界，使得在區域(x1,x2)內的機率密度分布為零。

(20)

將L_FP算子作用在 2-34 式，可以得到多個微分方程： 1 , ) ' , ( ) ' , ( ) (x P x x nP_₁ x x n L_FP _n _n 2-35 即 ... ) ' , ( 2 ) ' , ( ) ( ) ' ( ) ' , ( ) ( 1 2 1 x x P x x P x L x x x x P x L FP FP       知道了 2-35 的各階矩，就能夠得到首通時間的機率分布和各種統計性質。n1所代表的是首通時間的期望值，或稱平均首通時間。對於具有 Fig.2-1 位能形式的逃 逸問題而言，粒子的束縛區只在正方向有吸收邊界，x=A 可視為吸收壁。將初始 點 x’選為穩定點x_min，平均首通時間為 dx x x P x T



   0 ₁ _min 1( ') ( , ) 2-36 利用 2-35 的一階矩求解得P ，再帶入 2-36 則可得到平均首通時間： ₁ dx e dx e D x T dx e e D P A x D x U A D x U A x D x U D x U



      min min ) ( ) ( min 1 ) ( ) ( 1 1 ) ( 1 2-37 比較上節 2-25，可得到 r x T1( min)1，即首通時間與機率逃逸速率之間有密切的關係。在推導首通時間的時候我們只要求 A 處存在吸收壁，不一定要求噪音強度要比位能差小很多，這說明T 的計算結果比上一節的計算結果來得廣泛。即使不符₁ 合上一節的弱噪音條件，逃逸速率也可以 2-25 的同樣形式表示。

2.2.3 Fokker-Planck Equation 的本徵值

給定 Fokker-Planck Equation ( , ) L (x,t) t t x FP     ，其中 ₂ 2 ) ( x D x C x LFP        我們可以寫下L_FP算子的一組本徵值和本徵向量

(21)

將會看到式中的_n皆大於零。由於L_FP取共軛之後的L*_FP  L_FP，所以它並不是厄米算子；在量子力學中對於厄米算子已經有很充分的分析，我們可以通過某些變換來建立L_FP算子的本徵空間與厄米算子的本徵空間之關聯，透過這些連結後可將厄米算子的特性套用至經過變換後的LFP中。可用這個變換來達成： ) , ( ) , ( 2 ) ( t x e t x D x     2-39 將 2-39 帶回 Fokker-Planck Equation，可以得到 ) ( ' 2 1 ) ( 4 1 ) , ( ) , ( 2 2 2 x C D x C x D L t x L t t x s s         _ 2-40 這裡L 是厄米算子，透過 2-39 式的對應使得_s L 和_s L_FP有完全相同的本徵值。由於厄米算子的本徵值為實數，所以一維LFP算子的本徵值也均是實數。因



           [ ( )] 0 ) ( ) ( 2dx x x f D dx x f L x f _FP 2-41 所以L_FP的本徵值非正，它的最大本徵值為₀ 0,對應到的定態解| f₀〉(x)，而其他所有的本徵值皆為負、1 2...，接下來的討論我們按絕對值的大小依次排列本徵值的序號，即₁₂ ₃...。 FP L 的各個本徵值連繫著 Fokker-Planck Equation 的各種不同的時間尺度演化的過程。我們發現對逃逸問題而言₀不存在，因為此系統並不存在穩定態解。當噪音很小時，逃逸過程反映了 Fokker-Planck Equation 中最緩慢的時間過程，它應與最大時間尺度的本徵值₁有關。接著我們計算₁。 Fokker-Planck Equation 的本徵方程為 1 1 1 ) ( ) ( 1 e f f x e x D f L D x U D x U FP _       2-42 對 2-42 式從積分到 x，可得



      x D x U D x U y f dy e D f e x 1( ) / ) ( 1 1 / ) (  2-43

(22)

再對 2-43 式從x_min積分到 x，得到 ] ) ( ) ( [ ) ( ₁ ) ( 1 min 1 ) ( ) ( 1 min min



     x x x D x U D x U D x U y f dy e dx D x f e e x f  2-44 由於₁反映 Fokker-Planck Equation 緩慢變化的時間尺度，做為零級近似，我們取1 0代入 2-44，得 ) ( ) ( ₁ _min ) ( ) ( 1 min x f e e x f D x U D x U   2-45 再將 2-45 代入 2-44 並利用 f1(A)0，可以得到







    A x D y U D x U e dy e dx D min 1 ) ( ) ( 1  移項即可得到



     _A x D y U D x U e dy e dx D min ) ( ) ( 1  2-46 上式正是 2-25 之結果，即 Fokker-Planck Equation 不為零的最大本徵值是躍遷速率的倒數。

2.3 一維雙穩態之間的概率躍遷

具有 Fig.2-2 之雙穩態系統，和前面討論的逃逸問題不同的是在x處機率密度並不會到無窮大，而會滿足歸一化的條件。我們將要探討從某一個穩態出發到最終定態的機率密度變化，2.2.1 節的計算方法與結果仍可稍做改變用在雙穩態的情形，但在這節裡我們採用另一種方法來討論機率密度演化的問題，可以更簡易的推廣到多維多穩態系統，並且有更清楚的動態圖像。從准穩態發展到最終定態的過程即是兩個位能井之間進行機率躍遷交換的過程，此系 Fig. 2-2

(23)

統的 Fokker-Planck Equation 仍然用 2-19 來表示，此時系統位能在 1 s x 和 2 s x 兩處有兩個極小值，而x 為位能的極大值。在弱噪音的情況下，系統符合絕熱近似，即系_u 統從不穩定態到準穩態所需要的時間遠小於從准穩態達到整體穩定態所需的時間，所以兩個位能井內的機率密度分布( tx, )可以看成總是符合定態分布：                 u D x U u D x U x x e t N t x x x e t N t x t x , ) ( ) , ( , ) ( ) , ( ) , ( ) ( ) (    ， 2-47 當整體的最終穩定態尚未達到時，N_(t)、N_(t)會變化，但分布函數保持不變。可以分別在 1 s x 和 2 s x 採用高斯近似，得到                     u D x x x U D x U u D x x x U D x U x x e e t N t x x x e e t N t x t x S S S S S S , ) ( ) , ( , ) ( ) , ( ) , ( ₂ 2 2 2 2 1 1 1 ) )( ( '' 2 1 ) ( ) )( ( '' 2 1 ) (    ， 2-48 設p 分別為_ x 兩邊在 t 時刻所含的機率總量， u ) ( ' ' 2 ) ( ) , ( ) ( ) ( ' ' 2 ) ( ) , ( ) ( 2 2 1 1 ) ( ) ( S D x U x S D x U x x U D e t N dx t x t p x U D e t N dx t x t p S u S u               



    2-49 對 2-19 式分別在x 的左右兩邊做積分，並利用機率密度分布函數在_u x處約等於零，可以得到 u x x t x x D dt t dp      ( , )| ) ( _  2-50 從上式可以看出p 的變化率只與_ ( tx, )在x 附近的性質有關。 _u 我們將從 2-19 出發討論機率密度的流動。由於機率密度的改變發生在不穩定點 u x x 處，可以利用在不穩定點線性化取到 x 的一次項，即讓漂移項變成線性項，

使得 Fokker-Planck Equation 可寫成 Ornstein-Uhlenbeck process 的形式(由於

Fokker-Planck Equation 的解並不是很容易可得到，Ornstein-Uhlenbeck process 是少

(24)

) , ( )] , ( ) ( ' ' [ ) , ( ... ) ( ' ' ) ( ' ) ( ' 2 2 t x x D t x x U x t t x x xU x U x U u u u               2-51 我們可以利用將上式做傅立葉轉換至 k 空間，這將可使原來的二階微分方程變成 一階微分方程，從而解出( tx, )。首先若假設初始條件呈 函數的分布，則 Ornstein-Uhlenbeck process 的解可精確的解出為一個高斯分布函數： ) 1 ( 2 ) ' )( ( '' ) ' )( ( '' 2 ) ' )( ( '' 2 2 ) ' )( ( '' ) 1 ( 2 ) ( ' ' ) ' , ' | , ( U xu t t t t u x U u u e D x e x x U t t x U u e e D x U t x t x             _ 2-52 它描述了初始條件為 函數分布而發展到成高斯函數分布的機率密度變化(是一個

transition probability density)。考慮初始的分布在xxu中，即討論從 p(t)向p(t)的

流動。因為絕熱近似的關係，可以假設系統的初始分布函數為 D x U x U S u u D x x U in u S u e D x U k x x x x e t p k t x t x ) ( ) ( 2 ' ) ( '' 1 1 2 2 ) ( '' ' , 0 ' , ) ' ( ) ' , ' ( ) ' , ' (                  2-53 則由此一初始分布出發，在 t 時刻的機率密度函數為



   0 ' ) ' , ' ( ) ' , ' | , ( ) , (x t  x t x t in x t dx  _ 2-54 在條件 2 ''( )( ') 1 1     _U _x _t _t u e D 時，略出高於 x 二次方的項則可積出為 ] 2 ) ( ' ' 2 1 [ ) ' ( ) , ( x D x U k t p t x u    _ _  2-55 由上式可看出其解與時間間隔無關，即由初始條件達到局部穩定之後其機率密度的分布函數即和t無關。將 2-55 帶入 2-50 ) ( | ) ( '' | ) ( '' 2 ) ( 2 ) ( '' | ) , ( ) ( ) ( ) ( 1 1 U x e p t x U D t p D x U Dk t x x D dt t dp D x U x U u s u x x u S u               2-56 由於 dt t dp_( ) 可以寫成

(25)

) ( ) ( t p R dt t dp    __ _2-57 比較 2-56 與 2-57 兩式，可以得到從xx_u躍遷到xx_u的機率R ， _ D x U x U u s u S e x U x U R ) ( ) ( ₁ 1)| ''( )| ( '' 2 1    _ 2-58 同樣的若將初始條件設在x'xu處，經過上述的運算同樣的可求得從xxu要遷到 u x x 的機率R ， _ D x U x U u s u S e x U x U R ) ( ) ( 2 2)| ''( )| ( '' 2 1    _ 2-59 較一般化的機率密度躍遷過程可寫為 ) ( ) ( ) ( ) ( ) ( ) ( t p R t p R dt t dp t p R t p R dt t dp                 2-60 所以從准穩態出發到達穩定態的機率密度隨時間的變化完全由 2-58、2-59 與 2-60 所決定。而上式p 隨時間變化之間的作用矩陣 _             R R R R 有兩個本徵值：₀ 0,₁ R_R_。₀對應到最終定態的本徵值，₁反映了機率密度交換的最長時間尺度演化，在雙穩問題中，平均首通時間為    R R T₁ 1 2-61

2.4 多維多穩態之間的概率躍遷

在 2.3 節中所使用的方法關注的是在位能井之間的不穩定點處的性質，並不依賴於維度和穩態的數目，因此我們可以將上面的討論做推廣。首先考慮二維的雙穩態 系統：一個二維的系統具有位能U (x,y)，此系統忽略慣性項時的牛頓方程可寫成

(26)

                   ) , ( ) , ( 1 ) , ( ) , ( 1 y x U y y x U y y x U x y x U x y x   2-62 其對應到二維的 Fokker-Planck Equation 為 ) , , ( ) ( )] , , ( ) , ( [ )] , , ( ) , ( [ ) , , ( 2 2 2 2 t y x y x D t y x y x U y t y x y x U x t t y x y x                    2-63 假設此二維的位能面在 x>0 與 x<0 各有一個最低點(位能井)，原點處在 y 方向是一 吸引域而在 x 方向則是排斥域，即(0,0)是一鞍點；這兩個位能井之間的機率密度流 動必然在原點附近最為關鍵，而鞍點處亦是在某一方向的局域性極小值，所以可以藉由在鞍點處線性化使得變成 Ornstein-Uhlenbeck process 的形式： ) , , ( ) ( )] , , ( [ )] , , ( [ ) , , ( 2 2 2 2 2 1 x y t y x D t y x y y t y x x x t t y x _ _ _ _ _                 2-64 其中、₁ ₂皆大於零： ) 0 , 0 ( 2 2 2 ) 0 , 0 ( 2 2 1 | ) , ( 2 1 | ) , ( 2 1            y x y x y x U y y x U x   2-65 當初始分布呈 函數的分布，則可以精確的解出為 ) 1 ( 2 ) ( ) 1 ( 2 ) ( ) ( 2 ) ( 2 2 1 0 0 0 ) 0 ( 1 2 2 ) ( 2 0 2 ) 0 ( 1 2 2 ) ( 1 0 1 0 2 0 1 ₁₎₍₁ ₎ ( 2 ) , , | , , (                 tt o t t t t o t t e D e y y e D e x x t t t t e e e e D t y x t y x     _        2-66 現在我們要考慮從 x<0 處出發要遷至另一個位能井的機率密度。假定初始分布在 x<0 的位能井內，即         0 , 0 0 , ) ( ) , , ( ) ( 0 0 x x e t N t y x D x U ，  2-67 在原點展開到二次項得到

(27)

         0 , 0 0 , ) ( ) , , ( 2 2 2 1 ) 0 , 0 ( 0 0 x x e e t N t y x D y x D U    ， 2-68 將此初始條件對 2-66 積分，並且利用近似條件 ₀ 1 1   t t eD 並略去 x 平方以上 的高次項，可得到 D xx x D y D U e dx e e t N t y x t y x dy dx t y x 2 2 0 0 2 ) 0 , 0 ( 0 0 0 0 0 0 0 0 1 2 0 2 2 1 ) ( ) , , ( ) , , ( ) , , (                 



  2-69 對方程式 2-63 在x(,0)、y(,)區域進行積分，可得



               (0, ) (0, , ) ( , , ) | ₀ ) ( x x dy x y t dx d D t y y dyU dt t dp _ _ 2-70 將 2-66 帶入上式，由於第一項的貢獻比第二項小很多故可忽略，得到 D e t N dt t dp U_D 2 1 ) 0 , 0 ( ) ( ) (      2-71 在弱噪音的條件下，在 x<0 的位能井中的分布可以近似為高斯分布[3]， ) ( ) ( ) ( 2 ] 2 ) ( ) ( exp[ ) ( ) ( 2 1 ) ( 2 2 2 1 ) (                  



     D U D U De t N D w s dw ds e t N t p 2-72 在這裡₁_,₂()是在 x<0 的位能井處之矩陣                 ) ( 2 1 ) ( ) ( ) ( 2 1 yy xy xy xx U U U U 2-73 的兩個本徵值。s 和 w 為相對應的 U(x,y)的二次部分在位能井底的兩個主軸方向。 比較 2-71 和 2-72，可得到 ) ( ) ( t p R dt t dp   

(28)

D U U e R ) 0 , 0 ( ) ( 2 1 2 1( ) ( ) 2 1      _     2-74 用同樣的方法，假設初始條件都在 x>0 處，則可以算出由 x>0 躍遷到 x<0 的速率 D U U e R ) 0 , 0 ( ) ( 2 1 2 1( ) ( ) 2 1      _     2-75 將以上二維雙穩態的問題推廣到二維多穩態的問題：假設 n i y x U y x U_x( _i, _i) _y( _i, _i)0, 1,2,3... 2-76 是確定性方程 2-62 的ｎ個穩定點，以這些穩定點為中心的吸引域是A₁,A₂,A₃...，D_ij 是第 i 和第 j 個吸引域之鞍點。p(Ai,t)為ｔ時刻在A 吸引域內的機率，則i p(Ai,t)隨時間的演化應由以下的方程所確定： ) , ( ) , ( ) , ( 1 1 t A p R t A p R dt t A dp i n j ji i n j ij i



     2-77 其中             } ) ( ) , ( exp{ ) ( ) ( ) ( ) ( 2 1 , , 0 2 1 2 1 D D U y x U D D i i A A R _i _i _ij ij ij j i ij      不相鄰當 2-78 式中₁(i),₂(i)為漂移項在(x_i,y_i)處線性化矩陣的兩個本徵值，₁_,₂(D_ij)為漂移項在鞍點處線性化矩陣之正負本徵值的絕對值。在弱噪音的條件下，q 維多穩態系統的演化也可用 2-76 式來描述，這時候的躍遷速率修正為 ( ) ( 1)exp{ ( ) ( )} 2 1 2 2 1 1 D D U A U A R i ij q v i v q v ij  



      。

(29)

Chapter 3 Magnetic dipoles in two-dimensional

system with weak field

3.1 二維磁偶極系統的能量結構

考慮一個由兩個磁偶極組成的系統，以頭尾相接的方式沿著 z 方向排成一直線(如 Fig.3-1 所示)。設兩磁偶極的偶極矩分別為 1 m 、m ，兩者間的距離為 r，則此兩磁偶極的交互作用能量₂ )] ˆ )( ˆ ( 3 [ 1 4 3 1 2 1 2 0 r m r m m m r U        。設 I 為磁偶極的轉動慣量， 則系統所對應的 Hamilton 為 U I I       2 2 2 1 ( ) 2 1 ) ( 2 1 _ _ Fig. 3-1 兩磁偶極以頭尾相接的方式沿著 z 軸排列。

(30)

系統的交互能量如 Fig.3-2 所示，其中₁與₂分別為m 、₁ m 與₂  r 之間的夾角，能 量U 在二維角度平面上有連續重複的關係。將上圖的₁=0 與₁=360 度、相接₂=0 與₂=360 度相接，這就類似一個 torus 的幾何型態；可以看到在角度(0,0)、 (180,180)、(0,360)、(360,0)及(360,360)處有能量最低點，而(0,0)、(0,360)、(360,0) 及(360,360)所表示的是同一個位能井，所以此二維的兩個磁偶極系統為雙位能井的系統，而考慮系統狀態由一個位能井躍遷到另一個位能井所走的重要路徑有兩條：經過鞍點(90,270)與(270,90)、以及經過鞍點(90,90)與(270,270)。 當系統中加入一隨時間改變大小以及沿著 z 兩方向改變的磁場 B 時，會使原本的能量面隨著磁場變化而改變，Hamilton 變成'H Bm₁Bm₂。 Fig. 3-3 系統的能量面在隨時間變的外場 B=3sin(wt)下之變化。

(31)

23 如 Fig.3-3 所示，我們可以看到在外場為 B=3sin(wt)時，系統的能量面會隨之有明顯的變化，這是因為外場的強度夠大導致完全的改變了穩態的位置。對應於'之牛頓方程為 2 , 1 2 , 1 ; 2 , 1 , ' _ __      q j q H p j j 3-1 當系統中有噪音(t)的影響並考慮耗散項 1,2   ，3-1 之牛頓方程會變成二維的 Langevin 方程： 2 , 1 2 , 1 ; 2 , 1 ), ( ' _               q j t q H I j j j 3-2 為了簡化問題，在此考慮(t)為白噪音，即(t)符合 2-18 式。在沒有噪音的情況下，我們選取強度不足以改變系統的穩態結構且無法驅動系統在兩個最低的能量態中做躍遷的外場，利用式 3-2 模擬當噪音強度增強時系統躍遷的變化。

3.2 Molecular dynamics simulation

我們模擬上述系統隨著外場變化的動態過程，所使用的編譯器為 DevC++。為了之後能夠將此系統應用在模擬鳥類或昆蟲飛行之定向，在模擬的過程中所選用的外場是為赤道附近的地磁強度3105G，磁矩間的距離參考 Solov’yov 對鴿子喙部磁性物質的排列，而實際的微觀生物系統中的質量(或轉動慣量)確實是遠小於黏滯系數，符合理論中的假設。我們希望利用接近實際系統的參數來得到系統隨外場與噪音影響下變化的情形。在模擬的過程中，我們是取 100 組與 z 軸夾角不同的兩個磁偶極做為初始條件(50 個處於第一個狀態、50 個處於另一個狀態)，先經過一段時間達到局部的平衡後，在不同的噪音之下觀察這 100 個狀態受到隨時間改變的弱磁場作用，分別處在兩個狀態的系統個數如何變化，再對輸出的系統個數訊號做討論。在這裡所謂的兩個狀態指的是從鞍點及臨線劃分出的兩個能量狀態，其中最低點的能量對應到系統的兩個極小值(Fig. 3-2 的(0,0)及(180,180))。我們將系統簡化視為具有雙穩態的位

(32)

能，將整個角度的範圍畫分成由這兩個狀態所概括(如 Fig.3-4 所示)。 Fig. 3-4 系統兩個狀態的劃分示意圖。接著我們簡單的介紹模擬的方法與細節，並且討論得到的結果。

3.2.1 Verlet algorithm

Verlet algorithm 是將前一刻與下一刻的位置向量 r(t)在 t 時刻做泰勒展開至Δt 的三 次方項： ) ( ) ( ! 3 1 ) ( 2 1 ) ( ) ( ) ( ) ( ) ( ! 3 1 ) ( 2 1 ) ( ) ( ) ( 4 3 2 4 3 2 t O t t b t t a t t v t r t t r t O t t b t t a t t v t r t t r                       將上面兩個式子相加可以得出 ) ( ) ( ) ( ) ( 2 ) (t t r t r t t a t t2 O t4 r         3-3 其中 a(t)由牛頓力學的結果可以得到   m t a( ) 1 U(r(t))，U(r(t))為系統的位能；而 b(t)是位置對時間的三次微分項。由可以看到 r(t)的誤差值會是t4的數量級，因此對於計算位置向量是一種準確的方法。但是在計算速度的時候，由於

(33)

t t t r t t r t v        2 ) ( ) ( ) ( 3-4 我們會發現此時v(t)的誤差會變成 2 t  的數量級，要減少速度的誤差值最直觀的方 法變是將時間間隔Δt 減小，不過這樣子將會增加不少的運算量。為此許多人開始 對 Verlet algorithm 做改良以期能解決速度誤差值的問題。而在我們的模擬過程中 我所採用的是 velocity Verlet scheme，在交錯的時間點上計算位置與速度的值，。 其演算的程序如下： t t t a t t v t t v t t r U m t t a t t a t v t t v t t a t t v t r t t r                             ) ( 2 1 ) 2 ( ) ( )) ( ( 1 ) ( ) ( 2 1 ) ( ) 2 ( ) ( 2 1 ) ( ) ( ) ( 2 3-5 另外由於我所模擬的系統是一個磁偶極在二維空間裡的轉動過程，在這裡的位置向量會被磁偶極轉動的角度所取代、速度向量變為角速度、加速度是為角加速度進行運算與討論。

3.2.2 Langevin Thermostat

由於系統中含有耗散項，因此系統中的能量並不守恒。由於系統中的粒子並不會逃離此系統，我們可以說這是一個正則系綜(canonical ensemble)。正則系綜有溫度固定的性質，要實現定溫有很多種不同的模擬方式，我們所選用的是 Langevin Thermostat 的控溫方法。它的概念是從 Langevin 方程出發的，考慮一個在熱浴中的粒子，在此圖像裡我們假設此粒子的體積遠大於組成熱浴的小粒子；當大粒子朝某個方向移動時，小粒子會產生一個整體的耗散力v阻止大粒子前進，而小粒子本身也具有動能故會 在每個時刻以不固定的方向碰撞大粒子，此即隨機力 f’， Langevin 方程可寫成

(34)

3-6 耗散力的作用會減少大粒子的動能，而施於大粒子的隨機力則會增加大粒子的動能；藉由耗散力與隨機力的作用，可以使得系統的平均動能符合系統的溫度，這樣子的控溫方式稱為 Langevin Thermostat。系統達平衡時的動能與溫度之關係可以用來檢驗模擬的正確性，由式 3-3 的演算法看到加速度為 ( ) 1 U(r(t t)) m t t a     ，意味著當質量越小的時候加速度越大，而影響速度與位移的參數除了加速度之外還有t，當t越小時可以減少加速度的效應，這時候t與質量的選取便直接的影響到系統的演化，可以藉由溫度來規範動能以確定模擬系統的正確性。由於一個磁偶極具有一個轉動的自由度，滿足 I kT 2 1 2 1 2   〉〈  3-7 在沒有外場的情況下，當系統達到平衡時應滿足式 3-5；又因為我們假設系統變化的速率很慢符合絕熱近似，因此可利用系統動能的變化情形來判斷達平衡所需的步數。如 Fig.3-5 所示，當系統質量為 0.01 單位質量、T=1 時，t=0.01 即滿足 3-5 式；且由圖中可知當系統所經步數大於 14000 步即達平衡。利用動能與溫度之規範關係，我們可以選取合適的步長t以及步數來達到絕熱近似的條件。

(35)

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x 104 0 1000 2000 3000 4000 5000 6000 7000 8000 time(*0.01sec) k in e ti c e n e rg y (J )

3.3 路徑積分與最可能路徑

由 Fig.3-2 能量面的圖我們可以看到，粒子原本處於某一個穩定態欲躍遷到另一個穩定態時，所走的路徑主要有兩條，即兩個鞍點的附近(藍色部分與黃色部分)。由於粒子的平均首通時間是為躍遷速率的倒數，所以我們可以藉由計算躍遷速率進而得到粒子經過此兩個鞍點個別的平均首通時間。由第 2.4 節中我們得到了多維多穩態的躍遷速率，為方程式 D U U e R ) 0 , 0 ( ) ( 2 1 2 1( ) ( ) 2 1      _     與 D U U e R ) 0 , 0 ( ) ( 2 1 2 1( ) ( ) 2 1      _     其中1,2()是在兩個位能井底部的線性化矩陣之本徵值，而1,2在位能的鞍點處做線性化所得到的。之前理論的推導中並未涉及位能面隨時間而變化，所得到的 ) ( 2 , 1   和₁_,₂與時間無關，即在位能的局部高點與局部低點處的曲率並不隨時間而 Fig. 3-5 達平衡時動能隨時間的變化。 1.34 1.36 1.38 1.4 1.42 1.44 1.46 x 104 0 1 2 3 4 5 6 time(*0.01sec) ki n e ti c e n e rg y( J)

(36)

外場而改變(如 Fig.3-3 所示)，這時候很明顯的會看到位能面各處的曲率會隨著時間改變。因此我們利用 2-65 與 2-73 求出隨時間變的1,2()和1,2後，再詴著將躍遷速率寫成 2-5 的形式。當位能隨外場B0coswst變化時，每個時刻位能井的曲率皆不相同，又因為兩個位能井對稱於鞍點，其本徵值皆為與 ) cos cos 10 252 . 2 10 107 . 2 2 ( 2 1 ) cos cos 10 252 . 2 10 107 . 2 2 ( 2 1 0 0 15 8 0 0 15 8 t w B t w B t w B t w B s s s s               3-8 粒子由一個位能井躍遷到另一個位能井有可能經過的兩條路徑分開討論如下： A. 粒子行經最可能路徑，即經過鞍點(0.5Pi, 1.5Pi)附近到達另一個位能井處。躍遷機率為 ) cos 10 252 . 2 10 107 . 2 cos 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( 138 . 0 ) cos 10 252 . 2 10 107 . 2 cos 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( 138 . 0 0 15 8 0 0 15 8 0 cos 2 3 2 0 15 8 0 0 15 8 0 cos 2 3 1 0 0 t w B t w B t w B t w B e R t w B t w B t w B t w B e R s s s s D t w B s s s s D t w B s s                                   3-9 其中B₀ 0.003為外場強度的振幅。由 3-7 式我們看到 Exp 與根號裡的項皆包含隨

時間改變的項，將根號寫成 Exp 的形式，即 ...Exp[ln(…)]並且與前面的 Exp 項

合併得到 Exp[x]，我們會發現當B0 0.003時 x 隨外場B0coswst擁有近乎線性的關

(37)

0.003 0.002 0.001 0.001 0.002 0.003 0.048 0.049 0.050 0.051 0.052 x Fig. 3-6 D=1 時 x 之圖形。 可得到 x=a(B₀cosw_st)+c，因此躍遷速率 ) 0756 . 2 201446 . 0 ( cos 548 . 2 342457 . 0 0 2 1 2 1 ) (        D t w B D x s e e e t R   3-10 有 f(₀cos_st)之形式。利用 3-10 式可以得到躍遷速率隨時間與噪音強度 D 值之改變，如 Fig. 3-7 所示，將其倒數即可得到經過鞍點的平均首通時間(Fig. 3-8)。 Fig. 3-7 粒子經過最可能路徑到達另一個位能井之躍遷速率隨著時間與噪音強度

(38)

Fig. 3-8 粒子經過最可能路徑達到另一個位能井之平均首通時間。 B. 粒子行經第二路徑，即經過鞍點(0, Pi)附近到達另一個位能井處。躍遷機率為 )) ) cos ( 25 . 0 2 2 )( ) cos ( 25 . 0 2 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( 4 1 )) ) cos ( 25 . 0 2 2 )( ) cos ( 25 . 0 2 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( ) cos 10 252 . 2 10 107 . 2 cos 2 ( 4 1 2 0 2 0 0 15 8 0 0 15 8 0 cos 2 3 2 2 0 2 0 0 15 8 0 0 15 8 0 cos 2 3 1 0 0 t w B t w B t w B t w B t w B t w B e R t w B t w B t w B t w B t w B t w B e R s s s s s s D t w B s s s s s s D t w B s s                                             3-11 同樣的，將 Exp 項與後面的根號項合併成 Exp[x]，我們會發現當B₀ 0.003時 x 亦可寫成線性關係，如 Fig.3-9 所示。

(39)

0.003 0.002 0.001 0.001 0.002 0.003 1.900 1.895 x Fig. 3-9 D=1 時 x 之圖形。 改變不同的 D 值可以得到躍遷速率為 ) 39297 . 1 1033 . 0 ( cos 23268 . 1 192667 . 0 0 0 4 1 ) cos ( ) (         D t w B D s s e e t f t R     3-12 其相對應的平均首通時間如 Fig. 3-10 所示。 Fig. 3-10 粒子經過另一路徑之平均首通時間。

(40)

0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 noise intensity(dB) tr a n s it io n r a te (# /s e c ) path1 path2 Fig. 3-11 B0coswst=0.003 時，兩條路徑躍遷速率隨著噪音強度之變化。比較兩條路徑之平均首通時間，粒子經過最可能路徑到達另一個位能井所需的時間尺度較小，而經過另一個可能路徑所需的時間尺度較大，所以在 D 非常小的時候粒子會走最可能的路徑；而在 D 稍大後粒子才可能會走第二條路徑。在這裡必頇說明，以上對於第一、第二路徑之討論是在考慮只有兩個吸收邊界下的結果；上述當粒子越過第二條路徑到達另一個位能井時，原則上應該有四個吸收邊界、即此時亦有行經第一路徑的可能。因在此主要討論第一路徑，所以沒有考量當噪音強度 D 較高時之第二路徑之詳細結果；而當第二路徑發生時，和第一路徑耦合所產生之躍遷速率與首通時間該如何表示，這可以在之後做更深入的探討。

3.4 功率譜密度

經由上述模擬方法可得到在弱外場之下每個時刻處於某個狀態的系統個數，將這些數據求取自相關函數並做傅立葉轉換後，我們就可以得知系統資訊處於某個頻率下之強度，而上述的概念便是功率譜密度S()的結果。在實踐由數據求得功率

(41)

設現在有 N 個訊號點： ) 1 ( ) . . . . 1 ( ), 0 ( X X N X ，將這 N 個訊號分割成 k 組，每一組有 L 個訊號點，且每一組訊號點會有 D 個點產生重疊，如 Fig.3-12 所示。當 k=1 時，X₁(j) X(j) j 0,1,2...L1 當 k=2 時，X₂(j)X(jD) j0,1,2...L1 類推至 k=k 時，Xk(j)X(j(k1)D) j0,1,2...L1。將這 k 組訊號分別對其自相關函數做傅立葉變換，即可得到 k 個週期圖。將這 k 個週期圖做平均之後的結果即為該系統的功率譜密度[5]。

3.4.1 輸入噪音的功率譜

先前說過隨機共振不可或缺的三個因素是非線性系統、弱外場以及噪音，我們假設系統的雜訊是來自於分子的熱運動，因此所考慮的噪音是高斯型的白噪音(t)，滿足 ) ' ( 2 ) (t' (t) 0 (t) t t D        〉〈〉〈 3-13 其中 D 為分布函數之變異數。完美的白噪音在各個時間都是獨立的亂數，在實際 上很難達成。在這裡我們所使用產生噪音的方式是用 Matlab 內建的指令”rand”產 生一組平均為零、呈高斯分布的亂數序列，將它帶入模擬的系統中，再用時間序列做為變數從上述的噪音序列隨意取出作為不同時刻的值做為當下的噪音強度。輸入噪音的強度與功率譜密度如 Fig.3-13 上圖所示。此高斯白噪音的功率譜密度為 D e d D e t t d w S iw iw      



       ) ( ) ( ) ( ) ( 〈〉 3-14 可看出高斯白噪音的功率譜密度在每一個頻率之下的強度都一樣，而由 Fig.3-13 下圖可看出由 Matlab 所產生的高斯白噪音其頻譜密度之強度在各個頻率之下幾乎都相等，因此這樣子的噪音序列可是為白噪音來進行處理。 Fig. 3-12 Welch 對資訊點之切割示意圖。

(42)

0 1 2 3 4 5 6 7 8 9 10 x 104 0 0.5 1 time series(sec) n o ise a m p li tu d e (d B ) 10-4 10-3 10-2 10-1 100 -40 -20 0 20 40

Welch Power Spectral Density Estimate

Normalized Frequency ( rad/sample)

P o w e r/ fr e q u e n cy ( d B /r a d /sa m p le ) Fig. 3-13 由 Matlab 所產生的高斯白噪音與其功率譜。左圖為一段時間序列之下噪音的振幅，右圖是利用 Welch 方法所得到的功率譜密度，注意 x 軸是 Log 的尺度。

3.4.2 輸出訊號的功率譜密度

系統模擬所得到的輸出訊號是在弱外場的作用下每個時刻系統分別處在這兩個狀態(狀態的劃分如圖 3-4)的個數。在相同的隨時間變化之微弱外場下，改變噪音的強度會造成每個時刻處於兩個狀態的個數產生改變，如 Fig.3-14 所示。在無噪音的情形，由於外場的強度太弱無法將粒子由一個狀態踢到另一個狀態，即外場所施予的力無法使粒子越過位能井高度，因此這些一開始處於平衡態的一百個不同的初始狀態會平均的分布這兩個位能井區域，呈現兩個位能井區域的系統數為(50,50)之分布；當噪音小幅度的增加時，噪音給予系統的粒子一個隨機的力量，可能會使得原來無法被外場驅動的粒子能夠越過位能的高點而躍遷至另一個位能井中，因此系統處與某個狀態的粒子數會有跟著外場變化的趨勢。此時外場若再繼續的加大，會變成這個施予粒子的隨機力太大而造成多次的跳動，這個太過頻繁而隨機的跳躍使得我們所看到的個數變化變得雜亂。

二維磁偶極系統的隨機共振

國立交通大學

物理研究所

碩士論文