顧建宇，董小波

(1.中國(guó)科學(xué)院 云南天文臺(tái)，昆明650216；2.中國(guó)科學(xué)院大學(xué)，北京100049)

1 引言和理論基礎(chǔ)

銀河系研究的主要目標(biāo)之一，是測(cè)量并計(jì)算出銀河系的結(jié)構(gòu)，即各種組成成分(如恒星核球、恒星盤(pán)、暗物質(zhì)暈和氣體等)的空間分布，以及它們的速度空間信息。為了實(shí)現(xiàn)這個(gè)目標(biāo)，在觀測(cè)上，需要高精度的天體測(cè)量學(xué)及測(cè)光學(xué)數(shù)據(jù)(如Gaia)和光譜數(shù)據(jù)(如LAMOST)；在理論和計(jì)算上，星系動(dòng)力學(xué)建模是核心。星系動(dòng)力學(xué)建模的理論基礎(chǔ)是恒星系統(tǒng)的動(dòng)力學(xué)，即研究大量的點(diǎn)質(zhì)量粒子在它們共同的引力場(chǎng)中的運(yùn)動(dòng)[8]。關(guān)于銀河系建模，盡管星系作為自引力系統(tǒng)，沒(méi)有星系能處于嚴(yán)格的熱力學(xué)平衡態(tài)或完美的動(dòng)力學(xué)平衡態(tài)，然而平衡態(tài)動(dòng)力學(xué)模型(或稱(chēng)穩(wěn)態(tài)[steady state]，或準(zhǔn)定態(tài)[quasi-stationary state])，對(duì)于解釋銀河系和河外星系的觀測(cè)至關(guān)重要。從實(shí)際情況看，暗物質(zhì)暈作為星系引力勢(shì)的主要提供者，僅在短暫且罕見(jiàn)的主并合事件中，才遠(yuǎn)離動(dòng)態(tài)平衡，一般情況下該動(dòng)態(tài)平衡僅受吸積或次并合等事件的小擾動(dòng)。從模型需求的角度看，由動(dòng)力學(xué)測(cè)量數(shù)據(jù)推斷出物質(zhì)分布，該模型要求假設(shè)星系處于(準(zhǔn))平衡態(tài)；另外，平衡態(tài)模型是最簡(jiǎn)單的模型，更復(fù)雜的情況可以被建模為平衡態(tài)模型的“擾動(dòng)”、“變形”或迭加[20]。

運(yùn)動(dòng)積分是粒子在軌道運(yùn)動(dòng)中的不變量；作用量，則是可以與角變量形成正則坐標(biāo)的一類(lèi)特殊的運(yùn)動(dòng)積分。由于在“作用量-角變量”坐標(biāo)下運(yùn)動(dòng)方程形式的簡(jiǎn)單性，作用量成為很有吸引力的變量。并且作用量無(wú)論在微擾論方法上還是分布函數(shù)的建模上，都是十分必要的。但是，粒子在一般引力勢(shì)(例如真實(shí)的星系)中的運(yùn)動(dòng)，基本上不能得到解析的運(yùn)動(dòng)積分或作用量；僅當(dāng)引力勢(shì)為可分離勢(shì)時(shí)，能解析地得到。Stckel勢(shì)是橢球坐標(biāo)系下可分離的引力勢(shì)形式，該引力勢(shì)下的運(yùn)動(dòng)積分可被解析地計(jì)算[1]。在應(yīng)用上，通過(guò)將一般引力勢(shì)近似或擬合為Stckel勢(shì)，人們可估算一般引力勢(shì)下的運(yùn)動(dòng)積分。近十幾年來(lái)，一些研究者基于Stckel引力勢(shì)理論中的解析式，用(參見(jiàn)3.2節(jié))解析或數(shù)值的(例如Binney的工作)方法來(lái)作拓展，即基于Stckel勢(shì)來(lái)估算一般星系中的作用量(或運(yùn)動(dòng)積分)。特別是最近幾年來(lái)，Binney和他的合作者們把他們的方法都實(shí)現(xiàn)為計(jì)算機(jī)程序，并以公開(kāi)源碼的方式釋放了(詳見(jiàn)第2章)。

接下來(lái)，第1章的其他小節(jié)分別介紹：恒星系統(tǒng)的動(dòng)力學(xué)及統(tǒng)計(jì)力學(xué)的概要、運(yùn)動(dòng)積分以及Stckel勢(shì)的相關(guān)基本背景知識(shí)[1]。第2章介紹估算一般引力勢(shì)的運(yùn)動(dòng)積分和作用量的方法：2.1節(jié)介紹Stckel Fudge作弊法[13,20,22,24]；2.2節(jié)介紹Stckel引力勢(shì)擬合法[19]，以及更早期的工作；2.3節(jié)介紹一類(lèi)環(huán)面坐標(biāo)變換的作用量估算方法，基于Stckel玩具勢(shì)模型得到一般引力勢(shì)下的作用量[20,21,24]；2.4節(jié)是我們基于文獻(xiàn)[24]所介紹的TACT程序?qū)Ω鞣N方法做的比較；2.5節(jié)則是介紹一些其他求運(yùn)動(dòng)積分的方法和程序。第3章介紹近幾年來(lái)這些方法在銀河系數(shù)據(jù)上的應(yīng)用(如文獻(xiàn)[44])。第4章進(jìn)行討論和總結(jié)，特別是簡(jiǎn)介“作用量-角變量”方法在不可積的星系現(xiàn)象(見(jiàn)4.2節(jié))和非平衡的星系過(guò)程(見(jiàn)4.3節(jié))中應(yīng)用的最新進(jìn)展。

1.1 恒星系統(tǒng)的理論基礎(chǔ)

本文討論的單粒子在整體勢(shì)的運(yùn)動(dòng)，是從引力多體運(yùn)動(dòng)簡(jiǎn)化而來(lái)的，即單粒子近似或稱(chēng)平均場(chǎng)近似[60]。三維背景空間中一個(gè)包含有N個(gè)粒子的經(jīng)典系統(tǒng)，可用哈密頓方程描述：

對(duì)于該系統(tǒng)所有粒子的坐標(biāo)和動(dòng)量所構(gòu)成的全6N維相空間(一般稱(chēng)為Γ相空間)，通過(guò)相空間代表點(diǎn)個(gè)數(shù)的局部守恒可推導(dǎo)得，其概率密度ρ(q,p)滿足劉維爾方程：

其中，ρ=ρ(q i,p i,t),i=1,2,3,···,3N是系綜意義的概率分布；ρ,H的泊松括號(hào)式(2)中，左邊其實(shí)就是拉格朗日視角的隨流導(dǎo)數(shù)dρ/dt，即其相體積不變；這就是劉維爾保體積定理。

下面用x=x j,p=p j(j=1,2,3)代表單個(gè)粒子的坐標(biāo)和動(dòng)量；q k,p k代表第k個(gè)粒子的坐標(biāo)和速度；f n代表n粒子聯(lián)合分布函數(shù)。對(duì)上述劉維爾方程作單粒子近似f2(q1,p1,q2,p2,t)=f1(q1,p1,t)f1(q2,p2,t)①這是一種平均場(chǎng)近似，即將粒子間相互作用當(dāng)成總背景場(chǎng)，而不考慮粒子間相互作用和關(guān)聯(lián)；體現(xiàn)在分布函數(shù)上就是每個(gè)粒子的單粒子分布是無(wú)關(guān)聯(lián)的，因此乘積為它們的聯(lián)合分布。，可得到玻爾茲曼方程②由玻爾茲曼于1872年提出。從劉維爾到玻爾茲曼方程的約化，中間還有一個(gè)BBGKY級(jí)聯(lián)方程。：

其中，這里f=f(q,p,t)是單粒子分布函數(shù)，即處于自由度3的坐標(biāo)和速度所構(gòu)成的6維相空間(一般稱(chēng)為μ相空間)；v代表速度，F(xiàn)代表受力③該受力是對(duì)于該“單粒子”而言的，施力者為“其他粒子”(實(shí)際上是分布函數(shù))的作用(因而在這里直接體現(xiàn)平均場(chǎng)近似的思想)，以及可能存在的外場(chǎng)。玻爾茲曼方程可以用來(lái)描述非自洽的體系(自洽體系是指粒子所受力全部來(lái)自于研究系統(tǒng)本身，如自引力系統(tǒng))，即通過(guò)建立多種組分的玻爾茲曼方程來(lái)描述。；σ(?)是分子微分散射截面，這里考慮了二體碰撞；對(duì)于其他作用類(lèi)型(如剛性)的散射，將可能不再是這樣的形式。玻爾茲曼方程體系中所需要的假設(shè)非常繁雜(甚至冗余)，而且至今仍有難以嚴(yán)格數(shù)學(xué)證明的規(guī)律。本文對(duì)玻爾茲曼方程僅作介紹而不作分析(更多參見(jiàn)文獻(xiàn)[8]中的第4,7章和文獻(xiàn)[60]中的第3,4章)。此外可以驗(yàn)證，劉維爾方程在Γ相空間的分布函數(shù)的熵是不變的，而玻爾茲曼方程在μ相空間的分布函數(shù)的熵是不減的。

對(duì)玻爾茲曼方程進(jìn)行局域平衡假設(shè)，碰撞項(xiàng)將消失。其實(shí)，在星系這種具體的系統(tǒng)中，恒星的半徑與星系尺度相比很小，其平均自由程遠(yuǎn)大于系統(tǒng)尺度，碰撞時(shí)間遠(yuǎn)大于穿越整個(gè)星系的時(shí)間，比宇宙年齡還大很多量級(jí)，因此碰撞項(xiàng)本應(yīng)忽略不計(jì)。把玻爾茲曼方程兩邊乘上一個(gè)碰撞不變量ξ(x,v,t)(如質(zhì)量、動(dòng)量、能量等)，然后作局域平均(由于碰撞前后的某些對(duì)稱(chēng)性,局域平均后碰撞項(xiàng)不做貢獻(xiàn)，見(jiàn)文獻(xiàn)[60]第7A章)，可分別得到流體力學(xué)性質(zhì)的方程(注：流體力學(xué)過(guò)程是局域平衡過(guò)程)。例如對(duì)于動(dòng)量守恒，玻爾茲曼方程兩邊乘以ξ(x,v,t)=mv i(i=x,y,z)，通過(guò)對(duì)動(dòng)量空間積分并利用動(dòng)量空間的散度定理(無(wú)窮遠(yuǎn)處分布函數(shù)為0，limp→∞f=0)，可得到

此式實(shí)質(zhì)上等同于流體力學(xué)動(dòng)量守恒方程：

這樣，我們就有了“N體系統(tǒng)(力學(xué)描述)-劉維爾方程(6N維)-玻爾茲曼方程(6維)-流體力學(xué)方程(流場(chǎng))”的聯(lián)系。

另一方面，從無(wú)碰撞玻爾茲曼方程出發(fā)，兩邊乘以速度分量，再對(duì)速度空間積分并利用散度定理(與得到式(4)的方法類(lèi)似)，我們得到Jeans方程：

其中，ν代表位形空間數(shù)密度，σij為速度彌散張量。Jeans方程不僅通過(guò)平均簡(jiǎn)化了方程，還將分布函數(shù)與這些可觀測(cè)量聯(lián)系起來(lái)，從而對(duì)聯(lián)系觀測(cè)天文的工作很有幫助。Jeans方程還給人們提供了一種檢測(cè)分布函數(shù)正確性的方法。當(dāng)然,Jeans方程仍然不是閉合的。

在μ相空間中的粒子分布函數(shù)f=f(x,v,t)①由于p=mv(m為單粒子的質(zhì)量)，而且在星系動(dòng)力學(xué)的粗糙建模中常常將測(cè)試天體考慮為無(wú)窮多同質(zhì)粒子的分布，且為方便一般假設(shè)測(cè)試粒子的質(zhì)量為1，于是下面我們這里的習(xí)慣，以下討論的f都是關(guān)于x,v的分布函數(shù)。，跟粒子視角的單粒子近似下的運(yùn)動(dòng)(x(t),v(t))有著密切的聯(lián)系。這種單粒子分布函數(shù)，實(shí)際上假設(shè)體系是由無(wú)窮多的點(diǎn)粒子構(gòu)成，某種程度上與場(chǎng)論的意義比較像(場(chǎng)論就是處理無(wú)窮多自由度)。從場(chǎng)的角度來(lái)理解，這是把所有粒子之間的相互作用當(dāng)做背景場(chǎng)，然后對(duì)于任一粒子，只須研究它在這個(gè)場(chǎng)中的運(yùn)動(dòng)，不必考慮其他粒子對(duì)它的作用；此即前述平均場(chǎng)的思想。但是，實(shí)際上想詳細(xì)了解星系結(jié)構(gòu)與演化(如分布函數(shù))，直接從自引力系統(tǒng)的統(tǒng)計(jì)物理方程出發(fā)去求解分布函數(shù)并不容易，而且運(yùn)用最大熵原理之類(lèi)的方法也是不夠的(由于星系作為長(zhǎng)程相互作用系統(tǒng)的特性，見(jiàn)第4章)。為得到分布函數(shù)，一種方法是N體模擬，基于模擬的結(jié)果作統(tǒng)計(jì)。另一種方法是根據(jù)粒子在引力勢(shì)的運(yùn)動(dòng)特性來(lái)構(gòu)建分布函數(shù)模型[8]，以及其他理論手段(詳見(jiàn)文獻(xiàn)[8])；在這種方法中，對(duì)運(yùn)動(dòng)積分的研究是個(gè)基礎(chǔ)的環(huán)節(jié)。下面，我們主要從哈密頓力學(xué)以及單粒子近似下的恒星軌道運(yùn)動(dòng)的角度出發(fā)，介紹恒星在靜態(tài)引力勢(shì)下的運(yùn)動(dòng)積分的相關(guān)內(nèi)容。

1.2 運(yùn)動(dòng)積分和作用量-角變量體系簡(jiǎn)介

1.2.1 哈密頓-雅克比方程

力學(xué)系統(tǒng)在相空間的正則坐標(biāo)(q,p)的演化，可以看做是無(wú)窮多個(gè)連續(xù)以哈密頓量為生成函數(shù)的無(wú)窮小正則變換的結(jié)果。在正則變換中令新的哈密頓量為0，并選取第二類(lèi)生成函數(shù)F2②第二類(lèi)生成函數(shù)是指：在正則變換中，取獨(dú)立的廣義坐標(biāo)變量為老的坐標(biāo)和新的動(dòng)量，亦即生成函數(shù)為老坐標(biāo)、新動(dòng)量以及時(shí)間的函數(shù)?？墒鼓撤N新的正則坐標(biāo)Q,P均為常數(shù)，即Q=α,P=β，可得正則變換的方程：

被稱(chēng)為哈密頓-雅克比方程(哈雅方程)，其解的形式為S=S(q,β)就是生成函數(shù)F2。另外，第二類(lèi)正則變換中哈密頓函數(shù)之間的微分關(guān)系表示為：

1.2.2 運(yùn)動(dòng)積分

如果一個(gè)力學(xué)量不顯含時(shí)間而且保持不變，h(q,p)=const，就稱(chēng)之為系統(tǒng)的運(yùn)動(dòng)積分，即運(yùn)動(dòng)積分是粒子在軌道運(yùn)動(dòng)中的不變量。尋找系統(tǒng)的這種不變量，具有重要意義。對(duì)于一個(gè)自治哈密頓系統(tǒng)，哈密頓量不顯含時(shí)間。此時(shí)哈密頓量本身就是一個(gè)運(yùn)動(dòng)積分，其積分常數(shù)E等于系統(tǒng)總能量。利用公式(7)和(8)，哈雅方程的解可寫(xiě)成時(shí)間項(xiàng)和廣義坐標(biāo)項(xiàng)分離的形式，表示為：

式中的W稱(chēng)為哈密頓特征函數(shù)。對(duì)于自由度n的自洽哈密頓系統(tǒng)，最多有n個(gè)兩兩獨(dú)立且對(duì)合①如果兩個(gè)運(yùn)動(dòng)積分的泊松括號(hào)為0，那么它們是對(duì)合的。的運(yùn)動(dòng)積分，因此在上式的積分常數(shù)β={β1,...,βn}和h這(n+1)個(gè)常數(shù)中，至少有一個(gè)不獨(dú)立，即通過(guò)h=h(β1,...,βn)相聯(lián)系，我們?nèi)ˇ耼=h。對(duì)式(9)兩邊微分得到兩個(gè)微分關(guān)系，以及生成函數(shù)的微分關(guān)系，即式(8)；再將W代入S中，可得到關(guān)系W與新舊相空間坐標(biāo)的微分關(guān)系：

現(xiàn)在我們以哈密頓特征函數(shù)W作為生成函數(shù)作正則變換(q,p)→(Q,P)，其中，

此時(shí)的新哈密頓量仍是K=H，并且新哈密頓函數(shù)只含新廣義動(dòng)量而不含新廣義坐標(biāo)，該正則變換是作為后面討論運(yùn)動(dòng)積分坐標(biāo)的基礎(chǔ)。

更進(jìn)一步，對(duì)于自治(完全)可分離的生成函數(shù)W，可寫(xiě)為：

則相應(yīng)地可以分離出n個(gè)哈雅方程(不對(duì)i求和)，可表示為：

哈雅方程方法正可方便可分離系統(tǒng)的求解。

如上所述，若自由度n的自治哈密頓系統(tǒng)為(完全)可積的，則它存在n個(gè)互相對(duì)合的運(yùn)動(dòng)積分I i(i=1,2,3,···,n)(我們用I i代表已經(jīng)從E,β中選好的運(yùn)動(dòng)積分)；此即(劉維爾)可積性定理的內(nèi)容。另一方面，根據(jù)劉維爾保體積定理，以這些運(yùn)動(dòng)積分為某一組常數(shù)時(shí)的坐標(biāo)集合，M={(q,p)|I i(q,p)=βi}，它在以I i為哈密頓函數(shù)的相流下保持不變，并且同胚于n維環(huán)面T n(若不存在共振)，因此將(完全)可積系統(tǒng)的軌道轉(zhuǎn)化為一個(gè)n維準(zhǔn)周期運(yùn)動(dòng)[63]。

關(guān)于運(yùn)動(dòng)積分，值得注意兩點(diǎn)：第一，在星系動(dòng)力學(xué)語(yǔ)境中，運(yùn)動(dòng)積分和運(yùn)動(dòng)常量是作了嚴(yán)格區(qū)分的兩個(gè)概念；第二，關(guān)于運(yùn)動(dòng)積分還有更多細(xì)分的概念。對(duì)于恒星在星系引力勢(shì)中的運(yùn)動(dòng)，給定引力場(chǎng)，運(yùn)動(dòng)常量的定義是：關(guān)于相空間坐標(biāo)和時(shí)間的、沿恒星軌道不變的任意函數(shù)，即C(x(t),v(t),t)=const。運(yùn)動(dòng)積分也是在恒星軌道中保持常數(shù)的量(即運(yùn)動(dòng)常量)，但它進(jìn)一步要求是不顯含時(shí)間的函數(shù)，I(x(t),v(t))=const。運(yùn)動(dòng)積分可以分為孤立積分和非孤立積分兩類(lèi)[8]。孤立積分能影響恒星軌道在相空間的分布，即相空間軌跡(phase trajectory)所占據(jù)空間的維度；對(duì)于可積系統(tǒng)，每增加一個(gè)孤立積分，就把相空間軌跡所占據(jù)流形(即相空間的子空間)的維度降低一維。非孤立積分則對(duì)于恒星軌道在相空間的分布沒(méi)有影響，因而沒(méi)有實(shí)用價(jià)值。

另外，在完全可積的或部分可積的系統(tǒng)中，請(qǐng)注意對(duì)合積分與孤立積分這兩個(gè)概念的區(qū)別。對(duì)于自由度為n的自治哈密頓系統(tǒng)，相空間為2n維，最多有2n?1個(gè)獨(dú)立的(即線性無(wú)關(guān)的)運(yùn)動(dòng)積分，其中孤立運(yùn)動(dòng)積分最多可以是2n?1個(gè)；而對(duì)合是與可積性有關(guān)的概念，兩兩對(duì)合的運(yùn)動(dòng)積分不超過(guò)n個(gè)(當(dāng)對(duì)合積分的個(gè)數(shù)為n時(shí)，該系統(tǒng)就是劉維爾當(dāng)年所定義的完全可積的)。因此，若有一組運(yùn)動(dòng)積分是互相對(duì)合的，那么它們肯定是孤立積分；但反之不然。例如，對(duì)于球?qū)ΨQ(chēng)引力勢(shì)中的恒星軌道(6維相空間)，運(yùn)動(dòng)積分有5個(gè)，其中孤立積分可以是4個(gè)或者5個(gè)(僅當(dāng)軌道為開(kāi)普勒型時(shí))，但相互對(duì)合的積分的個(gè)數(shù)是3[8]。對(duì)于星系動(dòng)力學(xué)中常見(jiàn)的軸對(duì)稱(chēng)引力勢(shì)系統(tǒng)(例如盤(pán)星系)，一般有3個(gè)孤立運(yùn)動(dòng)積分，能量E、角動(dòng)量的垂向(z)分量L z、以及第三個(gè)運(yùn)動(dòng)積分I3，但I(xiàn)3往往是沒(méi)有解析形式的；對(duì)于三軸引力勢(shì)系統(tǒng)，一般有3個(gè)孤立運(yùn)動(dòng)積分，能量E和I2,I3，但I(xiàn)2,I3往往是沒(méi)有解析形式的。過(guò)去，前人有很多工作就是在找第2、第3運(yùn)動(dòng)積分，以及在研究穩(wěn)態(tài)系統(tǒng)的恒星分布時(shí)所用的、以此為自變量的分布函數(shù)f(I)(這也是本綜述文章的中心內(nèi)容)。

1.2.3 作用量-角變量體系

盡管運(yùn)動(dòng)積分I與其廣義坐標(biāo)(記為θ，角變量)已經(jīng)將運(yùn)動(dòng)化為準(zhǔn)周期運(yùn)動(dòng)，但是運(yùn)動(dòng)積分一般不是正則坐標(biāo)，即它與其廣義坐標(biāo)不滿足哈密頓共軛，運(yùn)動(dòng)積分經(jīng)過(guò)變換后能夠滿足共軛的叫做作用量J。下面我們對(duì)角變量和作用量的要求作一下簡(jiǎn)介(可見(jiàn)文獻(xiàn)[8]第3.5.1節(jié))。

如前所述，我們假定存在從I到J的可逆的變換(det(?J/?I)?=0)，使作用量J滿足：(1)運(yùn)動(dòng)積分，(2)正則坐標(biāo)。于是哈密頓方程(1)可用，關(guān)于J的部分有：

于是哈密頓量只是作用量的函數(shù)，表示為：

其中，?i(J)是角頻率。對(duì)于束縛軌道，角變量θ和運(yùn)動(dòng)積分所描述的就是一個(gè)n維準(zhǔn)周期運(yùn)動(dòng)，即前述的同胚于n維環(huán)面T n的流形M；M上某一維的任一閉合環(huán)路γi(同胚于圓周S1)，由某一組(θi,J i)所代表，并且要求θ沿γi的環(huán)路積分為2π。

于是我們可以將相空間軌道展開(kāi)成傅里葉級(jí)數(shù)：

角頻率之間的比例關(guān)系豐富了體系周期類(lèi)型和遍歷性的復(fù)雜性(見(jiàn)第4章)。天體動(dòng)力系統(tǒng)的周期特性一般來(lái)說(shuō)非常強(qiáng)，這也是作用量-角變量體系在天體力學(xué)中經(jīng)久不衰的重要原因。

其中第2個(gè)等式是由于正則變換的相體積不變性；第4個(gè)等式是轉(zhuǎn)化為邊界積分，并且刨除積分奇點(diǎn)的鄰域。將奇點(diǎn)置于作用量-角變量坐標(biāo)的中心，=0，此時(shí)等于作用量的不變量。于是，我們由正則變換的龐加萊不變量的不變性，可得：

其中，T是正則變換的矩陣，指在γi對(duì)應(yīng)的相空間所在區(qū)域。也由此可見(jiàn)，某個(gè)作用量的值J i代表子環(huán)γi所圍的面積(當(dāng)對(duì)q i,p i選取合適的標(biāo)度)，而θi代表轉(zhuǎn)過(guò)的角度。

此外，我們可用式(15)計(jì)算角頻率，用式(11)中第1行的哈密頓特征函數(shù)計(jì)算角變量的值。數(shù)值計(jì)算上見(jiàn)附錄A。

于是，從一般的相空間坐標(biāo)到作用量的關(guān)系為：相空間坐標(biāo)(q,p)→運(yùn)動(dòng)積分I和角坐標(biāo)θ→正則環(huán)面坐標(biāo)(J,θ)。

1.2.4 基于守恒量的分布函數(shù)

關(guān)于單粒子分布函數(shù)f(x,v)(此時(shí)為3個(gè)自由度的，由坐標(biāo)、速度所構(gòu)成的6維μ相空間)，對(duì)于(準(zhǔn))穩(wěn)態(tài)體系(指?f/?t=0)，由無(wú)碰撞玻爾茲曼方程可推導(dǎo)，分布函數(shù)可寫(xiě)成運(yùn)動(dòng)積分的函數(shù)f=f(I)，此即Jeans定理；而對(duì)于規(guī)則軌道具有非諧振頻率的準(zhǔn)穩(wěn)態(tài)恒星系統(tǒng)的分布函數(shù)，可寫(xiě)成3個(gè)作用量的函數(shù)f=f(J)，此即強(qiáng)Jeans定理。因此人們能更加不冗余、深刻地理解恒星的分布。

早先前人所用的分布函數(shù)模型，常常是f(E,···)的函數(shù)形式，即把能量作為分布函數(shù)的參數(shù)(參見(jiàn)3.2節(jié))。不幸的是，對(duì)于實(shí)際的建模工作而言，這是一個(gè)陷阱，從而阻滯了基于運(yùn)動(dòng)積分或作用量建模這項(xiàng)事業(yè)的進(jìn)展[15]。對(duì)于實(shí)際星系的動(dòng)力學(xué)建模，不能將E放入分布函數(shù)的參數(shù)列表中。原因很簡(jiǎn)單：實(shí)際星系是一個(gè)多組分系統(tǒng)，動(dòng)力學(xué)模型中至少要包括暗物質(zhì)暈和恒星組分。軌道的能量不是局域量，如果E出現(xiàn)在參數(shù)列表中，那么該軌道周?chē)南嗫臻g密度將以不受控制的方式變化。例如，最簡(jiǎn)單的軌道是恒星位于銀河系中心的規(guī)則軌道，當(dāng)添加暗物質(zhì)暈組分后，中心引力勢(shì)會(huì)變深，位于中心的恒星的能量會(huì)減少，而該軌道保持不變。但是，只有當(dāng)模型完全組裝好并且確定了其勢(shì)能之后，我們才能知道該軌道的能量E0=Φ(0)。因此我們無(wú)法在建模工作一開(kāi)始，就指定恒星的f(E,···)。而當(dāng)分布函數(shù)只是作用量的函數(shù)時(shí)，每個(gè)軌道周?chē)南嗫臻g密度則不會(huì)改變[15,51]。

總之，相比于一般運(yùn)動(dòng)積分，作用量-角變量體系的優(yōu)勢(shì)是：作用量是一種特殊的運(yùn)動(dòng)積分，而且作為正則坐標(biāo)，作用量有共軛坐標(biāo)即角變量；作用量體系適用于使用微擾方法，也正適合描述周期性強(qiáng)的天體運(yùn)動(dòng)；作用量是絕熱不變量，可用于緩慢變化的引力勢(shì)，以及鄰近勢(shì)能中的軌道識(shí)別；此外如上段所述，使用以作用量為參數(shù)的分布函數(shù)來(lái)為星系(多組分系統(tǒng))建模，才使得基于分布函數(shù)的建模成為可能。

1.3 Stckel勢(shì)簡(jiǎn)介

作坐標(biāo)變換(x,y,z)→(λ,μ,ν)，并令

其中，τ=λ,μ,ν即構(gòu)成共焦橢球坐標(biāo)。令常數(shù)α=?a2,β=?b2,γ=?c2，并且?γ≤ν≤?β≤μ≤?α≤λ。當(dāng)τ變化時(shí)，半軸長(zhǎng)變化但焦距不變，形成一系列的共焦二次曲面。該三維坐標(biāo)系的等坐標(biāo)面如圖1所示，如式(19)中令τ=λ，由于λ+α≥0,λ+β≥0,λ+γ≥0，此時(shí)每個(gè)等λ曲面代表一個(gè)橢球面(圖1中的紅色面)；類(lèi)似地，由于μ+α≤0,μ+β≥0,μ+γ≥0，等μ曲面代表單葉雙曲面(黃色)；等ν曲面代表雙葉雙曲面(綠色)。笛卡爾坐標(biāo)系、球坐標(biāo)系等，可看作共焦橢球坐標(biāo)系的極限情形。

圖1 共焦橢球坐標(biāo)系的等坐標(biāo)面示意圖

其中，ζ,η,κ分別是相應(yīng)橢球坐標(biāo)λ,μ,ν的任意函數(shù)(見(jiàn)文獻(xiàn)[2]中的式(6))①不同文獻(xiàn)對(duì)引力勢(shì)所采用的符號(hào)有時(shí)候不相同(V、Φ或其他)，相應(yīng)地有些文獻(xiàn)中哈密頓量記為H=p2+V，有些記為H=p2+Φ；類(lèi)似地，對(duì)其分離函數(shù)Fτ(τ)前面的正負(fù)符號(hào)的規(guī)定，也可能有所出入。本文采用如式(21)所示的方式。。令Fτ(λ)=4(λ+α)(λ+β)(λ+γ)ζ(λ)，F(xiàn)τ(μ)及Fτ(ν)依此類(lèi)推②為了書(shū)寫(xiě)方便，本文后面統(tǒng)一用F(τ)代指分離函數(shù)Fτ(λ)等3個(gè)分離函數(shù)，而實(shí)際上這3個(gè)函數(shù)可以采取不同的具體形式。，那么Stckel勢(shì)可以寫(xiě)成3個(gè)分量分離的形式。

其中，pλ=，pμ及pν依此類(lèi)推。再代入哈密頓-雅克比方程H=E，以及有pτ=?W/?τ，其中W為哈密頓特征函數(shù)。解哈密頓-雅克比方程，令兩邊乘以(λ?μ)(μ?ν)(ν ? λ)，可得

由于式(21)的特殊形式以及W的可分離性，可寫(xiě)為于是式(23)中括號(hào)內(nèi)的部分只是各自τ坐標(biāo)的函數(shù)，設(shè)之為Uτ(τ)并替換，可得：

將上式兩邊對(duì)每個(gè)坐標(biāo)τ={λ,μ,ν}作兩次微分，即將作用于式(24)兩邊。然后我們從其中的對(duì)同坐標(biāo)二次微分式中可得(τ)=0，即3個(gè)Uτ(τ)是線性函數(shù)Uτ(τ)=jττ ?kτ；利用不同坐標(biāo)的二次微分式，可得3個(gè)Uτ(τ)的線性系數(shù)相等，jτ=j；將Uτ(τ)代入式(24)，以及{λ,μ,ν}取值的任意性(或者再對(duì)每個(gè)坐標(biāo)τ=λ,μ,ν作微分)，可得kτ=k相等，kτ=k，即j,k就是可分離方程的分離常數(shù)。因此3個(gè)Uτ(τ)只能是相同的線性函數(shù)

由于Uτ(τ)即式(23)的中括號(hào)內(nèi)的部分，則可得：

更方便地，H,J,K可寫(xiě)為：

其中，

此外，還有更常用的另一套運(yùn)動(dòng)積分，兩者的關(guān)系是：

此即星系動(dòng)力學(xué)中常說(shuō)的“第二個(gè)、第三個(gè)運(yùn)動(dòng)積分”I2,I3。相關(guān)推導(dǎo)見(jiàn)文獻(xiàn)[2]。類(lèi)似地，由pτ=?W/?τ，得到動(dòng)量的表達(dá)式為

積分即得作用量，通過(guò)哈密頓特征函數(shù)則可求出角變量θ。附錄A介紹了如何具體數(shù)值求解作用量和角頻率(見(jiàn)文獻(xiàn)[1,22])。

2 基于Stckel勢(shì)近似的作用量估算方法

Sanders和Binney[22,24]還提供了相應(yīng)的開(kāi)源程序TACT①https://github.com/jls713/tact。該程序用C++所寫(xiě)，能夠用若干種方法計(jì)算出作用量、角頻率和角變量，并能給出誤差；另外，還提供了相應(yīng)的Python畫(huà)圖程序。關(guān)于算法的更多細(xì)節(jié)，請(qǐng)參見(jiàn)文獻(xiàn)[20,24]，以及TACT程序的源碼。

2.1 Stckel作弊法

關(guān)于估算一般引力勢(shì)系統(tǒng)的作用量-角變量，Binney[13]提出了高效的Stckel Fudge方法，可看作是一種作弊算法(以后稱(chēng)為Stckel作弊法)。其思想是：假設(shè)引力勢(shì)可以局部地寫(xiě)成Stckel勢(shì)的形式，那么以一些技巧局部地(即在恒星軌道的某一段之內(nèi))計(jì)算目標(biāo)勢(shì)的偽分離函數(shù)(參見(jiàn)式(21)的F(τ))，并且假定其他坐標(biāo)對(duì)分離函數(shù)誤差的交叉影響較小，然后以此計(jì)算一組新的運(yùn)動(dòng)積分，得到動(dòng)量，從而通過(guò)積分得到作用量。下面先簡(jiǎn)介Binney首先提出的軸對(duì)稱(chēng)情形的Stckel作弊法，然后介紹三軸Stckel作弊法[22]，后者是對(duì)前者的擴(kuò)展。

一般來(lái)說(shuō)，軸對(duì)稱(chēng)勢(shì)下有3個(gè)孤立運(yùn)動(dòng)積分，即總能量E，角動(dòng)量垂向分量L z，以及I3。軸對(duì)稱(chēng)橢球坐標(biāo)系相對(duì)于三軸橢球坐標(biāo)系而言，一個(gè)半軸長(zhǎng)參數(shù)與另外一個(gè)相等，如果α=β，則為長(zhǎng)等軸(prolate)橢球坐標(biāo)系(λ,?,ν)，其中?代表某個(gè)子午面的方位角；如果γ=β，則為扁等軸(oblate)橢球坐標(biāo)系(λ,μ,χ)，其中χ代表某個(gè)子午面的方位角。因此在描述粒子在軸對(duì)稱(chēng)勢(shì)中的運(yùn)動(dòng)時(shí)，可以消去一個(gè)自由度。

基于等軸橢球坐標(biāo)系作出新的(u,v,?)坐標(biāo)(參考文獻(xiàn)[8]第3.5.3節(jié))，其中R=?sinhusinv,z=?coshucosv。該坐標(biāo)系下軸對(duì)稱(chēng)Stckel勢(shì)為：

若目標(biāo)引力勢(shì)Φ(u,v)寫(xiě)成Stckel勢(shì)的形式，則F u(u)?F v(v)≈(sinh2u+sin2v)ΦS(u,v)；并且若目標(biāo)引力勢(shì)接近Stckel勢(shì)，則δF u(u)對(duì)v的依賴和δF v(v)對(duì)u的依賴是很小的。于是：

同樣，F(xiàn) v(v)可以重寫(xiě)為F v(π/2)+δF v(v)。

這里u0是一個(gè)不影響作用量計(jì)算的參考值，其選取是固定v而使δF u(u)最小值時(shí)的u(見(jiàn)文獻(xiàn)[13])。F u(u0)和F v(π/2)的計(jì)算為：

本小節(jié)主要介紹Sanders和Binney[22]提出的算法，其是對(duì)Binney[13]提出的Stckel作弊法的三軸推廣。對(duì)于一般軸對(duì)稱(chēng)勢(shì)系統(tǒng)，除了Stckel作弊法，Sanders[19]還提供了對(duì)一般軸對(duì)稱(chēng)勢(shì)的Stckel擬合的方法(下一小節(jié)介紹)。但是三軸Stckel勢(shì)比軸對(duì)稱(chēng)的多一維參數(shù)，對(duì)其擬合所需的具體擬合模型會(huì)復(fù)雜許多，于是Sanders和Binney推廣出了三軸勢(shì)的Stckel作弊法。該方法借助類(lèi)似于運(yùn)動(dòng)積分J,K的另一套運(yùn)動(dòng)積分(作弊法得到的新的不獨(dú)立的運(yùn)動(dòng)積分)，以及局部近似的Stckel勢(shì)，用動(dòng)量和運(yùn)動(dòng)積分的關(guān)系式得到每一點(diǎn)的動(dòng)量pτ，再積分得到作用量J；類(lèi)似地，可求出角變量θ。本方法以當(dāng)前真實(shí)的相空間點(diǎn)做初始條件(x0,v0)，來(lái)計(jì)算一系列動(dòng)量，這樣就能不需要Stckel近似勢(shì)的分離勢(shì)F(τ)的顯式表達(dá)式；然后積分得到作用量；最小化作用量的標(biāo)準(zhǔn)差(與每個(gè)初始(x0,v0)有關(guān))，調(diào)整焦距，從而再得到合適的作用量；對(duì)于軸對(duì)稱(chēng)勢(shì)系統(tǒng)，還可做迭代以減小誤差。下面介紹其具體步驟。

對(duì)于一個(gè)一般引力勢(shì)，定義三軸情形的輔助函數(shù)：

當(dāng)固定另外兩個(gè)橢球坐標(biāo)時(shí)，某個(gè)τ對(duì)應(yīng)的Cτ,Dτ可看作常數(shù)，從而使我們可以分開(kāi)來(lái)獨(dú)立計(jì)算所需的運(yùn)動(dòng)積分。

將式(35)代入到式(26)，得到動(dòng)量的另一個(gè)計(jì)算式

與式(37)相對(duì)比，就有了這些常數(shù)之間的關(guān)系式Jτ=j?Cτ,Kτ=k+Dτ(j,k正是前面式(26)的運(yùn)動(dòng)積分)。這樣，給定一個(gè)初始相空間點(diǎn)(x0,v0)和一套合適的橢球坐標(biāo)系(坐標(biāo)變換參數(shù)即(α,β,γ))，以及引力勢(shì)各點(diǎn)的值，便可計(jì)算某橢球坐標(biāo)的相空間點(diǎn)時(shí)的運(yùn)動(dòng)積分[22]。

下面介紹Jτ,Kτ的計(jì)算。由式(38)得：

Kτ依賴于Jτ，而其他量可以計(jì)算。當(dāng)固定另外兩組橢球坐標(biāo)和動(dòng)量時(shí)，由于Stckel勢(shì)的分離性，可以通過(guò)僅改變某個(gè)τ和pτ的值，便可使粒子仍保持在既定軌道上，從而認(rèn)為在Stckel勢(shì)下保持橢球相空間坐標(biāo)不變，哈密頓量對(duì)某個(gè)橢球坐標(biāo)的偏導(dǎo)為0(見(jiàn)文獻(xiàn)[22]中的第3章)，即：

分別將λ和λ0代入式(39)，并作差，除以Pλ，可推得

用式(41)來(lái)計(jì)算式(40)第一項(xiàng)，可見(jiàn)同時(shí)可以消去?Φ/?λ，因此得到Jλ以及類(lèi)似的Jμ,Jν的表達(dá)式：

于是得到6個(gè)估算的運(yùn)動(dòng)積分Jτ,Kτ的表達(dá)式，其中只有3個(gè)是獨(dú)立的[22]。

(2)合適焦距的選取

有了這些新的運(yùn)動(dòng)積分，可用式(38)求得相應(yīng)τ點(diǎn)的動(dòng)量，之后積分得到作用量。軸對(duì)稱(chēng)情形也類(lèi)似。然而，坐標(biāo)變換的參數(shù)的焦距任取γ=?1kpc2)，其大小非常影響該方法所估算的作用量的方差(認(rèn)為是誤差)。Sanders和Binney給出了數(shù)值調(diào)整合適焦距的方法。例如對(duì)于短軸閉環(huán)軌道(short axis closed loop orbits)，軌道形成在z=0面，焦距為y=±?1的橢圓。在這種軌道，只有作用量Jν不為0，利用這一點(diǎn)，給每個(gè)能量值找出一個(gè)合適的焦距。

其步驟是，先給定一個(gè)一般勢(shì)，最大√和最小能量值以及期間的能量網(wǎng)格點(diǎn)；在每個(gè)網(wǎng)格點(diǎn)E i，畫(huà)一個(gè)從坐標(biāo)(0,y k,0)以速度為垂直出發(fā)的半周期的軌道；記錄所畫(huà)軌道的相空間點(diǎn)并計(jì)算作用量Jν，求出使作用量的方差最小的焦距?1。?2的計(jì)算也可用上述方法求出。

(3)環(huán)面迭代增加精度

2.2 Stckel擬合法

這里介紹Sanders[19]提出的長(zhǎng)等軸橢球坐標(biāo)系下的一般軸對(duì)稱(chēng)勢(shì)的Stckel擬合法。該方法是對(duì)一般引力勢(shì)進(jìn)行Stckel形式的擬合，即采樣軌道來(lái)擬合出一個(gè)合適的Stckel引力勢(shì)的分離函數(shù)F(τ)，來(lái)得到運(yùn)動(dòng)積分。Sanders和Binney[24]把Stckel作弊法和這種Stckel擬合法歸類(lèi)為“非收斂方法”，因?yàn)檫@兩種方法依賴于做一定的近似，而近似所導(dǎo)致的誤差不能通過(guò)增加計(jì)算量來(lái)不斷減小。

在某個(gè)子午面上，運(yùn)動(dòng)的位形空間從柱坐標(biāo)(R,?,z)變換到扁球坐標(biāo)(λ,?,ν)依賴于

其中，λ,ν是上式作為關(guān)于τ的方程根，a2=α,b2=β,c2=γ，類(lèi)似于三軸情形下的式(19)。此時(shí)Stckel勢(shì)的形式簡(jiǎn)化為：

軸對(duì)稱(chēng)時(shí)第三運(yùn)動(dòng)積分的解析表達(dá)式為：

動(dòng)量和運(yùn)動(dòng)積分的關(guān)系表示為：

(1)給定初始的(x,v)，以及真實(shí)的勢(shì)(軌道積分中所感受到的勢(shì))，并在真實(shí)的勢(shì)中計(jì)算實(shí)際的運(yùn)動(dòng)積分E,L z。

(2)焦距的選取。作上述變換到橢球坐標(biāo)。然而，變換的參數(shù)，即焦距?=a2?c2(因?yàn)閍=b，只要一個(gè)a2?c2便可確定橢球坐標(biāo)系，可取c2=1)，非常影響作用量的標(biāo)準(zhǔn)差，因此通過(guò)一定的手續(xù)來(lái)估算一個(gè)更好的焦距?。Sanders[19]采用的估計(jì)方法是，由于Stckel勢(shì)的分離性，可得

由上式以及兩種坐標(biāo)之間的變換關(guān)系，可推得[19]：

既然假設(shè)所估算的局部勢(shì)Φ就是該Stckel的勢(shì)ΦS，于是就用該式來(lái)估計(jì)a2的值，其中的偏導(dǎo)用已知的真實(shí)Φ來(lái)數(shù)值計(jì)算梯度而得；每個(gè)a2值，需要在軌道中每個(gè)點(diǎn)的相應(yīng)值來(lái)計(jì)算，最后做平均，這些點(diǎn)在下一步的軌道采樣中獲得。

(3)軌道采樣以及運(yùn)動(dòng)區(qū)域的估算。進(jìn)行軌道積分，即在(x,v)坐標(biāo)下通過(guò)每一點(diǎn)的勢(shì)和時(shí)間積分來(lái)算下一步的坐標(biāo)(x′,v′)，即相當(dāng)于數(shù)值法畫(huà)軌道。Sanders和Binney的TACT程序中用守恒格式的數(shù)值方法畫(huà)出軌道，并且每個(gè)軌道周期離散為300個(gè)點(diǎn)。

在畫(huà)軌道時(shí)：1)在前幾周期的一些點(diǎn)中，隔5個(gè)點(diǎn)用式(48)估算a2(c2已取1)并平均，最后得到那個(gè)最合適的扁球坐標(biāo)；2)通過(guò)前后?τ/?t異號(hào)等條件(其中每次需要x到τ的坐標(biāo)變換)來(lái)判斷是否達(dá)到軌道的τ的端點(diǎn)，找到后就記錄，得到λ+?,ν+,ν?=c2)，從而得到軌道在橢球坐標(biāo)的所在區(qū)域。

假如所擬合的勢(shì)Φ是Stckel的時(shí)，會(huì)有?(λ,ν)Φ(λ,ν)=F(λ)?F(ν)，正是據(jù)此得到若干插值點(diǎn)局部的假設(shè)是Stckel形式的所求引力勢(shì)的值。關(guān)于F(τ)的選取，其要求也正是使其與Stckel勢(shì)的偏差

最小。歸一化的權(quán)重函數(shù)形式為：

這個(gè)權(quán)重函數(shù)是提前選的，相關(guān)內(nèi)容可見(jiàn)文獻(xiàn)[1,4]。于是擬合中所選用的Stckel勢(shì)的具體形式為：

其中，

這是因?yàn)樵谶@種情況下偏差是最小的。證明方法是，式(49)取最小時(shí)，對(duì)于λ,ν分別成立，即將其分別對(duì)λ,ν求導(dǎo)應(yīng)為0，得

再做一些代換即可得到式(50)。

接下來(lái)，有了式(50)的形式后，便可結(jié)合式(51)和(52)直接積分得到某個(gè)點(diǎn)(λ,ν)的F(τ)，如此選取若干個(gè)坐標(biāo)點(diǎn)(TACT程序選取了40個(gè)點(diǎn))，用來(lái)三次樣條插值得到所擬合的勢(shì)(樣條插值使擬合的勢(shì)會(huì)更平滑)。

(5)計(jì)算運(yùn)動(dòng)積分。借用這個(gè)擬合勢(shì)，以及之前的E,L z，用第三運(yùn)動(dòng)積分的式(45)在3個(gè)端點(diǎn)τ+?計(jì)算I3并平均，從而得到3個(gè)相應(yīng)運(yùn)動(dòng)積分。其中計(jì)算時(shí)需用真實(shí)性條件>0進(jìn)行判斷，如果不符合，就不用端點(diǎn)而用初始相空間點(diǎn)。

(6)計(jì)算作用量-角變量。計(jì)算的F(τ)除了上一步用來(lái)算I3，還用來(lái)在式(46)中以此得到當(dāng)前點(diǎn)的pτ，從而計(jì)算作用量J和角變量θ(數(shù)值計(jì)算步驟具體見(jiàn)附錄A)。

早在1985年，de Zeeuw和Lynden-Bell[2]就已提出了把Stckel勢(shì)作為模型，來(lái)局部擬合和全局?jǐn)M合星系引力勢(shì)的思想和方法。前者在具體星系引力勢(shì)的平衡點(diǎn)附近展開(kāi)引力勢(shì)，然后通過(guò)對(duì)應(yīng)Stckel的橢球坐標(biāo)形勢(shì)來(lái)確定參數(shù)；后者通過(guò)以平均的方法來(lái)近似Stckel勢(shì)的3個(gè)分離函數(shù)，從而計(jì)算運(yùn)動(dòng)積分。下面簡(jiǎn)單介紹。

(1)局部擬合

在平衡位置附近將V=V(x2,y2,z2)形式的引力勢(shì)展開(kāi)為：

其中，指標(biāo)k,l,m只取偶數(shù)值。另外，Stckel勢(shì)可寫(xiě)為式(21)的分離形式，以(λ,μ,ν)=(?α,?β,?γ)作為展開(kāi)源點(diǎn)，展開(kāi)分離的函數(shù)選取在展開(kāi)源點(diǎn)附近的引力勢(shì)為0，即ζ?1=η?1=κ?1=0。于是可得：

推導(dǎo)中用到了(λ+α)的展開(kāi)式，比如：

文獻(xiàn)[2]還討論了這些系數(shù)的約束條件、存在性和合適形式的選取。

(2)全局?jǐn)M合

文獻(xiàn)[2]的第4章對(duì)橢球勢(shì)做了討論，表明與橢圓星系有關(guān)的引力勢(shì)不僅可以局部地，而且還可以在全局地近似為St¨ckel形式，并且給出了橢球勢(shì)(具體為鏡像對(duì)稱(chēng)的引力勢(shì)V=V(x2,y2,z2))的全局?jǐn)M合方法。該方法要求所擬合的勢(shì)在整體上離Stckel勢(shì)最近，并且當(dāng)V恰好是Stckel勢(shì)時(shí)應(yīng)得到精確的結(jié)果。下面列出其方法的主要步驟。

對(duì)于任意函數(shù)Q=Q(λ,μ,ν)定義下列平均

其中，Λ(λ),M(μ),N(ν)是根據(jù)積分的約束選擇的權(quán)重函數(shù)，相應(yīng)的歸一化常數(shù)為：

定義輔助函數(shù)

那么，分離函數(shù)可計(jì)算為：

2.3 基于環(huán)面映射的方法

本節(jié)介紹的內(nèi)容包括：(1)McGill和Binney[11]提出、Binney及其合作者后來(lái)完善的，通過(guò)正則變換構(gòu)造環(huán)面(作用量-角變量)來(lái)計(jì)算作用量的方法，即基于玩具模型的環(huán)面映射(torus mapping,TM)思想及其TM程序[12]；(2)基于此思想的兩個(gè)擴(kuò)展：迭代環(huán)面構(gòu)造法(iterative torus construction,ItTC)[22,24]和軌道積分?jǐn)M合生成函數(shù)法(orbit integration to generation function,O2GF)[21,24]。這類(lèi)方法的特點(diǎn)是精度最高，計(jì)算花費(fèi)大，屬于收斂方法，即可以通過(guò)充分大的計(jì)算量，把作用量的誤差一直不斷地減小下去。

將環(huán)面(n-Torus)視為幾何物體，而將解析的哈密頓量的環(huán)面變形為所求的目標(biāo)哈密頓量的環(huán)面，如同微擾理論，通過(guò)尋求從玩具環(huán)面坐標(biāo)到目標(biāo)環(huán)面坐標(biāo)的正則變換來(lái)實(shí)現(xiàn)[11]。即用已知的解析引力模型作為玩具模型，來(lái)逼近一般的目標(biāo)引力系統(tǒng)，而這是在環(huán)面坐標(biāo)的意義上進(jìn)行的，并且通過(guò)正則變換來(lái)實(shí)現(xiàn)，此即環(huán)面映射的關(guān)鍵思想。一般所用的玩具模型為等時(shí)勢(shì)或諧振子勢(shì)，如式(70)和(71)，這些都是Stckel勢(shì)，可解析求解。盡管環(huán)面映射法不同于前述的通過(guò)近似Stckel的分離函數(shù)Fτ(τ)的Stckel作弊法，但其需要從某種作為玩具模型的Stckel勢(shì)出發(fā)，本質(zhì)上相當(dāng)于某種簡(jiǎn)化版的微擾法，因而仍可看作是Stckel勢(shì)的應(yīng)用拓展。

本節(jié)中，若無(wú)特殊說(shuō)明，用(Jtoy,θtoy)代表玩具環(huán)面(即玩具勢(shì)下粒子運(yùn)動(dòng)的環(huán)面坐標(biāo))，(J,θ)代表目標(biāo)環(huán)面(即所求的一般引力勢(shì)下粒子運(yùn)動(dòng)的環(huán)面坐標(biāo))。從玩具環(huán)面映射到目標(biāo)環(huán)面，將此看作一個(gè)正則變換，選取第二類(lèi)生成函數(shù)S=S(J,θtoy)，將這個(gè)生成函數(shù)對(duì)角變量傅里葉展開(kāi)，其一般形式為：

其中，這里的虛單位i是為后面求導(dǎo)表達(dá)式的簡(jiǎn)便而加的。當(dāng)沿著玩具環(huán)面在θ的任一分量增加一個(gè)周期2π，相應(yīng)地，目標(biāo)圓環(huán)上的像點(diǎn)也應(yīng)完成一個(gè)周期回路，最終這要求A=J。于是

通過(guò)作用量-角變量的實(shí)數(shù)性以及哈密頓量的時(shí)間反演不變性，可以對(duì)生成函數(shù)的系數(shù)進(jìn)行約束[11]：

于是有[21]：

其中，向量N={(i,j,k)}代表一系列精度的坐標(biāo)(i,j,k)集，(i,j,k)滿足(k>0)或(k=0,j>0)或(k=0,j=0,i>0)[12,21]。對(duì)這個(gè)生成函數(shù)S(θtoy,J)求微分，得

這樣就通過(guò)生成函數(shù)聯(lián)系起了玩具環(huán)面和目標(biāo)環(huán)面。

2.3.1 環(huán)面映射程序

Binney和McMillan[12]提供了TM程序(Torus Mapper)①TM程序的源碼可下載于：https://github.com/PaulMc Millan-Astro/Torus。，該程序?qū)崿F(xiàn)三個(gè)方面的功能，具體功能與輸入輸出如下：用戶寫(xiě)出目標(biāo)引力勢(shì)及其一階微分，以及給定作用量(J r,J?,J z)，1)對(duì)于特定θ，可計(jì)算(x,v)；2)計(jì)算雅克比?(x)/?(θ)；3)給定(x,J)，計(jì)算θ。與傳統(tǒng)的恒星軌道相比，環(huán)面的優(yōu)點(diǎn)是：一個(gè)環(huán)面是由100個(gè)數(shù)字指定的，即S n,?S n/?J和一些其他參數(shù)，而不是由沿時(shí)間序列的數(shù)千個(gè)相空間坐標(biāo)指定，因此可以在給定分辨率下顯著縮小軌道庫(kù)的大小[12]。

環(huán)面映射方法的流程見(jiàn)式(68)[12]，通過(guò)三步正則變換完成，變換形式為：

已知的真實(shí)勢(shì)和軌道，給定目標(biāo)作用量J和系數(shù)S n的試探值，由玩具坐標(biāo)θtoy的網(wǎng)格點(diǎn)，根據(jù)式(66)可得到相應(yīng)的玩具作用量Jtoy。在玩具勢(shì)下，可直接通過(guò)解析的方法從環(huán)面坐標(biāo)得到相空間坐標(biāo)(xtoy,vtoy)。有了這些相空間坐標(biāo)可計(jì)算相應(yīng)采樣點(diǎn)的能量，之后用Levenberg-Marquardt算法迭代以最小化哈密頓量偏差δH，得到最終的S n。類(lèi)似地，可得到?S n/?J，再通過(guò)式(67)可得到目標(biāo)角變量θ[24]。

這樣，通過(guò)真實(shí)的環(huán)面坐標(biāo)(目標(biāo)環(huán)面坐標(biāo)作為已知)，用合適的系數(shù)S n(J)和?S n/?J，轉(zhuǎn)化到玩具環(huán)面坐標(biāo)；再通過(guò)玩具勢(shì)下的解析性就得到相空間坐標(biāo)(xtoy,vtoy)，并由點(diǎn)變換得到(x,v)。點(diǎn)變換也是一類(lèi)正則變換，在程序中僅當(dāng)J z/J r>0.05且嘗試在不進(jìn)行點(diǎn)變換的情況下擬合環(huán)面失敗時(shí)，再使用點(diǎn)變換[12]。

文獻(xiàn)[12]中的附錄圖A1給出了環(huán)面映射方法的原理示意，對(duì)于一個(gè)好的玩具模型，更高階的生成函數(shù)傅里葉系數(shù)能越來(lái)越好地?cái)M合軌道。

2.3.2 迭代環(huán)面構(gòu)造法

如上所述，TM程序[11,12]給出了一般軸對(duì)稱(chēng)系統(tǒng)的相空間坐標(biāo)和環(huán)面坐標(biāo)轉(zhuǎn)換的功能，主要是由(J,θ)計(jì)算(x,v)。對(duì)于環(huán)面坐標(biāo)的計(jì)算，除了按照上述方法，重復(fù)地構(gòu)造環(huán)面，直到得到的相空間坐標(biāo)穿過(guò)給定相空間坐標(biāo)，TACT程序的其中一個(gè)方法對(duì)其進(jìn)行了改進(jìn)，即迭代環(huán)面構(gòu)造法(ItTC)[22,24]。該方法是“Stckel作弊法+TM程序”的迭代算法，即用Stckel作弊法計(jì)算試探的環(huán)面坐標(biāo)，然后尋找環(huán)面上與其最近的相空間點(diǎn)，最后補(bǔ)上殘差并迭代[22]。

給定一個(gè)初始相空間點(diǎn)(x,v)，先用Stckel作弊法計(jì)算一個(gè)環(huán)面坐標(biāo)(J S,θS)；然后用環(huán)面映射得到一個(gè)以此為作用量的環(huán)面，該環(huán)面的頻率為?，在該環(huán)面上找到一個(gè)與給定點(diǎn)最近的相空間點(diǎn)(x S,v S)，即

使其最小。之后，計(jì)算所找到的這一點(diǎn)的環(huán)面坐標(biāo)J P(x S,v S),θP(x S,v S)，得到?J=J P?J S作為誤差，假設(shè)這個(gè)誤差變化很慢，于是將其作為殘差補(bǔ)回去，得到一個(gè)更好的作用量J2=J S+?J。

可以再次構(gòu)造一個(gè)以J2為作用量的環(huán)面，然后繼續(xù)從中找到一個(gè)與給定初始點(diǎn)最近的相空間點(diǎn)，……，如此往復(fù)迭代，得到更精確的作用量，這個(gè)過(guò)程是收斂的。

2.3.3 軌道積分?jǐn)M合生成函數(shù)法

這里介紹Sanders和Binney[21,24]給出的軌道積分?jǐn)M合生成函數(shù)法。該方法過(guò)程為：軌道時(shí)序采樣(x,v)(t i)→計(jì)算玩具環(huán)面坐標(biāo)(Jtoy,θtoy)(t i)→解生成函數(shù)系數(shù)矩陣。它和TM程序的區(qū)別是，該方法將環(huán)面映射思想應(yīng)用于一段運(yùn)動(dòng)過(guò)程(即一段時(shí)間序列)，計(jì)算出在時(shí)間t上沿軌道采樣的一系列相空間坐標(biāo)的作用量、角變量和角頻率，即構(gòu)造生成函數(shù)是基于軌道積分，而不是TM所基于的某個(gè)時(shí)間點(diǎn)的哈密頓量偏差δH的最小二乘。

具體步驟為，給定一個(gè)真實(shí)相空間點(diǎn)作為初始點(diǎn)，首先計(jì)算一個(gè)時(shí)間N TT的軌道積分(畫(huà)軌道)，并采集Nsamp個(gè)相空間采樣點(diǎn)，其中T是具有相同能量的大圓軌道的周期；然后在每個(gè)相空間點(diǎn)，解析計(jì)算其玩具勢(shì)中的(θtoy,Jtoy)；之后將這些玩具環(huán)面坐標(biāo)點(diǎn)代入到式(66)和(67)來(lái)求解整個(gè)方程組的未知量(目標(biāo)作用量J，生成函數(shù)傅里葉系數(shù)S n及其導(dǎo)數(shù)?S n/?J，以及目標(biāo)角頻率?和常數(shù)θ(0))。

(1)軌道積分與軌道分類(lèi)

進(jìn)行軌道積分。在沒(méi)有轉(zhuǎn)動(dòng)的情況下，一般三軸勢(shì)允許兩類(lèi)基本的非共振軌道：環(huán)軌道和盒軌道[1,21]。同一真實(shí)勢(shì)，不同的初始條件可能導(dǎo)致不同類(lèi)型的軌道，于是TACT為不同類(lèi)型的軌道匹配相應(yīng)的玩具勢(shì)，與給定軌道類(lèi)的環(huán)面相同的幾何結(jié)構(gòu)。

(2)玩具勢(shì)的擬合

對(duì)于環(huán)軌道，玩具勢(shì)為等時(shí)勢(shì)：

對(duì)于盒軌道，玩具勢(shì)為諧振子勢(shì)：

其判定依據(jù)為軌道類(lèi)型，即看L x和L z的方向變化。在用玩具引力勢(shì)之前，需要在軌道上采集樣點(diǎn)，擬合出相應(yīng)玩具引力勢(shì)的質(zhì)量和標(biāo)長(zhǎng)參數(shù)。

(3)生成函數(shù)

在每個(gè)相空間樣點(diǎn)t i計(jì)算相應(yīng)的玩具作用量-角變量(θtoy(t i),Jtoy(t i))，并產(chǎn)生相應(yīng)的式(66)，其中各J和S n未知。Sanders和Binney軌道積分?jǐn)M合生成函數(shù)法的實(shí)現(xiàn)技巧為：“我們不能精確地解這些方程，因?yàn)樘幚淼姆匠逃袩o(wú)窮多個(gè)未知數(shù)。我們可以在每個(gè)方程的右邊只包含有限個(gè)項(xiàng)，則方程右邊和左邊不能精確一致，因此可以采取的正確方法是：對(duì)每個(gè)方程的左右邊的差作平方，然后對(duì)所有方程的差平方求和(即軌道積分)，并最小化這個(gè)總的平方和”[21]。其總平方和的表達(dá)式是：

其中，k=1,2,3代表坐標(biāo)維數(shù)，N限制為有限維向量(有限精度)，|n|≤Nmax，Nmax≈6；(t i)代表樣點(diǎn)，最外層的求和(關(guān)于指標(biāo)i)是遍歷所有采樣點(diǎn)，此即“軌道積分”的含義。更多關(guān)于樣點(diǎn)個(gè)數(shù)的選取原則請(qǐng)?jiān)斠?jiàn)參考文獻(xiàn)[21]中的2.2和2.3節(jié)。

下面最小化E1，即令的偏導(dǎo)為0，

這些方程可轉(zhuǎn)化為一個(gè)矩陣方程組

其中，

其中，c nk(t i)為N乘3的矩陣，x Jtoy和b Jtoy為(N+3)維的列向量，I3為單位矩陣，A Jtoy為對(duì)稱(chēng)矩陣。最終求解矩陣方程得到目標(biāo)作用量J和生成函數(shù)傅里葉系數(shù)S n。

類(lèi)似地，對(duì)于求角變量θ′，則最小化

得到角頻率?′以及θ(0)和?S n/?J。

此外，對(duì)式(66)中的玩具作用量Jtoy在玩具角變量θtoy空間作平均，得到目標(biāo)作用量

便是一種更簡(jiǎn)單的方法——平均生成函數(shù)法(average generation function,AvGF)[24,29]。

2.4 TACT所提供的方法的效果對(duì)比

我們運(yùn)行TACT程序[23,24]用幾種方法分別估算出作用量，如圖2所示。我們所使用的引力勢(shì)為某種軸對(duì)稱(chēng)引力勢(shì)，即McMillan[57]擬合的最佳銀河系質(zhì)量分布模型，分為薄厚兩個(gè)恒星盤(pán)、HI和H2兩個(gè)氣體盤(pán)以及核球和暗物質(zhì)暈；各個(gè)星系成分的模型輸入?yún)?shù)(適合TACT程序的格式)，已在作者的github主頁(yè)①https://github.com/PaulMc Millan-Astro/GalPot/blob/master/pot/PJM17_best.Tp ot。給出。我們采用的恒星軌道為典型的銀河系薄盤(pán)軌道，初始條件為(x,v)=(8.178,0.1,0.225;20.13,187.01,14.95)，單位為kpc和km/s。圖中的幾種方法分別是：Fudge v1和Fudge v2為Stckel Fudge作弊法，其中，F(xiàn)udge v1用式(48)估算焦距，F(xiàn)udge v2基于殼軌道(shell orbits)來(lái)估算焦距；Fit對(duì)應(yīng)Stckel勢(shì)擬合法，ItTC對(duì)應(yīng)環(huán)面迭代法，O2GF對(duì)應(yīng)軌道積分?jǐn)M合生成函數(shù)法，AvGF為平均生成函數(shù)法；CAA和SAA分別為柱絕熱近似和球絕熱近似方法，下一節(jié)作簡(jiǎn)介。軌道中的每個(gè)點(diǎn)都獨(dú)立計(jì)算一次作用量，因而作用量的數(shù)值隨時(shí)間有所波動(dòng)。

圖2 作用量數(shù)值計(jì)算的效果對(duì)比

Sanders和Binney[24]對(duì)幾種方法的在軌道類(lèi)型、精度和時(shí)間等方面的性能做了更詳細(xì)的對(duì)比(見(jiàn)其中的圖2到圖6)?？傮w上看，收斂方法比非收斂方法有更高的精度，但耗時(shí)也高1～3個(gè)量級(jí)。其中，軌道積分?jǐn)M合生成函數(shù)法精度最高，耗時(shí)較長(zhǎng)；Stckel作弊法速度很快，是大量計(jì)算的最經(jīng)濟(jì)方法，但精度不如收斂性方法高，特別是對(duì)于盒軌道的誤差大。另外，對(duì)于軸對(duì)稱(chēng)勢(shì)經(jīng)過(guò)環(huán)面迭代后(即ItTC法，目前TM程序只提供了軸對(duì)稱(chēng)的環(huán)面迭代)，Stckel作弊法的精度會(huì)提高，而耗時(shí)也變長(zhǎng)。此外，TACT程序包中，軌道積分?jǐn)M合生成函數(shù)法和Stckel作弊法還能用于三軸勢(shì)。

2.5 其他一些求運(yùn)動(dòng)積分的方法或程序

前面兩節(jié)我們主要介紹了Binney及其合作者們近年來(lái)的工作。接下來(lái)我們介紹其他一些基于Stckel勢(shì)的估算運(yùn)動(dòng)積分或作用量的方法。

2.5.1 “近似的運(yùn)動(dòng)積分”解析式

若Φ在(?α,?β,?γ)連續(xù)，應(yīng)有Φ(?α,?β,?γ)=0。讓勢(shì)兩邊取(λ,?β,?γ)位置的值，則有

同理有

將上面的分離函數(shù)代回到式(78)，可得

這樣就能分別用一個(gè)“橢球坐標(biāo)軸”(如(λ,?β,?γ)中后兩個(gè)坐標(biāo)固定)上的勢(shì)值來(lái)表示任意點(diǎn)的引力勢(shì)。

采取類(lèi)似的手法，這個(gè)引力勢(shì)下的運(yùn)動(dòng)積分可以解析性地表達(dá)出來(lái)。運(yùn)動(dòng)積分J,K的表達(dá)式為

相關(guān)推導(dǎo)見(jiàn)參考文獻(xiàn)[2]中的式(9)和(14)。將式(79)和(80)代入，可計(jì)算得到與式(78)同樣方式的運(yùn)動(dòng)積分J,K的表達(dá)式。

類(lèi)似地，由式(15)，I2,I3(第一個(gè)運(yùn)動(dòng)積分是能量E)的表達(dá)式可寫(xiě)為：

即參考文獻(xiàn)[54]中的式(9)和(10)。文獻(xiàn)[54]推導(dǎo)了上式中的橢球坐標(biāo)的函數(shù)項(xiàng)，得到

其中，Φ(λ,μ,ν)即為式(78)中的形式，其簡(jiǎn)化形式見(jiàn)文獻(xiàn)[54]中的式(16)和(17)。

該方法計(jì)算“每一點(diǎn)的”運(yùn)動(dòng)積分，需要當(dāng)前點(diǎn)的坐標(biāo)、速度和角動(dòng)量等軌道信息。文獻(xiàn)[54]在后文用高度精確的數(shù)值方法計(jì)算出軌道，并研究了這種方法所計(jì)算的運(yùn)動(dòng)積分的誤差，最后做了Jeans方程的檢驗(yàn)。

2.5.2 絕熱近似

對(duì)于軸對(duì)稱(chēng)勢(shì)，文獻(xiàn)[24]還介紹了用近似公式來(lái)估算作用量的方法。通過(guò)假設(shè)在柱坐標(biāo)系或球坐標(biāo)系中引力勢(shì)是可分離的，可求出這些作用量。

(1)柱絕熱近似[45]。這種觀點(diǎn)認(rèn)為，盤(pán)軌道恒星的垂向(即z方向)頻率?z比其徑向頻率?r大得多，因此當(dāng)恒星徑向振蕩時(shí)，主導(dǎo)其垂向振蕩的勢(shì)可被認(rèn)為變化緩慢，其結(jié)果就是J z是絕熱不變的[24]。用這些假設(shè)推導(dǎo)出作用量的近似公式

以及計(jì)算積分上下限時(shí)需要用到的垂直能量(即垂直動(dòng)能加垂直勢(shì)能，其中前者是由于速度的平方按正交坐標(biāo)方向分為三個(gè)分量平方的相加；后者是保守力按照某方向的分解所得的勢(shì)能)和有效勢(shì)：

其中，R c是導(dǎo)心半徑(可參見(jiàn)有關(guān)“本輪近似”的文獻(xiàn))。

(2)球絕熱近似假設(shè)沿長(zhǎng)球坐標(biāo)的球面存在絕熱不變振蕩。作用量的近似公式為：

其中，

具體計(jì)算可見(jiàn)參考文獻(xiàn)[24]。

2.5.3 一些相關(guān)程序

除了Sanders和Binney[24]的TACT程序包以及其環(huán)面迭代功能所依賴的TM程序包，星系動(dòng)力學(xué)還有一些程序可供使用，這些程序現(xiàn)在都可以公開(kāi)獲取。主要程序如下所述。

(1)Gadget，Gadget①https://wwwmpa.mpa-garching.mpg.de/galform/gadget/index.shtml以及其后續(xù)發(fā)展Arepo②https://arep o-code.org程序，是“引力N體+流體力學(xué)”模擬程序。在星系動(dòng)力學(xué)領(lǐng)域，我們可以用它來(lái)做無(wú)碰撞N體模擬(特別是當(dāng)我們考察暗物質(zhì)粒子構(gòu)成的暗物質(zhì)暈的動(dòng)力學(xué)時(shí))。

(2)Galpy，Galpy(a Python library for galactic dynamics)③https://github.com/jobovy/galpy是Bovy[30]提供的Python程序包，能用來(lái)進(jìn)行數(shù)值軌道積分、作用量-角變量、分布函數(shù)等相關(guān)計(jì)算。

(3)AGAMA，AGAMA(action-based galaxy modelling architecture)④https://github.com/Galactic Dynamics-Oxford/Agama是Vasiliev[31]提供的星系動(dòng)力學(xué)程序，能提供任意解析密度剖面或N體模型的引力勢(shì)的計(jì)算、軌道積分和分析、位置-速度與作用-角度變量之間的轉(zhuǎn)換、分布函數(shù)的計(jì)算以及迭代構(gòu)造自洽多組分星系模型等功能。

2.5.4 深度學(xué)習(xí)在作用量計(jì)算中的應(yīng)用

在這一小節(jié)中，我們介紹最新的一個(gè)研究動(dòng)態(tài)。Ibata等人[32]將人工智能(AI)中的深度學(xué)習(xí)方法應(yīng)用到作用量計(jì)算上。該文用AI算法來(lái)實(shí)現(xiàn)在靜態(tài)引力勢(shì)下從粒子在相空間的“位置-速度”到“作用量-角變量”的映射，并且得到加速度場(chǎng)。該算法的核心原理與前述“生成函數(shù)法”類(lèi)似，先從數(shù)據(jù)點(diǎn)中擬合玩具勢(shì)；然后通過(guò)從玩具環(huán)面(θ,J)到目標(biāo)環(huán)面(θ′,J′)的正則變換表達(dá)式，迭代出最合適的生成函數(shù)G；之后是正則點(diǎn)變換的生成函數(shù)P，最終得到所求的目標(biāo)環(huán)面(θ′′,J′′)。其中兩個(gè)正則變換所用網(wǎng)絡(luò)的形式為：

該算法的程序?yàn)锳CTIONFINDER，基于Python的Pytorch所寫(xiě)。據(jù)文獻(xiàn)[32]介紹，經(jīng)測(cè)試ACTIONFINDER得到的作用量精度要比Stckel作弊法要高，而且它還有如下優(yōu)點(diǎn)：相比于早期的Torus Mapper對(duì)初值要求更不敏感；不需要向它提供系統(tǒng)內(nèi)軌道的公平樣本，幾乎可覆蓋任何感興趣區(qū)域的樣本；不需要預(yù)先知道所研究系統(tǒng)的哈密頓量或引力勢(shì)；此外，該算法所用的深層神經(jīng)網(wǎng)絡(luò)，相比前述“生成函數(shù)法”線性分解成的傅里葉系數(shù)，能夠生成更加靈活的非線性函數(shù)，并且這種靈活性可能使其更容易擬合更一般的動(dòng)力系統(tǒng)。看來(lái)，將來(lái)經(jīng)過(guò)一定的發(fā)展完善后，深度學(xué)習(xí)或其他AI技術(shù)用來(lái)計(jì)算作用量，適用性會(huì)比傳統(tǒng)方法更廣泛、性能更佳，應(yīng)該有一定的發(fā)展前景，未來(lái)或成為分析更復(fù)雜體系的另一種嘗試。

3 應(yīng)用

在星系研究中，作用量的一個(gè)重要應(yīng)用在分布函數(shù)方面，即使用以作用量作為自變量的分布函數(shù)，為星系建模[15]。一般來(lái)說(shuō)，我們難以直接獲得分布函數(shù)的形式。在建模以及數(shù)據(jù)擬合的實(shí)際工作中，人們通常是為不同的星系組分(如恒星盤(pán)、核球、暗物質(zhì)暈等)采用某種先驗(yàn)的分布函數(shù)形式，然后根據(jù)數(shù)據(jù)擬合來(lái)確定其中的具體參數(shù)的取值。

大約10年前，基于分布函數(shù)的建模方法實(shí)際上很難應(yīng)用(見(jiàn)Read[61]和Binney[15]的綜述)，星系的建模工作往往是基于Jeans方程，即對(duì)分布函數(shù)所滿足的無(wú)碰撞玻爾茲曼方程取矩(即所謂“矩方法”)。假定星系處于穩(wěn)態(tài)(即?f/?t=0)，從無(wú)碰撞玻爾茲曼方程出發(fā)，分別乘上速度分量并且積分掉速度空間(即取矩)，即得到Jeans方程。對(duì)于實(shí)際的銀河系建模工作，研究者還會(huì)假設(shè)軸對(duì)稱(chēng)(即在柱極坐標(biāo)下?f/??=0,?f/?v?=0)，那么Jeans方程的表達(dá)式為(柱極坐標(biāo))：

其中，ν為示蹤者恒星的數(shù)密度，為分布函數(shù)的0階矩；平均速度為1階矩；速度彌散張量為 ∫

為2階矩。在觀測(cè)上，速度彌散可由某一個(gè)位置處諸恒星的速度方差來(lái)計(jì)算。

正如Read所總結(jié)[61]，Jeans方程和矩方法的優(yōu)點(diǎn)是：與其他方法相比非?？?，可以探索較大的參數(shù)空間；并且不需要假設(shè)分布函數(shù)的形式，因?yàn)橹皇窍拗屏怂木亍Ｈ秉c(diǎn)是：必須對(duì)數(shù)據(jù)進(jìn)行整理以計(jì)算速度的2階矩(但實(shí)際上恒星的速度數(shù)據(jù)往往不是高斯分布，因此不能完全為二階矩所刻畫(huà))；分布函數(shù)的形式不能使用(因此未能使用恒星完全的分布信息)；Jeans方程組一般情況下不能閉合(即不能求解)，在某些情況下可能得到一個(gè)實(shí)際分布函數(shù)不存在的解。

如果這時(shí)候我們進(jìn)一步假設(shè)該星系中的恒星運(yùn)動(dòng)只遵守兩個(gè)運(yùn)動(dòng)積分(即文獻(xiàn)中常表述的f(E,L z))，那么σR=σz，而且σRz=0。此時(shí)，z方向的Jeans方程，就僅與各量的z方向的梯度有關(guān)，與R方向的運(yùn)動(dòng)解除耦合了，即：

在給定ν,Φ之后，我們可以得到σR(或σz)的表達(dá)式以及的表達(dá)式，也就是說(shuō)Jeans方程閉合了[8]。

在更早之前，就連上述2個(gè)運(yùn)動(dòng)積分情形下的Jeans方程都難以得到銀河系觀測(cè)數(shù)據(jù)(恒星的密度分布數(shù)據(jù)以及運(yùn)動(dòng)學(xué)數(shù)據(jù))的支持，研究者把引力勢(shì)的表達(dá)式(即泊松方程)也“一維化”了，簡(jiǎn)化為：

式(92)和(93)構(gòu)成了所謂的“一維近似”[61]。由于??Φ/?z即垂向力K z(引力的z方向分量)，早期文獻(xiàn)中又把這種一維近似稱(chēng)作“K zforce”方法[39]。如果用更方便的觀測(cè)量恒星面密度(Σ)來(lái)代替密度ρ，那么泊松方程可以表示為K z=2πGΣ。

在下面的幾個(gè)小節(jié)中，本文介紹近年來(lái)有代表性的建模工作。這幾個(gè)工作都是基于3個(gè)運(yùn)動(dòng)積分(或作用量)的分布函數(shù)來(lái)為銀河系建模。這些工作中，對(duì)于觀測(cè)到的恒星盤(pán)成分，一般是使用基于作用量的準(zhǔn)等溫球分布函數(shù)模型，或者基于運(yùn)動(dòng)積分的準(zhǔn)等溫球分布函數(shù)。對(duì)于暗物質(zhì)暈，這些工作使用質(zhì)量分布模型(如NFW模型)而不是分布函數(shù)模型。后來(lái)，Piffl等人[51]及Binney和Piffl[49]開(kāi)始嘗試使用分布函數(shù)模型來(lái)直接為暗物質(zhì)暈建模并擬合銀河系數(shù)據(jù)；這項(xiàng)暗物質(zhì)暈建模工作由于數(shù)據(jù)有限以及軟件上的問(wèn)題，目前的效果不太理想[15]。在本節(jié)的末尾，我們還介紹前人提出的針對(duì)恒星盤(pán)的另一種分布函數(shù)模型。

前人的這些工作，主要是基于RAVE和SDSS的銀河系巡天數(shù)據(jù)?，F(xiàn)在，在Gaia的天測(cè)和測(cè)光數(shù)據(jù)以及大規(guī)模光譜巡天(例如中國(guó)的LAMOST)的基礎(chǔ)上，我們預(yù)計(jì)這些基于分布函數(shù)建模的方法將大放異彩。

3.1 基于f(J R,L z,J z)的銀河系動(dòng)力學(xué)建模

3.1.1 Bovy和Rix的模型

Bovy和Rix[44]用SEGUE的G矮星樣本(范圍為5

(1)恒星的分布函數(shù)模型

準(zhǔn)等溫分布模型(quasi-distribution function,qDF)[45]是一種最常見(jiàn)的恒星分布模型。將示蹤星的分布函數(shù)建模為一個(gè)或多個(gè)等溫分量；根據(jù)定義，等溫分量具有恒定的速度色散，與z無(wú)關(guān)，從而速度分布必須是高斯分布[39]。Bovy和Rix[44]文中的每個(gè)MAP都是由準(zhǔn)等溫分布擬合的恒星種群；分布函數(shù)的自由參數(shù)，與示蹤劑群體的徑向密度分布，徑向和垂向速度色散，以及它們?nèi)绾坞SR變化有關(guān)。Bovy和Rix[44]所用的基于作用量的恒星分布模型來(lái)自文獻(xiàn)[46]：

(2)銀河系引力勢(shì)模型

為計(jì)算軌道作用量，Bovy和Rix[44]提出一個(gè)完整的三維模型來(lái)計(jì)算銀河系的勢(shì)：一個(gè)先驗(yàn)有許多自由參數(shù)的“核球+恒星盤(pán)+氣體盤(pán)+暗物質(zhì)暈”模型。

指數(shù)盤(pán)密度為雙冪律盤(pán)模型：

引力勢(shì)通過(guò)該密度的泊松方程的解析解得到。

核球模型為：

其中，α=1,β=4時(shí)，稱(chēng)為Hernquist模型；α=2,β=4時(shí)，稱(chēng)為Jaffe模型。暗物質(zhì)暈為冪律模型：

它們的共同引力勢(shì)為總引力勢(shì)。

在Bovy和Rix[44]對(duì)單個(gè)MAP的qDF擬合中，除了恒星盤(pán)的標(biāo)長(zhǎng)R d和暈對(duì)盤(pán)在R0處的徑向力的貢獻(xiàn)以外，其他基本參數(shù)(恒星盤(pán)的標(biāo)高z h，以及R0附近圓周速度v c(R0)和旋轉(zhuǎn)曲線的“平坦度”d lgv C(R0)/d lg(R))都被固定。

(3)擬合手續(xù)

Bovy和Rix[44]所使用的擬合方法為最大似然法，其擬合步驟如下，流程示意圖見(jiàn)圖3。

圖3 Bovy和Rix擬合流程的示意圖[44]

觀測(cè)量下的分布為：

其中，D為觀測(cè)量距離，l為銀經(jīng)，b為銀緯，g?r為顏色，[Fe/H]為觀測(cè)量金屬豐度，v為速度，由觀測(cè)得到，(R,z,Φ)為銀心柱坐標(biāo)系，|J|為坐標(biāo)變換的雅克比，S為選擇函數(shù)。

歸一化后，每個(gè)數(shù)據(jù)點(diǎn)的似然為：

加上參數(shù)條件是表示設(shè)參數(shù)pΦ,pDF取某一組值的條件，最大似然就是在參數(shù)取值空間中找到似然最大時(shí)的參數(shù)值。作用量的計(jì)算會(huì)用到總引力勢(shì)Φ(R,z)，這些引力勢(shì)的參數(shù)實(shí)際是待定的基準(zhǔn)勢(shì)；正是在給定分布函數(shù)模型和引力勢(shì)模型后，用最大似然法依靠數(shù)據(jù)從參數(shù)空間找到最合適的參數(shù)后，才能得到最合適的模型。因此，引力勢(shì)的參數(shù)，正是通過(guò)數(shù)據(jù)對(duì)作用量的依賴，輸入到觀測(cè)到的密度當(dāng)中[44]。

數(shù)據(jù)點(diǎn)是獨(dú)立的，考慮異常數(shù)據(jù)點(diǎn)(outliers)后，總的似然為：

最后根據(jù)最大似然法估計(jì)，使這個(gè)總似然最大，即找到所求參數(shù)pΦ,pDF的最合適值。

(4)結(jié)果

利用這些模型的所有參數(shù)的最佳結(jié)果，Bovy和Rix[44]得到了星系盤(pán)的徑向標(biāo)長(zhǎng)R d=(2.15±0.14)kpc，與通過(guò)光度法推斷的恒星質(zhì)量的空間分布一致。

來(lái)自不同的MAP數(shù)據(jù)，以不同的半徑約束Σ1.1(R)；因而B(niǎo)ovy和Rix[44]用這43個(gè)MAP數(shù)據(jù)的平均半徑和Σ1.1的平均值的關(guān)系，研究了Σ1.1隨半徑的變化，這些MAP的結(jié)果也一致[44]。類(lèi)似地，還得到了1.1kpc的垂向力K z,1.1(R)。使用從文獻(xiàn)[41]距離標(biāo)度的后驗(yàn)分布函數(shù)所導(dǎo)出的最佳半徑，將這些分別擬合為：

Bovy和Rix[44]在計(jì)算恒星盤(pán)對(duì)圓周速度的貢獻(xiàn)、暗物質(zhì)暈的冪律指數(shù)等結(jié)果方面，對(duì)完整Jeans方程和垂向力兩種方法進(jìn)行了對(duì)比(見(jiàn)文獻(xiàn)[44]中的圖18,19,20)；還展示了恒星盤(pán)和暗物質(zhì)暈分別對(duì)旋轉(zhuǎn)速度的局部貢獻(xiàn)(見(jiàn)文獻(xiàn)[44]中的圖21)；得到的對(duì)暗物質(zhì)暈的徑向輪廓的限制為：在95%置信度下，ρDM(r;UR0)∝1/rα,α<1.53。

3.1.2 Piffl等人的模型

與Bovy和Rix[44]提出的模型類(lèi)似，同樣用數(shù)據(jù)限制銀河系的參數(shù)，Piffl等人[50]用距銀河系平面約1.5kpc內(nèi)約200000顆巨星的運(yùn)動(dòng)學(xué)，來(lái)測(cè)量太陽(yáng)附近質(zhì)量密度的垂向分布，以及局部暗物質(zhì)暈貢獻(xiàn)。其所用的模型的組分更多，由氣體盤(pán)、薄盤(pán)、厚盤(pán)、核球和暗物質(zhì)暈五個(gè)部分組成，下面介紹其質(zhì)量模型和分布函數(shù)。

(1)質(zhì)量模型

盤(pán)的密度模型為：

其中，氣體盤(pán)非零參數(shù)Rhole是為了在圓盤(pán)中創(chuàng)建中心腔，而其他兩個(gè)盤(pán)則設(shè)置為零。

核球和暗物質(zhì)暈的密度模型為：

因此Piffl等人[50]的質(zhì)量模型具有8個(gè)參數(shù)，薄盤(pán)和厚盤(pán)分別具有3個(gè)自由參數(shù)(Σ0,R d,z d)，而暗物質(zhì)暈具有2個(gè)自由參數(shù)ρ0,r0，其中薄盤(pán)和厚盤(pán)的徑向標(biāo)長(zhǎng)取相同。文獻(xiàn)[50]中的表1給出了銀河系模型的固定參數(shù)取值。

(2)分布函數(shù)模型

分布函數(shù)分為盤(pán)和暈兩大部分，

其中，fdisk為盤(pán)的分布函數(shù)，fhalo為暈的分布函數(shù)。Fhalo是暈與盤(pán)的質(zhì)量之比。

盤(pán)的分布函數(shù)為：

其中，fthin為薄盤(pán)的分布函數(shù)，fthick為厚盤(pán)的分布函數(shù)，F(xiàn)thick是厚盤(pán)與薄盤(pán)的質(zhì)量之比。

其中，ν仍代表垂向頻率而不是橢球坐標(biāo)。但是Piffl等人[50]的模型中，薄盤(pán)還考慮了星族年齡，其速度彌散除了L z，還依賴于年齡τ，

進(jìn)一步假設(shè)薄盤(pán)中的恒星形成率隨時(shí)間指數(shù)下降，特征時(shí)間標(biāo)度為t0，于是薄盤(pán)的分布函數(shù)為：

Piffl等人[50]還在分布函數(shù)中添加了用于恒星暈的組分fhalo，以防止擬合過(guò)程使厚恒星盤(pán)變形，從而解釋樣本中的暈星的存在。恒星暈的分布函數(shù)為Posti等人[27]的冪律模型

Piffl等人[50]采用最小二乘法來(lái)建立這些模型，先給出了最佳擬合的待定參數(shù)，然后給出了各組分的密度輪廓，他們還分別給出了重子物質(zhì)與暗物質(zhì)暈對(duì)圓周速度和中盤(pán)面總面密度隨徑向的局部貢獻(xiàn)變化。

3.2 基于f(E,L z,I3)的銀河系動(dòng)力學(xué)建模

前面幾個(gè)準(zhǔn)等溫模型及其變形都是基于作用量的分布函數(shù)模型。除了基于作用量的，以前也有直接基于E,L z,I3三個(gè)運(yùn)動(dòng)積分的分布函數(shù)模型。關(guān)于第三運(yùn)動(dòng)積分，大體上也是以Stckel引力勢(shì)中的運(yùn)動(dòng)積分的形式為藍(lán)本，作近似假設(shè)。下面簡(jiǎn)介其中的一些模型。

這里，(r,z,Φ)為柱坐標(biāo)，為z軸上的焦距。

恒星盤(pán)分布函數(shù)假定為一個(gè)三運(yùn)動(dòng)積分模型：

并用Jeans方程進(jìn)行了檢驗(yàn)。

3.2.2 一個(gè)更多參數(shù)的分布函數(shù)模型

下面介紹一種參數(shù)更多、復(fù)雜的分布函數(shù)模型[40,42]。Famaey等人[42]構(gòu)建基于3個(gè)運(yùn)動(dòng)積分的分布函數(shù)f(E,L z,I3)，從而實(shí)際星系的恒星盤(pán)通過(guò)這種形式的分布函數(shù)(或者若干此種分布函數(shù)的線性組合)來(lái)正確地?cái)M合出。他們的方法是，(1)對(duì)于I3的計(jì)算，簡(jiǎn)單套用Stckel勢(shì)的I3表達(dá)式；(2)對(duì)于分布函數(shù)的形式，在前人的一個(gè)基于2個(gè)運(yùn)動(dòng)積分的分布函數(shù)[40]的基礎(chǔ)上，引入更多自由參數(shù)作拓展。

分布函數(shù)的具體形式為：

通過(guò)計(jì)算分布函數(shù)下若干階速度的平均，即可得到該分布函數(shù)(運(yùn)動(dòng)積分空間)的矩：

從而0階矩為質(zhì)量密度如ρ(λ,ν)=μ0,0,0(λ,ν)，二階矩為速度彌散如(λ,ν)=μ0,0,2(λ,ν)。文獻(xiàn)[42]的第4章給出了一種方法，通過(guò)由積分區(qū)域邊界的線性組合表示積分區(qū)域中的點(diǎn)，來(lái)計(jì)算其各階矩中的積分；本文附錄B簡(jiǎn)介了其基本思想。有了矩的計(jì)算，文獻(xiàn)[42]的第5章介紹并分析了分布函數(shù)中各種參數(shù)所體現(xiàn)的物理特征。

文獻(xiàn)[42]最后介紹了對(duì)實(shí)際的恒星盤(pán)建模的方法：通過(guò)將分布函數(shù)組分作為基礎(chǔ)函數(shù)線性組合。這需要在配置空間中引入一個(gè)網(wǎng)格，并最小化方差(如對(duì)于分量FΛ代表0階矩，即密度時(shí))：

其中，Λ=(α1,α2,β,δ,η,z0,s,a)為分布函數(shù)的參數(shù)集，cΛ為系數(shù)。擬合程序中不斷加入滿足要求特征的新組分，直至χ2的最小值不顯著變化。

4 討論和總結(jié)

前面內(nèi)容聚焦于準(zhǔn)穩(wěn)態(tài)時(shí)星系的分布函數(shù)特點(diǎn)及建模。但是，與常見(jiàn)的穩(wěn)態(tài)、準(zhǔn)穩(wěn)態(tài)的熱力學(xué)系統(tǒng)相比，自引力系統(tǒng)有些特點(diǎn)本質(zhì)不同，這些特點(diǎn)也許會(huì)限制前述方法的適用范圍。另一方面，當(dāng)考慮處于非穩(wěn)態(tài)(遠(yuǎn)離平衡時(shí))的星系結(jié)構(gòu)和演化過(guò)程，或者考慮少體引力系統(tǒng)(粒子數(shù)較少的引力系統(tǒng))中的具體恒星的運(yùn)動(dòng)，上述基于作用量的建模方法是否有效，或者說(shuō)如何使得上述方法仍能應(yīng)用于某些具體場(chǎng)合，這是非常值得探索的課題。

在本章里，我們將首先在動(dòng)力系統(tǒng)理論的視角下(dynamical-systems theory)，研究自引力系統(tǒng)的動(dòng)力學(xué)(見(jiàn)4.1節(jié))；然后介紹和評(píng)述近年來(lái)研究者把“作用量-角變量”方法應(yīng)用到完全不可積的星系現(xiàn)象和星系結(jié)構(gòu)的進(jìn)展(共振和星系棒，見(jiàn)4.2節(jié))，應(yīng)用于非平衡的星系過(guò)程——潮汐星流(tidal streams)的研究，包括星流的形成過(guò)程，在作用量空間、角變量空間對(duì)星流的刻畫(huà)，以及把星流用來(lái)擬合主星系的引力勢(shì)分布(見(jiàn)4.3節(jié))；最后是全文的總結(jié)(見(jiàn)4.4節(jié))。

4.1 動(dòng)力系統(tǒng)理論視角下的研究

從統(tǒng)計(jì)物理角度，引力N體系統(tǒng)是一種長(zhǎng)程相互作用系統(tǒng)。與一般的近程作用系統(tǒng)不同，這類(lèi)系統(tǒng)中任何粒子之間的作用幾乎都不能忽略(粒子間的作用勢(shì)ψ∝r?α，且α小于位形空間維數(shù)d；假設(shè)粒子在全空間均勻分布，那么總作用勢(shì)將是發(fā)散的)，并且與勢(shì)能有關(guān)的宏觀量不再具有統(tǒng)計(jì)物理所經(jīng)常要求的可加性。另一方面，自引力系統(tǒng)由于可形成核-暈結(jié)構(gòu)使其熵可能不具有最大值，因此并沒(méi)有一個(gè)真正的平衡態(tài)；一般來(lái)說(shuō)，相比于近程作用系統(tǒng)，自引力系統(tǒng)的弛豫過(guò)程可以更快，例如劇烈弛豫機(jī)制(violent relaxation)[55,58]。這些特點(diǎn)導(dǎo)致自引力系統(tǒng)與統(tǒng)計(jì)物理中常見(jiàn)的近程作用系統(tǒng)顯著不同，往往不能用平衡態(tài)統(tǒng)計(jì)方法來(lái)很好地研究。

從動(dòng)力學(xué)角度，作用量-角變量坐標(biāo)將n自由度可積系統(tǒng)的2n維∑相空間變成n維環(huán)面T n，如果存在一組非零整數(shù)k1,k2,k3,···,k n，使得角頻率?滿足時(shí)，則存在共振。共振關(guān)系可有不同多種。龐加萊(Poincar)截面法就是在相空間(即坐標(biāo)維+速度維)中取一個(gè)低于相空間維度的“面”，記錄系統(tǒng)的軌道穿過(guò)這個(gè)面的點(diǎn)，并研究這些落點(diǎn)規(guī)律，該方法是降維研究運(yùn)動(dòng)的常見(jiàn)且直觀的方法。以3自由度可積系統(tǒng)為例(如在Stckel勢(shì)形式的星系中的任意一個(gè)恒星軌道)，該系統(tǒng)在相空間中呈現(xiàn)為一個(gè)三維環(huán)體(3-torus)；若系統(tǒng)的其中兩個(gè)頻率比?1/?2為有理數(shù)，則在其相應(yīng)的二維環(huán)面的運(yùn)動(dòng)是周期的，其在環(huán)面上的運(yùn)動(dòng)軌跡為周期軌線，在龐加萊截面上則是一些離散的不動(dòng)點(diǎn)；若運(yùn)動(dòng)是準(zhǔn)周期的，則是布滿環(huán)面的準(zhǔn)周期軌線，落點(diǎn)在截面上將形成一條閉合曲線。

當(dāng)可積系統(tǒng)受不可積微擾且成為近可積系統(tǒng)(近可積系統(tǒng)是大部分的相空間仍被規(guī)則軌道占據(jù)的受擾可積系統(tǒng))，KAM定理給出了近可積系統(tǒng)中仍然穩(wěn)定存在局部的不變環(huán)面的條件。KAM定理是20世紀(jì)中葉Kolmogorov,Arnold和Moser提出的重要定理，后續(xù)的工作放寬了定理的條件。對(duì)于截面上的不動(dòng)點(diǎn)，當(dāng)施加擾動(dòng)，簡(jiǎn)單來(lái)說(shuō)，經(jīng)過(guò)多次運(yùn)動(dòng)后，雙曲(不穩(wěn)定)不動(dòng)點(diǎn)將逐漸偏離，橢圓(穩(wěn)定)不動(dòng)點(diǎn)則將圍繞不動(dòng)點(diǎn)形成閉合環(huán)面，或在滿足KAM定理的一定區(qū)域內(nèi)(即KAM環(huán)面)穩(wěn)定存在，或形成更高階的新的兩種不動(dòng)點(diǎn)。如此迭代，將導(dǎo)致更復(fù)雜和不可預(yù)知的演化，從而產(chǎn)生混沌。此外，對(duì)于n>2自由度的近可積系統(tǒng)，相空間為2n維，等能面(或稱(chēng)能量流形)為(2n?1)維，而n維的KAM環(huán)面不再能將等能面分割成閉區(qū)域的集合，即存在阿諾德擴(kuò)散現(xiàn)象(Arnold diffusion)。阿諾德擴(kuò)散使得n≥3自由度的動(dòng)力系統(tǒng)的行為更為復(fù)雜，在穩(wěn)定的軌道之間存在著高度不穩(wěn)定的軌道，這樣運(yùn)動(dòng)軌跡就可以繞過(guò)KAM環(huán)面到其他區(qū)域，仍會(huì)出現(xiàn)混沌；人們對(duì)于阿諾德擴(kuò)散的理解，目前仍然非常有限。大致上，我們可以說(shuō)，有理環(huán)面的破壞源于運(yùn)動(dòng)模式的非線性共振，非線性共振進(jìn)一步導(dǎo)致混沌運(yùn)動(dòng)[63]?，F(xiàn)代的“動(dòng)力系統(tǒng)”學(xué)科本身也在拓?fù)浜头治錾嫌辛烁鼑?yán)密的定義和更多的理論，我們期待著這些現(xiàn)代理論能在自引力系統(tǒng)的動(dòng)力學(xué)研究中大放異彩。

4.2 應(yīng)用于“不可積”的星系現(xiàn)象

絕大多數(shù)自引力系統(tǒng)(如星系)都不是可積系統(tǒng)；即使數(shù)值計(jì)算，由于混沌的存在，也幾乎不可能精確計(jì)算某個(gè)具體天體(如某個(gè)恒星)實(shí)際的長(zhǎng)期運(yùn)動(dòng)。盡管如此，對(duì)于處于準(zhǔn)穩(wěn)態(tài)的恒星集體(即星系)而言，前述章節(jié)所介紹的Binney等人的工作表明：恒星集體的統(tǒng)計(jì)性質(zhì)，例如各種分布，是可以利用“作用量-角變量”方法來(lái)可靠地獲得的。我們進(jìn)一步感興趣的是，“作用量-角變量”方法是否可以應(yīng)用于共振、混沌等所謂“不可積的”的恒星運(yùn)動(dòng)，可以應(yīng)用于哪些具體場(chǎng)合？最近Binney[16,17]研究了把環(huán)面映射方法應(yīng)用到棒旋星系的共振軌道的情況(以銀河系觀測(cè)數(shù)據(jù)為實(shí)例)。他在文獻(xiàn)[16]中，分別研究了在銀河系棒的這種擾動(dòng)強(qiáng)度下，“作用量-角變量”工具對(duì)于外林德布拉德共振、共轉(zhuǎn)共振、內(nèi)林德布拉德共振的適用性。他發(fā)現(xiàn)，對(duì)于銀河系被外林德布拉德共振和共轉(zhuǎn)共振所俘獲的恒星軌道，環(huán)面映射法能很好地?cái)M合出；但是更靠近棒內(nèi)部的內(nèi)林德布拉德共振俘獲的軌道，環(huán)面映射法就顯得不夠好了。Binney[17]更進(jìn)一步把上述方法應(yīng)用到太陽(yáng)鄰域恒星的速度空間數(shù)據(jù)中，發(fā)現(xiàn)在銀河系的棒模型下，環(huán)面映射法所擬合的共轉(zhuǎn)共振情形與觀測(cè)數(shù)據(jù)一致，而外林德布拉德共振情形不被數(shù)據(jù)所支持；也就是說(shuō)，基于環(huán)面映射法和微擾理論的計(jì)算，對(duì)于棒旋星系中的共振俘獲現(xiàn)象，可以有很好的實(shí)用性了。

4.3 應(yīng)用于遠(yuǎn)離平衡的星系過(guò)程：潮汐星流

當(dāng)一個(gè)星團(tuán)或矮星系(稱(chēng)為前身天體)靠近一個(gè)大星系(稱(chēng)為主星系)時(shí)，前身天體各個(gè)不同位置所受的引力不同，從而發(fā)生形狀和速度的變化；當(dāng)所處的引力場(chǎng)梯度足夠大時(shí)，前身天體的恒星們會(huì)被主星系巨大潮汐力逐步拉扯，或多或少的恒星脫離前身天體，常形成長(zhǎng)條狀或稱(chēng)纖維狀的結(jié)構(gòu)，此即潮汐星流(詳見(jiàn)文獻(xiàn)[33])。

研究者發(fā)現(xiàn)，對(duì)于星流的這種形成過(guò)程，以及對(duì)于星流特性的刻畫(huà)，最方便的方式是基于“作用量-角變量”體系[34,35]。當(dāng)前身矮星系或星團(tuán)的成員恒星沒(méi)有被潮汐剝離，隨著前身天體在主星系引力勢(shì)中運(yùn)動(dòng)時(shí)，在主星系參考系中，這些恒星的作用量是不守恒的。但一旦某恒星被潮汐力剝離前身天體之后，就可以忽略前身天體以及星流中的其他恒星對(duì)它的引力作用，認(rèn)為只是在主星系的引力勢(shì)下運(yùn)動(dòng)。那么，自剝離時(shí)刻起，該恒星的作用量就一直保持剝離前的數(shù)值，不再隨時(shí)間變化了。下面我們沿用Eyre和Binney[36]的記號(hào)和表述方式，來(lái)簡(jiǎn)介這個(gè)方法。設(shè)前身天體(小星系或星團(tuán))位于固定背景勢(shì)下的規(guī)則軌道上，其作用量為J0。背景勢(shì)具有哈密頓量H(J)，其中J是潮汐星流中的某恒星的作用量，一般與J0相差不遠(yuǎn)。對(duì)于潮汐星流中的某個(gè)恒星，其作用量是常量，角變量隨時(shí)間線性增加。在J0附近，對(duì)哈密頓量作泰勒展開(kāi)：

其中δJ=J?J0；?0是前身天體的軌道的頻率，是哈密頓量的海森矩陣(即H的二階微分)，則位于前身天體軌道附近的某恒星的軌道頻率為：

關(guān)于星流的成員恒星，它們的作用量和角變量都是有著某種分布，即某種程度的彌散，?J和?θ0。由前述可知，成員恒星的作用量彌散?J自從星流形成后就不變了；但角變量彌散不會(huì)保持在初始的?θ0，而是將隨時(shí)間線性增大。角變量彌散的關(guān)系式為[34,36]：

其中約等號(hào)是由于潮汐流形成后，θ0將遠(yuǎn)小于正比于時(shí)間的那一項(xiàng)。于是有

上式表明海森矩陣D的作用是：把恒星們?cè)谧饔昧靠臻g的分布(參見(jiàn)文獻(xiàn)[36]的圖10)線性映射到角變量空間的分布(參見(jiàn)文獻(xiàn)[36]的圖13)。D是一個(gè)對(duì)稱(chēng)矩陣，因此其特征取決于其相互正交的特征向量(n=1,2,3)，和實(shí)特征值λn[36]。海森矩陣決定了星流的拉伸方向，即星流沿著最大特征方向拉長(zhǎng)[20]：

什么樣的引力勢(shì)分布形式(必要條件)可以形成星流等技術(shù)細(xì)節(jié)，這里不作詳細(xì)介紹，讀者可以參考Sanders的博士論文(第5章和第6章)[20]。特別是該論文第6章詳細(xì)介紹了如何應(yīng)用星流數(shù)據(jù)來(lái)擬合限定主星系的引力勢(shì)分布。

最后，我們要指出的是：當(dāng)前身天體為很小的星系或者星團(tuán)時(shí)，星流的軌跡可以粗略地(1階近似)看作是前身天體在主星系中的軌道；因此，過(guò)去研究者常簡(jiǎn)單假定星流為一條軌道，從而擬合主星系引力勢(shì)，此即文獻(xiàn)中所說(shuō)的軌道擬合法(例如文獻(xiàn)[38])。但是，基于上面的分析，星流的軌跡并不必然與前身天體的軌道或者任何其他軌道重合(參見(jiàn)文獻(xiàn)[36]的圖13)；因此軌道擬合法會(huì)導(dǎo)致較大的主星系引力勢(shì)偏差(詳見(jiàn)文獻(xiàn)[20]的第6章)。

4.4 結(jié)束語(yǔ)

本文主要目的是介紹星系動(dòng)力學(xué)建模方面的新進(jìn)展，即作用量計(jì)算方法，以及近年來(lái)基于作用量(或運(yùn)動(dòng)積分)的分布函數(shù)建模在銀河系數(shù)據(jù)上的應(yīng)用，另外還簡(jiǎn)介了“作用量-角變量”方法在不可積的星系現(xiàn)象(共振和星系棒)和非平衡過(guò)程(例如潮汐星流的形成)的應(yīng)用，以期為中國(guó)在“Gaia+LAMOST”黃金時(shí)代的銀河系研究做些貢獻(xiàn)。目的之二，厘清上述方法及其應(yīng)用的理論基礎(chǔ)和脈絡(luò)。在實(shí)際或理想的星系中，唯一的一類(lèi)可積系統(tǒng)是具備Stckel形式引力勢(shì)的星系；相應(yīng)地，幾乎所有星系建模工作的解析性基礎(chǔ)，都是Stckel勢(shì)理論。因此，我們還在本文第1章花不少篇幅從N個(gè)粒子的經(jīng)典力學(xué)(哈密頓方程和劉維爾方程)出發(fā)，依次引入平均場(chǎng)近似，介紹無(wú)碰撞玻爾茲曼方程；然后在經(jīng)典力學(xué)的正則變換以及相空間的語(yǔ)境下，引入運(yùn)動(dòng)積分、正則環(huán)面坐標(biāo)即作用量和角變量；接著就介紹唯一的嚴(yán)格可積的(理想的)“星系”，即介紹Stckel勢(shì)理論的各種解析式。介紹這些來(lái)龍去脈，期望我國(guó)學(xué)者不僅能應(yīng)用Binney等國(guó)外專(zhuān)家發(fā)明發(fā)展的方法來(lái)分析銀河系數(shù)據(jù)，還能夠在理論基礎(chǔ)和方法上同樣有所創(chuàng)新。

附錄A 計(jì)算作用量、角頻率和角變量的步驟

這里介紹Sanders和Binney的程序中數(shù)值計(jì)算作用量、角頻率和角變量的步驟(文獻(xiàn)[1]中第8節(jié)及附錄；文獻(xiàn)[22]的附錄A)。程序中的這些中間量的計(jì)算以及軌道分類(lèi)，來(lái)自Sanders和Binney的TACT程序[23]的/tact/aa/inc/stackel aa.h中的函數(shù)定義。

(1)作用量J

作用量的計(jì)算公式為

其中，τ?,τ+是=0時(shí)的根。注意對(duì)于有些環(huán)軌道，τ的積分上下限可能是?α,?β,?γ，具體可見(jiàn)文獻(xiàn)[22]的表1，其詳細(xì)討論可見(jiàn)文獻(xiàn)[1]中的軌道分類(lèi)內(nèi)容。

(2)角頻率?

由哈密頓方程(?τ==?H/?Jτ)以及微分關(guān)系，可以得到[1]：

將這個(gè)關(guān)于角頻率的線性方程反解可得角頻率，如

可用式(125)與動(dòng)量和運(yùn)動(dòng)積分的關(guān)系計(jì)算：

注意上式右端作為分母的pτ，在積分邊界會(huì)消失。Sanders和Binney[22]作了下述變換：

使被積函數(shù)在積分邊界更平滑地變到零。這樣，每一點(diǎn)積分變量的值都得到了，再數(shù)值積分得到角頻率。程序中用的是高斯-勒讓德正交格式來(lái)進(jìn)行積分。

(3)角變量θ

用哈密頓特征函數(shù)計(jì)算角變量，表示為：

計(jì)算?W/?I所用的W(τ,J)為：

其中，被積函數(shù)的處理類(lèi)似于前述所示；函數(shù)Fτ(pτ,x)是考慮軌道類(lèi)型的特點(diǎn)，為消除橢球坐標(biāo)中的簡(jiǎn)并性以使角變量覆蓋笛卡爾坐標(biāo)中一個(gè)振動(dòng)的全部范圍而引入的(文獻(xiàn)[22]附錄中的式(A6))，即計(jì)算中還要進(jìn)行軌道方面的判斷和修正。

附錄B Famaey等人[42]的運(yùn)動(dòng)積分分布函數(shù)模型的矩的計(jì)算處理方法簡(jiǎn)介

對(duì)于式(118)，從被積的速度空間轉(zhuǎn)換到3個(gè)運(yùn)動(dòng)積分所構(gòu)成的空間

通過(guò)對(duì)積分邊界線(上式取等號(hào)時(shí))的線性組合，來(lái)表示式(118)中的積分區(qū)域中的點(diǎn)，

之后，將(E,I3)的積分區(qū)間轉(zhuǎn)換到新的積分變量空間(固定L z)，然后對(duì)比變換前后的每一項(xiàng)，解出上式中的組合系數(shù)v,u,t的關(guān)系，從而用一維數(shù)值積分計(jì)算其矩(更多細(xì)節(jié)見(jiàn)文獻(xiàn)[42]中的4.2和4.3節(jié))。