陳云龍 伍歆

(南昌大學(xué)理學(xué)院,南昌 330031)

(2013年1月11日收到;2013年3月27日收到修改稿)

1 引言

辛積分器[1]和流形改正方法[2,3]都屬于幾何積分算法[4],兩者同樣都能確保能量積分誤差在長時(shí)間積分中沒有長期變化,但差別在于前者能夠保持Hamilton系統(tǒng)的辛結(jié)構(gòu),而后者卻不能.辛方法因具有這兩個(gè)優(yōu)點(diǎn),從而成為求解Hamilton系統(tǒng)長期演化的最佳數(shù)值積分方法,它已在分子動力學(xué)[5]、量子力學(xué)[6-8]、聲波[9]、等離子體[10]、太陽系動力學(xué)和相對論天體物理[11-14]等領(lǐng)域得到廣泛應(yīng)用.

辛積分器分為顯式[15]和隱式[1]兩種形式.有時(shí)也會用到它們的混合形式[11-14,16-18],但這種混合型本質(zhì)上還是歸結(jié)于隱式方法.顯辛方法比隱辛方法具有計(jì)算效率方面的優(yōu)勢,故通常受到學(xué)者們的關(guān)注.應(yīng)當(dāng)注意它的使用要求Hamilton至少要分解成兩個(gè)可積部分,例如由動量和坐標(biāo)分別組成的動能與勢能就是一種Hamilton分解形式.二階Verlet辛積分器[4]盡管是這種分解形式下的一種簡單的顯辛方法,但在辛方法研究中卻發(fā)揮重要作用,被認(rèn)為是獲得高階或高精度辛積分算法的基礎(chǔ).直接消去其三階截?cái)嗾`差項(xiàng),若按不對稱組合方式可以得到三階辛方法[15];若按對稱組合方式則可以取得四階Forest-Ruth辛積分器[19],還能得到Y(jié)oshida高階辛積分器[20].當(dāng)Hamilton分解成主要的未受攝部分和次要的攝動部分并且兩部分均可積時(shí)[21],Verlet辛積分器的三階截?cái)嗾`差中有一項(xiàng)會比另一項(xiàng)大很多.若將該階甚至各階截?cái)嗾`差中的主要項(xiàng)消去但次要項(xiàng)仍保留,則得到類高階(pseudo-high-order)辛方法[22,23],還可以取得辛校正器[24].在太陽系動力學(xué)中,這種攝動分解形式及構(gòu)造的有關(guān)算法與動-勢能分解的同階算法來比顯著提高了數(shù)值精度.應(yīng)該指出攝動分解情形下Verlet辛積分器的三階截?cái)嗾`差中的所謂次要項(xiàng)當(dāng)實(shí)施動-勢能分解時(shí)與另一項(xiàng)不再有明顯的主次之分,亦即沒有很大的量級差異,將會在算法構(gòu)造中很有用,它本質(zhì)上就是力的梯度.Ruth[15]正是認(rèn)識到這一點(diǎn),把它加入勢能中再與動能組合求解構(gòu)造了一個(gè)三階力梯度辛方法.沿著該思路,可以設(shè)計(jì)一系列高階力梯度辛方法[25-30]及其截?cái)嗾`差各項(xiàng)系數(shù)平方和最小的最優(yōu)化型算法[31,32].不超過四階的力梯度辛方法各組合系數(shù)可以全為正數(shù),意味著積分過程中的每一子步都用正步長,這有益于求解不可逆系統(tǒng).大量文獻(xiàn)[25-32]表明力梯度辛方法在Kepler二體問題、攝動Kepler二體問題、太陽系N體問題、分子動力學(xué)多粒子問題、量子力學(xué)薛定諤方程以及H′enon-Heiles系統(tǒng)等模擬中比同階非力梯度辛方法要好一到幾個(gè)數(shù)量級精度,并且最優(yōu)化型總要優(yōu)于沒有經(jīng)過優(yōu)化處理的算法.類高階辛方法、辛校正器和力梯度辛方法實(shí)際上都是直接將誤差項(xiàng)中的某些部分提取出來加入組合來構(gòu)造新的算法,以便使所得算法誤差里不再含有該提取部分,從而改善數(shù)值精度.從這一角度來看,它們均含有誤差校正的意思.

考慮到四階力梯度辛方法與非力梯度Forest-Ruth辛積分器來比具有精度好的優(yōu)點(diǎn),文獻(xiàn)[26]將其應(yīng)用于求解慣性坐標(biāo)系下的圓型限制性三體問題.這里第三個(gè)小天體受到兩主天體所給引力的梯度連同引力被一起用來構(gòu)造算法.這種處理是合理的,因?yàn)閼T性系下的第三體受力在物理上是非常清晰的.然而,圓型限制性三體問題的研究大都在旋轉(zhuǎn)坐標(biāo)系下進(jìn)行,這時(shí)第三體除受到兩主天體引力外還受到旋轉(zhuǎn)坐標(biāo)系所產(chǎn)生的非慣性力,使得第三體受力變得復(fù)雜起來,進(jìn)而導(dǎo)致力梯度算法的應(yīng)用變得困難,因?yàn)槲覀儾磺宄μ荻染烤故悄姆N力的梯度,即是關(guān)于引力的梯度還是非慣性力的?也許有人認(rèn)為應(yīng)當(dāng)是針對引力與非慣性力的合力來實(shí)施梯度運(yùn)算,但這樣會使引力與非慣性力對應(yīng)的Hamilton之和不可積,也就是顯辛方法不再適合,故顯辛方法在實(shí)際應(yīng)用中需要將引力與非慣性力分開計(jì)算.因此,我們不直接從第三體受力分析出發(fā)而是改由有關(guān)Lie算子運(yùn)算來探討力梯度辛方法是否適合求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題.這是本文的一個(gè)主要目的.具體來說,通過評估比較四階力梯度辛方法[25]、最優(yōu)化四階力梯度辛方法[31]和Forest-Ruth辛方法[19],我們希望找到適合求解旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題精度最好的辛積分器.然后,利用所挑出的算法計(jì)算兩鄰近軌道的Lyapunov指數(shù)[33]和快速Lyapunov指標(biāo)[34],達(dá)到完全確保高精度辛方法能夠貫穿于這些混沌指標(biāo)計(jì)算的全過程,可為正確揭示系統(tǒng)動力學(xué)定性性質(zhì)服務(wù).

2 算法描述

設(shè)動能T(p)是關(guān)于n維動量p的二次型T(p)=p2/2,而勢能V(q)僅是n維坐標(biāo)q的函數(shù),它們組成的Hamilton函數(shù)為

依次定義T和V的Lie導(dǎo)數(shù)算子為

上式中,Ti=?/?pi,Vi=?/?qi,以下類同.Lie導(dǎo)數(shù)A分別作用于坐標(biāo)qi和動量pi,實(shí)際上就是關(guān)于坐標(biāo)和動量的普通導(dǎo)數(shù),即q˙i=Aqi=Ti,p˙i=Api≡0,表明A僅僅對位置坐標(biāo)有作用,故稱為位置型算子;B自然稱為動量型算子,它對pi作用就得到力或加速度.利用這兩個(gè)算子可以構(gòu)造二階Verlet辛積分器[4]

這里,τ為步長,而

在三階截?cái)嗾`差中,互逆子

容易驗(yàn)證

但

于是,可得

其中,力 f=(-V1,···,-Vn).顯然,C 就是力梯度算子,與B屬于同類的動量型算子.

既然如此,可將截?cái)嗾`差中的C提取出來并和B等同看待再跟A組合構(gòu)造一類更高階算法,稱為力梯度辛方法.例如,一個(gè)簡單四階力梯度辛方法[25]為

如還要求五階截?cái)嗾`差各項(xiàng)系數(shù)的平方和最小,則得一個(gè)最優(yōu)化四階力梯度辛方法[31]

其中,各步進(jìn)算子的系數(shù)約定如下：

當(dāng)不考慮C加入組合式中,則是通常的辛方法.常見的例子就是四階Forest-Ruth辛積分器[19]

3 圓型限制性三體問題

本節(jié)致力于力梯度型算法在旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題中的具體實(shí)現(xiàn),并以非力梯度型為參考評估這些算法的精度優(yōu)劣,最后利用所挑選的精度最佳算法計(jì)算一些混沌指標(biāo).

3.1 物理模型

經(jīng)典平面圓型限制性三體問題[35]揭示的是第三個(gè)無限小質(zhì)量天體在做勻速圓周運(yùn)動的兩主天體的引力作用下的運(yùn)動情形.約定兩主天體總質(zhì)量為1,質(zhì)量為m1=1-μ的一體位于質(zhì)心旋轉(zhuǎn)坐標(biāo)系內(nèi)的x軸上的點(diǎn)(-μ,0),而質(zhì)量為m2=μ≤1/2的另一體位于點(diǎn)(1-μ,0).這樣兩主天體固聯(lián)在x軸上始終保持1不變,但它們以勻角速度1做圓運(yùn)動,即x軸圍繞原點(diǎn)以勻角速度1旋轉(zhuǎn)并不受小天體引力的影響.在這兩主天體的引力作用下小天體的坐標(biāo)q=(x,y)和動量p=(px,py)演化以無量綱Hamilton函數(shù)表為

力函數(shù)U(x,y)具有如下形式

該系統(tǒng)具有正則方程

容易驗(yàn)證存在Jacobi常數(shù)

顯辛方法的應(yīng)用并不針對整個(gè)系統(tǒng)(10)或(11)而是將其分離成兩個(gè)可積部分分別求解.系統(tǒng)(10)的一種分解形式可以是如下的“動能T”和“勢能V”兩部分

勢能V來自于兩主天體施加給小天體的引力作用部分,而ypx-xpy是因旋轉(zhuǎn)坐標(biāo)系為非慣性系而給小天體運(yùn)動附加的影響部分.

從(13)式來看,T不是動量p的二次型,這使得力梯度算法的應(yīng)用似乎遇到了困難,下面力圖解決這一問題.

3.2 算法實(shí)現(xiàn)

由于T不是動量p的嚴(yán)格二次型,故其Lie導(dǎo)數(shù)算子A不再有方程(2)形式,而是由位置型算子變?yōu)閯恿颗c位置混合型算子

它作用于坐標(biāo)和動量就得到微分方程組

從初始狀態(tài)量(x0,y0,px0,py0)出發(fā)經(jīng)過時(shí)間t后獲得該方程組的分析解為這意味著混合算子A是可積可解的.關(guān)于V的動量型算子B也是如此,但若非慣性系所施加的影響部分ypx-xpy與V結(jié)合便是不可積、不可解的.

當(dāng)(14)式所給算子A代替(2)式中的算子A,我們?nèi)菀鬃C明(4)和(5)兩式仍然成立,從而(6)式也仍然適用,即動量型算子

這表明算子C還是兩主天體引力的梯度,并不是兩主天體引力與旋轉(zhuǎn)坐標(biāo)系所附加非慣性力的合力之梯度.因此,旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題仍可如(1)式所示Hamilton系統(tǒng)那樣用力梯度方法求解.為了更清楚地說明力梯度辛方法的具體實(shí)現(xiàn)過程,我們以算法S1為例列舉其從l步到(l+1)步的差分格式

我們采用Jacobi常數(shù)CJ=3.12及大約相當(dāng)于木星與太陽質(zhì)量比參數(shù)μ=0.001和初值x=0.29,y=px=0,而初值py從能量積分(10)式或Jacobi積分(12)式確定并在開方運(yùn)算中取正根.該積分軌道是圖1(a)所示Poincare′截面y=0(y˙＞0)上的一個(gè)環(huán),表明是擬周期有序軌道.讓步長τ從0.01到0.1逐次變大,可得圖1(b)所示的每個(gè)算法對每個(gè)給定步長積分100000步后的最大Jacobi常數(shù)誤差.該圖揭示了非力梯度Forest-Ruth辛方法S3的精度要顯著劣于力梯度算法S1精度.還可以看出最優(yōu)化力梯度算法S2是最好的.

圖1 (a)Jacobi常數(shù)C J=3.12的Poincar′e截面y=0上的擬周期有序軌道;(b)三個(gè)算法對每個(gè)給定步長積分100000步后的最大Jacobi常數(shù)誤差

圖2 (a)Jacobi常數(shù)C J=3.06的Poincar′e截面y=0上的混沌軌道;(b)三個(gè)算法對每個(gè)給定步長積分100000步后的最大Jacobi常數(shù)誤差

當(dāng)僅改變Jacobi常數(shù)CJ=3.06及相應(yīng)初值py時(shí),則得圖2(a)所示Poincar′e截面上的一些隨機(jī)分布的點(diǎn),表明是混沌軌道.圖2(b)的最大Jacobi常數(shù)誤差再次證實(shí)三個(gè)算法的數(shù)值性能優(yōu)劣關(guān)系與有序軌道的情形是完全一致的.就圖1(b)和2(b)比較來看,當(dāng)步長τ沒有達(dá)到0.1前,對于有序和混沌軌道兩情形每個(gè)算法精度基本上差不多;但當(dāng)步長τ接近0.1時(shí),針對混沌軌道情況,兩算法S1和S3已失去數(shù)值穩(wěn)定性,而算法S2仍幾乎保持有序軌道的精度.

因此,無論是對有序軌道還是對混沌軌道進(jìn)行數(shù)值模擬,我們總是發(fā)現(xiàn)力梯度類算法精度要大大優(yōu)于非力梯度類的,而最優(yōu)化型算法能夠取得最好精度.這樣,高精度辛方法的使用能確保所得動力學(xué)定性結(jié)果的可靠.下面我們考慮方法S2在混沌動力學(xué)方面的應(yīng)用.

3.3 混沌指標(biāo)

除Poincar′e截面方法外,Lyapunov指數(shù)和快速Lyapunov指標(biāo)等都是區(qū)分系統(tǒng)有序與混沌的常用方法.這些方法克服了Poincare′截面限于相空間維數(shù)為3的應(yīng)用局限,并且均以計(jì)算兩鄰近軌道的分離演化為基礎(chǔ)來獲得.

Lyapunov指數(shù)是衡量兩鄰近軌道平均指數(shù)分離比的指標(biāo).通常采用的是最大Lyapunov指數(shù),定義為

式中的ξ(t)和ξ(0)分別是t時(shí)刻與初始時(shí)刻的切向量,即變分方程之解.(18)式稱為計(jì)算Lyapunov指數(shù)的變分法[33].當(dāng)然,它的計(jì)算也可用兩鄰近軌道的差來代替變分方程之解,稱為計(jì)算Lyapunov指數(shù)的兩粒子法[33],定義如下：

d(t)和d(0)是兩鄰近軌道分別在t時(shí)刻與初始時(shí)刻的距離.應(yīng)當(dāng)注意初始距離必須適當(dāng)選擇,如果d(0)太大,那么兩鄰近軌道的差將與變分方程的解差別太大;但若d(0)太小,則舍入誤差增大.也就是說,d(0)太大和太小都會導(dǎo)致不正確的Lyapunov指數(shù).文獻(xiàn)[36]認(rèn)為就雙精度而言初始距離的適當(dāng)選擇為d(0)=10-8,還指出應(yīng)適當(dāng)考慮重整化次數(shù).對于有界軌道,若λ＞0,則系統(tǒng)混沌;當(dāng)λ=0時(shí),系統(tǒng)被認(rèn)為是有序的.這是Lyapunov指數(shù)區(qū)分有序和混沌的準(zhǔn)則.

注意(19)式代替(18)式計(jì)算Lyapunov指數(shù)的優(yōu)點(diǎn)有二：對于復(fù)雜的動力系統(tǒng),如相對論引力系統(tǒng),可以避免推導(dǎo)復(fù)雜變分方程的麻煩;針對Hamilton系統(tǒng)來說,可以完全確保辛方法貫穿于Lyapunov指數(shù)計(jì)算的全過程,因?yàn)閮闪Ｗ臃ň褪怯眯练椒ǚe分運(yùn)動方程兩次來計(jì)算Lyapunov指數(shù).基于此因,我們采用兩粒子法計(jì)算Lyapunov指數(shù).選取步長τ=0.01,參數(shù)和參考軌道初始條件同上,而鄰近軌道初始條件為x=0.29+10-8,其他初始條件與參考軌道的相同但初值py相應(yīng)改變.每積分10步,實(shí)施一次重整化.于是,我們獲得圖3中的Lyapunov指數(shù).圖1軌道的Lyapunov指數(shù)λ趨于0故為有序的,而圖2軌道的λ=0.023則表明其混沌性.

一個(gè)比Lyapunov指數(shù)區(qū)分有序與混沌更靈敏、更快的方法是快速Lyapunov指標(biāo),其定義如下[37]：

也可用兩鄰近軌道的差來代替變分方程之解來計(jì)算快速Lyapunov指標(biāo),稱為兩鄰近軌道(或兩粒子法)的快速Lyapunov指標(biāo)[34],定義為

不同于(20)式,(21)式計(jì)算快速Lyapunov指標(biāo)需要重整化,否則,混沌情形的兩鄰近軌道距離d(t)達(dá)到1后出現(xiàn)軌道飽和使其不再增長、計(jì)算無法繼續(xù).為避免這一問題,每當(dāng)d(t)=1實(shí)施重整化,即不改變方向?qū)⑧徑壍览氐脚c參考軌道相距d(0)處進(jìn)行它的下次積分.假設(shè)重整化次數(shù)為k,那么,兩粒子法的快速Lyapunov指標(biāo)實(shí)際按下式計(jì)算

對于有序和混沌兩種情況,切向量長度隨時(shí)間演化的速度是完全不同的,即有序系統(tǒng)的切向量長度隨時(shí)間代數(shù)式變大,而混沌系統(tǒng)的切向量長度隨時(shí)間指數(shù)式增長.這一特性可以作為有序與混沌的區(qū)分標(biāo)準(zhǔn).

由于兩鄰近軌道的快速Lyapunov指標(biāo)具有兩粒子法計(jì)算Lyapunov指數(shù)那樣的優(yōu)點(diǎn),我們考慮它的應(yīng)用.選擇鄰近軌道初始條件為x=0.29+10-9,其他初始條件和參考軌道的與上面相同但初值py相應(yīng)改變.實(shí)際上,初始距離為d(0)≈10-9,可以比兩粒子法的Lyapunov指數(shù)的初始距離要小些,而且重整化時(shí)間間隔要大很多.如圖4所示,圖1軌道的快速Lyapunov指標(biāo)隨時(shí)間增長緩慢,屬于有序軌道類型;圖2軌道的快速Lyapunov指標(biāo)隨時(shí)間增長迅速,應(yīng)歸之于混沌軌道一類.

對比圖3和圖4還可以看出利用Lyapunov指數(shù)區(qū)分有序與混沌需要積分的時(shí)間至少為t=10000,而根據(jù)快速Lyapunov指標(biāo)只需積分到t=1000即可分辨兩軌道類型.快速Lyapunov指標(biāo)由于擁有識別混沌速度快的優(yōu)點(diǎn),因而常用來追蹤小行星軌道穩(wěn)定性和研究其他動力學(xué)問題[38-40].

圖3 兩軌道的Lyapunov指數(shù)

圖4 兩軌道的快速Lyapunov指標(biāo)

4 結(jié)論

旋轉(zhuǎn)坐標(biāo)系下的圓型限制性三體問題因含非慣性系所附加的影響部分使得動能不是動量的嚴(yán)格二次型,意味著力梯度辛積分算法的應(yīng)用似乎遇到了理論方面的困難.本文從Lie算子運(yùn)算出發(fā)嚴(yán)格推導(dǎo)了力梯度算子的物理意義仍然是兩主天體對第三個(gè)小天體引力的梯度,并不是兩主天體引力與旋轉(zhuǎn)坐標(biāo)系所附加非慣性力的合力之梯度.這樣,解決了力梯度辛方法應(yīng)用的理論問題.數(shù)值實(shí)驗(yàn)揭示力梯度類算法精度總要大大優(yōu)于非力梯度類的,而最優(yōu)化型算法能夠取得最好精度.將最優(yōu)化型算法計(jì)算兩鄰近軌道的Lyapunov指數(shù)和快速Lyapunov指標(biāo)可以完全確保高精度辛方法能夠貫穿于這些混沌指標(biāo)計(jì)算的全過程,以便準(zhǔn)確刻畫動力學(xué)定性性質(zhì).

未來的工作希望借助這種最優(yōu)化型辛積分算法和恰當(dāng)?shù)幕煦缰笜?biāo)探討含輻射壓的平面圓型限制性三體問題、平面橢圓型限制性三體問題、空間圓型限制性三體問題和空間橢圓型限制性三體問題等動力學(xué)性質(zhì).

[1]Feng K,Qin M Z 2009 Symplectic Geometric Algorithmsfor Hamiltonian Systems(Hangzhou：Zhejiang Scienceand Technology Publishing House)

[2]Zhong SY,Wu X 2010 Phys.Rev.D 81 104037

[3]Mei L J,Wu X,Liu FY 2012 Chin.Phys.Lett.29 050201

[4]Hairer E,Lubich C,Wanner G 1999 Geometric Numerical Integration.(Berlin：Springer)

[5]Chi Y H,Liu X S,Ding PZ 2006 Acta Phys.Sin.55 6320(in Chinese)[匙玉華,劉學(xué)深,丁培柱2006物理學(xué)報(bào)55 6320]

[6]Luo X Y,Liu X S,Ding PZ 2007 Acta Phys.Sin.56 604(in Chinese)[羅香怡,劉學(xué)深,丁培柱2007物理學(xué)報(bào)56 604]

[7]Liu X S,Wei JY,Ding PZ 2005 Chin.Phys.14 231

[8]Bian X B,Qiao H X,Shi T Y 2007 Chin.Phys.16 1822

[9]Cao Y,Yang K Q 2003 Acta Phys.Sin.52 1984(in Chinese)[曹禹,楊孔慶2003物理學(xué)報(bào)52 1984]

[10]Hu WP,Deng Z C 2008 Chin.Phys.B 17 3923

[11]Zhong SY,Wu X,Liu SQ,Deng X F 2010 Phys.Rev.D 82 124040

[12]Zhong SY,Wu X 2011 Acta Phys.Sin.60 090402(in Chinese)[鐘雙英,伍歆2011物理學(xué)報(bào)60 090402]

[13]Zhong SY,Liu S 2012 Acta Phys.Sin.61 120401(in Chinese)[鐘雙英,劉崧2012物理學(xué)報(bào)61 120401]

[14]Wu X,Zhong SY 2011 Gen.Relat.Gravit.43 2185

[15]Ruth RD 1983 IEEETran.Nucl.Sci.30 2669

[16]Preto M,Saha P 2009 Astrophy.J.703 1743

[17]Liao X H 1997 Celest.Mech.Dyn.Astron.66 243

[18]Lubich C,Walther B,Braugmann B 2010 Phys.Rev.D 81 104025

[19]Forest E,Ruth RD 1990 Physica D 43 105

[20]Yoshida H 1990 Phys.Lett.A 150 262

[21]Wisdom J,Holman M 1991 Astron.J.102 1528

[22]Preto M,Tremaine S 1999 Astron.J.118 2532

[23]Laskar J,Robutel P 2001 Celest.Mech.Dyn.Astron.80 39

[24]Wisdom J,Holman M,Touma J1996 Fields Inst.Commun.10 217

[25]Chin SA 1997 Phys.Lett.A 75 226

[26]Chin SA,Chen CR 2005 Celest.Mech.Dyn.Astron.91 301

[27]Chin SA 2007 Phys.Rev.E 75 036701

[28]Xu J,Wu X 2010 Res.Astron.Astrophys.10 173

[29]Sun W,Wu X,Huang G Q 2011 Res.Astron.Astrophys.11 353

[30]Li R,Wu X 2010 Science China：Physics,Mechanics&Astronomy 53 1600

[31]Omelyan IP,Mryglod IM,Folk R 2003 Comput.Phys.Commun.151 272

[32]Li R,Wu X 2010 Acta Phys.Sin.59 7135(in Chinese)[李榮,伍歆2010物理學(xué)報(bào)59 7135]

[33]Wu X,Huang T Y 2003 Phys.Lett.A 313 77

[34]Wu X,Huang T Y,Zhang H 2006 Phys.Rev.D 74 083001

[35]Murray C D,Dermott SF 1999 Solar System Dynamics(Cambridge,UK：Cambridge Univ.Press)

[36]Tancredi G,S′anchez A,Roig F 2001 Astron.J.121 1171

[37]Froeschl′e C,Lega E 2000 Celest.Mech.Dyn.Astron.78 167

[38]Wu X,Xie Y 2008 Phys.Rev.D 77 103012

[39]Wang Y,Wu X 2012 Chin.Phys.B 21 050504

[40]Wang Y Z,Wu X,Zhong S Y 2012 Acta Phys.Sin.61 160401(in Chinese)[王玉詔,伍歆,鐘雙英2012物理學(xué)報(bào)61 160401]