Stckel引力勢理論的拓展及其應用

2021-07-13 01:34:28顧建宇董小波

天文學進展 2021年2期

關鍵詞：方法

顧建宇，董小波

(1.中國科學院云南天文臺，昆明650216；2.中國科學院大學，北京100049)

1 引言和理論基礎

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

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

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

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

本文討論的單粒子在整體勢的運動，是從引力多體運動簡化而來的，即單粒子近似或稱平均場近似[60]。三維背景空間中一個包含有N個粒子的經典系統(tǒng)，可用哈密頓方程描述：

對于該系統(tǒng)所有粒子的坐標和動量所構成的全6N維相空間(一般稱為Γ相空間)，通過相空間代表點個數的局部守恒可推導得，其概率密度ρ(q,p)滿足劉維爾方程：

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

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

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

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

此式實質上等同于流體力學動量守恒方程：

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

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

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

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

1.2 運動積分和作用量-角變量體系簡介

1.2.1 哈密頓-雅克比方程

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

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

1.2.2 運動積分

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

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

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

此時的新哈密頓量仍是K=H，并且新哈密頓函數只含新廣義動量而不含新廣義坐標，該正則變換是作為后面討論運動積分坐標的基礎。

更進一步，對于自治(完全)可分離的生成函數W，可寫為：

則相應地可以分離出n個哈雅方程(不對i求和)，可表示為：

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

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

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

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

1.2.3 作用量-角變量體系

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

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

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

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

于是我們可以將相空間軌道展開成傅里葉級數：

角頻率之間的比例關系豐富了體系周期類型和遍歷性的復雜性(見第4章)。天體動力系統(tǒng)的周期特性一般來說非常強，這也是作用量-角變量體系在天體力學中經久不衰的重要原因。

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

其中，T是正則變換的矩陣，指在γi對應的相空間所在區(qū)域。也由此可見，某個作用量的值J i代表子環(huán)γi所圍的面積(當對q i,p i選取合適的標度)，而θi代表轉過的角度。

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

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

1.2.4 基于守恒量的分布函數

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

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

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

1.3 Stckel勢簡介

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

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

圖1 共焦橢球坐標系的等坐標面示意圖

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

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

由于式(21)的特殊形式以及W的可分離性，可寫為于是式(23)中括號內的部分只是各自τ坐標的函數，設之為Uτ(τ)并替換，可得：

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

由于Uτ(τ)即式(23)的中括號內的部分，則可得：

更方便地，H,J,K可寫為：

其中，

此外，還有更常用的另一套運動積分，兩者的關系是：

此即星系動力學中常說的“第二個、第三個運動積分”I2,I3。相關推導見文獻[2]。類似地，由pτ=?W/?τ，得到動量的表達式為

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

2 基于Stckel勢近似的作用量估算方法

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

2.1 Stckel作弊法

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

一般來說，軸對稱勢下有3個孤立運動積分，即總能量E，角動量垂向分量L z，以及I3。軸對稱橢球坐標系相對于三軸橢球坐標系而言，一個半軸長參數與另外一個相等，如果α=β，則為長等軸(prolate)橢球坐標系(λ,?,ν)，其中?代表某個子午面的方位角；如果γ=β，則為扁等軸(oblate)橢球坐標系(λ,μ,χ)，其中χ代表某個子午面的方位角。因此在描述粒子在軸對稱勢中的運動時，可以消去一個自由度。

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

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

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

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

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

對于一個一般引力勢，定義三軸情形的輔助函數：

當固定另外兩個橢球坐標時，某個τ對應的Cτ,Dτ可看作常數，從而使我們可以分開來獨立計算所需的運動積分。

將式(35)代入到式(26)，得到動量的另一個計算式

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

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

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

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

用式(41)來計算式(40)第一項，可見同時可以消去?Φ/?λ，因此得到Jλ以及類似的Jμ,Jν的表達式：

于是得到6個估算的運動積分Jτ,Kτ的表達式，其中只有3個是獨立的[22]。

(2)合適焦距的選取

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

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

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

2.2 Stckel擬合法

這里介紹Sanders[19]提出的長等軸橢球坐標系下的一般軸對稱勢的Stckel擬合法。該方法是對一般引力勢進行Stckel形式的擬合，即采樣軌道來擬合出一個合適的Stckel引力勢的分離函數F(τ)，來得到運動積分。Sanders和Binney[24]把Stckel作弊法和這種Stckel擬合法歸類為“非收斂方法”，因為這兩種方法依賴于做一定的近似，而近似所導致的誤差不能通過增加計算量來不斷減小。

在某個子午面上，運動的位形空間從柱坐標(R,?,z)變換到扁球坐標(λ,?,ν)依賴于

其中，λ,ν是上式作為關于τ的方程根，a2=α,b2=β,c2=γ，類似于三軸情形下的式(19)。此時Stckel勢的形式簡化為：

軸對稱時第三運動積分的解析表達式為：

動量和運動積分的關系表示為：

(1)給定初始的(x,v)，以及真實的勢(軌道積分中所感受到的勢)，并在真實的勢中計算實際的運動積分E,L z。

(2)焦距的選取。作上述變換到橢球坐標。然而，變換的參數，即焦距?=a2?c2(因為a=b，只要一個a2?c2便可確定橢球坐標系，可取c2=1)，非常影響作用量的標準差，因此通過一定的手續(xù)來估算一個更好的焦距?。Sanders[19]采用的估計方法是，由于Stckel勢的分離性，可得

由上式以及兩種坐標之間的變換關系，可推得[19]：

既然假設所估算的局部勢Φ就是該Stckel的勢ΦS，于是就用該式來估計a2的值，其中的偏導用已知的真實Φ來數值計算梯度而得；每個a2值，需要在軌道中每個點的相應值來計算，最后做平均，這些點在下一步的軌道采樣中獲得。

(3)軌道采樣以及運動區(qū)域的估算。進行軌道積分，即在(x,v)坐標下通過每一點的勢和時間積分來算下一步的坐標(x′,v′)，即相當于數值法畫軌道。Sanders和Binney的TACT程序中用守恒格式的數值方法畫出軌道，并且每個軌道周期離散為300個點。

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

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

最小。歸一化的權重函數形式為：

這個權重函數是提前選的，相關內容可見文獻[1,4]。于是擬合中所選用的Stckel勢的具體形式為：

其中，

這是因為在這種情況下偏差是最小的。證明方法是，式(49)取最小時，對于λ,ν分別成立，即將其分別對λ,ν求導應為0，得

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

接下來，有了式(50)的形式后，便可結合式(51)和(52)直接積分得到某個點(λ,ν)的F(τ)，如此選取若干個坐標點(TACT程序選取了40個點)，用來三次樣條插值得到所擬合的勢(樣條插值使擬合的勢會更平滑)。

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

(6)計算作用量-角變量。計算的F(τ)除了上一步用來算I3，還用來在式(46)中以此得到當前點的pτ，從而計算作用量J和角變量θ(數值計算步驟具體見附錄A)。

早在1985年，de Zeeuw和Lynden-Bell[2]就已提出了把Stckel勢作為模型，來局部擬合和全局擬合星系引力勢的思想和方法。前者在具體星系引力勢的平衡點附近展開引力勢，然后通過對應Stckel的橢球坐標形勢來確定參數；后者通過以平均的方法來近似Stckel勢的3個分離函數，從而計算運動積分。下面簡單介紹。

(1)局部擬合

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

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

推導中用到了(λ+α)的展開式，比如：

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

(2)全局擬合

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

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

其中，Λ(λ),M(μ),N(ν)是根據積分的約束選擇的權重函數，相應的歸一化常數為：

定義輔助函數

那么，分離函數可計算為：

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

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

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

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

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

通過作用量-角變量的實數性以及哈密頓量的時間反演不變性，可以對生成函數的系數進行約束[11]：

于是有[21]：

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

這樣就通過生成函數聯(lián)系起了玩具環(huán)面和目標環(huán)面。

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

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

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

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

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

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

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

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

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

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

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

2.3.3 軌道積分擬合生成函數法

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

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

(1)軌道積分與軌道分類

進行軌道積分。在沒有轉動的情況下，一般三軸勢允許兩類基本的非共振軌道：環(huán)軌道和盒軌道[1,21]。同一真實勢，不同的初始條件可能導致不同類型的軌道，于是TACT為不同類型的軌道匹配相應的玩具勢，與給定軌道類的環(huán)面相同的幾何結構。

(2)玩具勢的擬合

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

對于盒軌道，玩具勢為諧振子勢：

其判定依據為軌道類型，即看L x和L z的方向變化。在用玩具引力勢之前，需要在軌道上采集樣點，擬合出相應玩具引力勢的質量和標長參數。

(3)生成函數

在每個相空間樣點t i計算相應的玩具作用量-角變量(θtoy(t i),Jtoy(t i))，并產生相應的式(66)，其中各J和S n未知。Sanders和Binney軌道積分擬合生成函數法的實現(xiàn)技巧為：“我們不能精確地解這些方程，因為處理的方程有無窮多個未知數。我們可以在每個方程的右邊只包含有限個項，則方程右邊和左邊不能精確一致，因此可以采取的正確方法是：對每個方程的左右邊的差作平方，然后對所有方程的差平方求和(即軌道積分)，并最小化這個總的平方和”[21]。其總平方和的表達式是：

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

下面最小化E1，即令的偏導為0，

這些方程可轉化為一個矩陣方程組

其中，

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

類似地，對于求角變量θ′，則最小化

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

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

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

2.4 TACT所提供的方法的效果對比

我們運行TACT程序[23,24]用幾種方法分別估算出作用量，如圖2所示。我們所使用的引力勢為某種軸對稱引力勢，即McMillan[57]擬合的最佳銀河系質量分布模型，分為薄厚兩個恒星盤、HI和H2兩個氣體盤以及核球和暗物質暈；各個星系成分的模型輸入參數(適合TACT程序的格式)，已在作者的github主頁①https://github.com/PaulMc Millan-Astro/GalPot/blob/master/pot/PJM17_best.Tp ot。給出。我們采用的恒星軌道為典型的銀河系薄盤軌道，初始條件為(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)來估算焦距；Fit對應Stckel勢擬合法，ItTC對應環(huán)面迭代法，O2GF對應軌道積分擬合生成函數法，AvGF為平均生成函數法；CAA和SAA分別為柱絕熱近似和球絕熱近似方法，下一節(jié)作簡介。軌道中的每個點都獨立計算一次作用量，因而作用量的數值隨時間有所波動。

圖2 作用量數值計算的效果對比

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

2.5 其他一些求運動積分的方法或程序

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

2.5.1 “近似的運動積分”解析式

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

同理有

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

這樣就能分別用一個“橢球坐標軸”(如(λ,?β,?γ)中后兩個坐標固定)上的勢值來表示任意點的引力勢。

采取類似的手法，這個引力勢下的運動積分可以解析性地表達出來。運動積分J,K的表達式為

相關推導見參考文獻[2]中的式(9)和(14)。將式(79)和(80)代入，可計算得到與式(78)同樣方式的運動積分J,K的表達式。

類似地，由式(15)，I2,I3(第一個運動積分是能量E)的表達式可寫為：

即參考文獻[54]中的式(9)和(10)。文獻[54]推導了上式中的橢球坐標的函數項，得到

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

該方法計算“每一點的”運動積分，需要當前點的坐標、速度和角動量等軌道信息。文獻[54]在后文用高度精確的數值方法計算出軌道，并研究了這種方法所計算的運動積分的誤差，最后做了Jeans方程的檢驗。

2.5.2 絕熱近似

對于軸對稱勢，文獻[24]還介紹了用近似公式來估算作用量的方法。通過假設在柱坐標系或球坐標系中引力勢是可分離的，可求出這些作用量。

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

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

其中，R c是導心半徑(可參見有關“本輪近似”的文獻)。

(2)球絕熱近似假設沿長球坐標的球面存在絕熱不變振蕩。作用量的近似公式為：

其中，

具體計算可見參考文獻[24]。

2.5.3 一些相關程序

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

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

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

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

2.5.4 深度學習在作用量計算中的應用

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

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

3 應用

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

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

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

為2階矩。在觀測上，速度彌散可由某一個位置處諸恒星的速度方差來計算。

正如Read所總結[61]，Jeans方程和矩方法的優(yōu)點是：與其他方法相比非?？?，可以探索較大的參數空間；并且不需要假設分布函數的形式，因為只是限制了它的矩。缺點是：必須對數據進行整理以計算速度的2階矩(但實際上恒星的速度數據往往不是高斯分布，因此不能完全為二階矩所刻畫)；分布函數的形式不能使用(因此未能使用恒星完全的分布信息)；Jeans方程組一般情況下不能閉合(即不能求解)，在某些情況下可能得到一個實際分布函數不存在的解。

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

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

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

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

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

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

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

3.1.1 Bovy和Rix的模型

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

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

準等溫分布模型(quasi-distribution function,qDF)[45]是一種最常見的恒星分布模型。將示蹤星的分布函數建模為一個或多個等溫分量；根據定義，等溫分量具有恒定的速度色散，與z無關，從而速度分布必須是高斯分布[39]。Bovy和Rix[44]文中的每個MAP都是由準等溫分布擬合的恒星種群；分布函數的自由參數，與示蹤劑群體的徑向密度分布，徑向和垂向速度色散，以及它們如何隨R變化有關。Bovy和Rix[44]所用的基于作用量的恒星分布模型來自文獻[46]：

(2)銀河系引力勢模型

為計算軌道作用量，Bovy和Rix[44]提出一個完整的三維模型來計算銀河系的勢：一個先驗有許多自由參數的“核球+恒星盤+氣體盤+暗物質暈”模型。

指數盤密度為雙冪律盤模型：

引力勢通過該密度的泊松方程的解析解得到。

核球模型為：

其中，α=1,β=4時，稱為Hernquist模型；α=2,β=4時，稱為Jaffe模型。暗物質暈為冪律模型：

它們的共同引力勢為總引力勢。

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

(3)擬合手續(xù)

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

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

觀測量下的分布為：

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

歸一化后，每個數據點的似然為：

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

數據點是獨立的，考慮異常數據點(outliers)后，總的似然為：

最后根據最大似然法估計，使這個總似然最大，即找到所求參數pΦ,pDF的最合適值。

(4)結果

利用這些模型的所有參數的最佳結果，Bovy和Rix[44]得到了星系盤的徑向標長R d=(2.15±0.14)kpc，與通過光度法推斷的恒星質量的空間分布一致。

來自不同的MAP數據，以不同的半徑約束Σ1.1(R)；因而Bovy和Rix[44]用這43個MAP數據的平均半徑和Σ1.1的平均值的關系，研究了Σ1.1隨半徑的變化，這些MAP的結果也一致[44]。類似地，還得到了1.1kpc的垂向力K z,1.1(R)。使用從文獻[41]距離標度的后驗分布函數所導出的最佳半徑，將這些分別擬合為：

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

3.1.2 Piffl等人的模型

與Bovy和Rix[44]提出的模型類似，同樣用數據限制銀河系的參數，Piffl等人[50]用距銀河系平面約1.5kpc內約200000顆巨星的運動學，來測量太陽附近質量密度的垂向分布，以及局部暗物質暈貢獻。其所用的模型的組分更多，由氣體盤、薄盤、厚盤、核球和暗物質暈五個部分組成，下面介紹其質量模型和分布函數。

(1)質量模型

盤的密度模型為：

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

核球和暗物質暈的密度模型為：

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

(2)分布函數模型

分布函數分為盤和暈兩大部分，

其中，fdisk為盤的分布函數，fhalo為暈的分布函數。Fhalo是暈與盤的質量之比。

盤的分布函數為：

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

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

進一步假設薄盤中的恒星形成率隨時間指數下降，特征時間標度為t0，于是薄盤的分布函數為：

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

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

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

前面幾個準等溫模型及其變形都是基于作用量的分布函數模型。除了基于作用量的，以前也有直接基于E,L z,I3三個運動積分的分布函數模型。關于第三運動積分，大體上也是以Stckel引力勢中的運動積分的形式為藍本，作近似假設。下面簡介其中的一些模型。

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

恒星盤分布函數假定為一個三運動積分模型：

并用Jeans方程進行了檢驗。

3.2.2 一個更多參數的分布函數模型

下面介紹一種參數更多、復雜的分布函數模型[40,42]。Famaey等人[42]構建基于3個運動積分的分布函數f(E,L z,I3)，從而實際星系的恒星盤通過這種形式的分布函數(或者若干此種分布函數的線性組合)來正確地擬合出。他們的方法是，(1)對于I3的計算，簡單套用Stckel勢的I3表達式；(2)對于分布函數的形式，在前人的一個基于2個運動積分的分布函數[40]的基礎上，引入更多自由參數作拓展。

分布函數的具體形式為：

通過計算分布函數下若干階速度的平均，即可得到該分布函數(運動積分空間)的矩：

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

文獻[42]最后介紹了對實際的恒星盤建模的方法：通過將分布函數組分作為基礎函數線性組合。這需要在配置空間中引入一個網格，并最小化方差(如對于分量FΛ代表0階矩，即密度時)：

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

4 討論和總結

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

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

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

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

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

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

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

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

4.3 應用于遠離平衡的星系過程：潮汐星流

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

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

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

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

其中約等號是由于潮汐流形成后，θ0將遠小于正比于時間的那一項。于是有

上式表明海森矩陣D的作用是：把恒星們在作用量空間的分布(參見文獻[36]的圖10)線性映射到角變量空間的分布(參見文獻[36]的圖13)。D是一個對稱矩陣，因此其特征取決于其相互正交的特征向量(n=1,2,3)，和實特征值λn[36]。海森矩陣決定了星流的拉伸方向，即星流沿著最大特征方向拉長[20]：

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

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

4.4 結束語

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

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

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

(1)作用量J

作用量的計算公式為