Stckel引力势理论的拓展及其应用

2021-07-13顾建宇董小波

天文学进展 2021年2期

顾建宇，董小波

(1.中国科学院云南天文台，昆明650216；2.中国科学院大学，北京100049)

1 引言和理论基础

银河系研究的主要目标之一，是测量并计算出银河系的结构，即各种组成成分(如恒星核球、恒星盘、暗物质晕和气体等)的空间分布，以及它们的速度空间信息。为了实现这个目标，在观测上，需要高精度的天体测量学及测光学数据(如Gaia)和光谱数据(如LAMOST)；在理论和计算上，星系动力学建模是核心。星系动力学建模的理论基础是恒星系统的动力学，即研究大量的点质量粒子在它们共同的引力场中的运动[8]。关于银河系建模，尽管星系作为自引力系统，没有星系能处于严格的热力学平衡态或完美的动力学平衡态，然而平衡态动力学模型(或称稳态[steady state]，或准定态[quasi-stationary state])，对于解释银河系和河外星系的观测至关重要。从实际情况看，暗物质晕作为星系引力势的主要提供者，仅在短暂且罕见的主并合事件中，才远离动态平衡，一般情况下该动态平衡仅受吸积或次并合等事件的小扰动。从模型需求的角度看，由动力学测量数据推断出物质分布，该模型要求假设星系处于(准)平衡态；另外，平衡态模型是最简单的模型，更复杂的情况可以被建模为平衡态模型的“扰动”、“变形”或迭加[20]。

运动积分是粒子在轨道运动中的不变量；作用量，则是可以与角变量形成正则坐标的一类特殊的运动积分。由于在“作用量-角变量”坐标下运动方程形式的简单性，作用量成为很有吸引力的变量。并且作用量无论在微扰论方法上还是分布函数的建模上，都是十分必要的。但是，粒子在一般引力势(例如真实的星系)中的运动，基本上不能得到解析的运动积分或作用量；仅当引力势为可分离势时，能解析地得到。Stckel势是椭球坐标系下可分离的引力势形式，该引力势下的运动积分可被解析地计算[1]。在应用上，通过将一般引力势近似或拟合为Stckel势，人们可估算一般引力势下的运动积分。近十几年来，一些研究者基于Stckel引力势理论中的解析式，用(参见3.2节)解析或数值的(例如Binney的工作)方法来作拓展，即基于Stckel势来估算一般星系中的作用量(或运动积分)。特别是最近几年来，Binney和他的合作者们把他们的方法都实现为计算机程序，并以公开源码的方式释放了(详见第2章)。

接下来，第1章的其他小节分别介绍：恒星系统的动力学及统计力学的概要、运动积分以及Stckel势的相关基本背景知识[1]。第2章介绍估算一般引力势的运动积分和作用量的方法：2.1节介绍Stckel Fudge作弊法[13,20,22,24]；2.2节介绍Stckel引力势拟合法[19]，以及更早期的工作；2.3节介绍一类环面坐标变换的作用量估算方法，基于Stckel玩具势模型得到一般引力势下的作用量[20,21,24]；2.4节是我们基于文献[24]所介绍的TACT程序对各种方法做的比较；2.5节则是介绍一些其他求运动积分的方法和程序。第3章介绍近几年来这些方法在银河系数据上的应用(如文献[44])。第4章进行讨论和总结，特别是简介“作用量-角变量”方法在不可积的星系现象(见4.2节)和非平衡的星系过程(见4.3节)中应用的最新进展。

1.1 恒星系统的理论基础

本文讨论的单粒子在整体势的运动，是从引力多体运动简化而来的，即单粒子近似或称平均场近似[60]。三维背景空间中一个包含有N个粒子的经典系统，可用哈密顿方程描述：

对于该系统所有粒子的坐标和动量所构成的全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粒子联合分布函数。对上述刘维尔方程作单粒子近似f2(q1,p1,q2,p2,t)=f1(q1,p1,t)f1(q2,p2,t)①这是一种平均场近似，即将粒子间相互作用当成总背景场，而不考虑粒子间相互作用和关联；体现在分布函数上就是每个粒子的单粒子分布是无关联的，因此乘积为它们的联合分布。，可得到玻尔兹曼方程②由玻尔兹曼于1872年提出。从刘维尔到玻尔兹曼方程的约化，中间还有一个BBGKY级联方程。：

其中，这里f=f(q,p,t)是单粒子分布函数，即处于自由度3的坐标和速度所构成的6维相空间(一般称为µ相空间)；v代表速度，F代表受力③该受力是对于该“单粒子”而言的，施力者为“其他粒子”(实际上是分布函数)的作用(因而在这里直接体现平均场近似的思想)，以及可能存在的外场。玻尔兹曼方程可以用来描述非自洽的体系(自洽体系是指粒子所受力全部来自于研究系统本身，如自引力系统)，即通过建立多种组分的玻尔兹曼方程来描述。；σ(Ω)是分子微分散射截面，这里考虑了二体碰撞；对于其他作用类型(如刚性)的散射，将可能不再是这样的形式。玻尔兹曼方程体系中所需要的假设非常繁杂(甚至冗余)，而且至今仍有难以严格数学证明的规律。本文对玻尔兹曼方程仅作介绍而不作分析(更多参见文献[8]中的第4,7章和文献[60]中的第3,4章)。此外可以验证，刘维尔方程在Γ相空间的分布函数的熵是不变的，而玻尔兹曼方程在µ相空间的分布函数的熵是不减的。

对玻尔兹曼方程进行局域平衡假设，碰撞项将消失。其实，在星系这种具体的系统中，恒星的半径与星系尺度相比很小，其平均自由程远大于系统尺度，碰撞时间远大于穿越整个星系的时间，比宇宙年龄还大很多量级，因此碰撞项本应忽略不计。把玻尔兹曼方程两边乘上一个碰撞不变量ξ(x,v,t)(如质量、动量、能量等)，然后作局域平均(由于碰撞前后的某些对称性,局域平均后碰撞项不做贡献，见文献[60]第7A章)，可分别得到流体力学性质的方程(注：流体力学过程是局域平衡过程)。例如对于动量守恒，玻尔兹曼方程两边乘以ξ(x,v,t)=mv i(i=x,y,z)，通过对动量空间积分并利用动量空间的散度定理(无穷远处分布函数为0，limp→∞f=0)，可得到

此式实质上等同于流体力学动量守恒方程：

这样，我们就有了“N体系统(力学描述)-刘维尔方程(6N维)-玻尔兹曼方程(6维)-流体力学方程(流场)”的联系。

另一方面，从无碰撞玻尔兹曼方程出发，两边乘以速度分量，再对速度空间积分并利用散度定理(与得到式(4)的方法类似)，我们得到Jeans方程：

其中，ν代表位形空间数密度，σij为速度弥散张量。Jeans方程不仅通过平均简化了方程，还将分布函数与这些可观测量联系起来，从而对联系观测天文的工作很有帮助。Jeans方程还给人们提供了一种检测分布函数正确性的方法。当然,Jeans方程仍然不是闭合的。

在µ相空间中的粒子分布函数f=f(x,v,t)①由于p=mv(m为单粒子的质量)，而且在星系动力学的粗糙建模中常常将测试天体考虑为无穷多同质粒子的分布，且为方便一般假设测试粒子的质量为1，于是下面我们这里的习惯，以下讨论的f都是关于x,v的分布函数。，跟粒子视角的单粒子近似下的运动(x(t),v(t))有着密切的联系。这种单粒子分布函数，实际上假设体系是由无穷多的点粒子构成，某种程度上与场论的意义比较像(场论就是处理无穷多自由度)。从场的角度来理解，这是把所有粒子之间的相互作用当做背景场，然后对于任一粒子，只须研究它在这个场中的运动，不必考虑其他粒子对它的作用；此即前述平均场的思想。但是，实际上想详细了解星系结构与演化(如分布函数)，直接从自引力系统的统计物理方程出发去求解分布函数并不容易，而且运用最大熵原理之类的方法也是不够的(由于星系作为长程相互作用系统的特性，见第4章)。为得到分布函数，一种方法是N体模拟，基于模拟的结果作统计。另一种方法是根据粒子在引力势的运动特性来构建分布函数模型[8]，以及其他理论手段(详见文献[8])；在这种方法中，对运动积分的研究是个基础的环节。下面，我们主要从哈密顿力学以及单粒子近似下的恒星轨道运动的角度出发，介绍恒星在静态引力势下的运动积分的相关内容。

1.2 运动积分和作用量-角变量体系简介

1.2.1 哈密顿-雅克比方程

力学系统在相空间的正则坐标(q,p)的演化，可以看做是无穷多个连续以哈密顿量为生成函数的无穷小正则变换的结果。在正则变换中令新的哈密顿量为0，并选取第二类生成函数F2②第二类生成函数是指：在正则变换中，取独立的广义坐标变量为老的坐标和新的动量，亦即生成函数为老坐标、新动量以及时间的函数。可使某种新的正则坐标Q,P均为常数，即Q=α,P=β，可得正则变换的方程：

被称为哈密顿-雅克比方程(哈雅方程)，其解的形式为S=S(q,β)就是生成函数F2。另外，第二类正则变换中哈密顿函数之间的微分关系表示为：

1.2.2 运动积分

如果一个力学量不显含时间而且保持不变，h(q,p)=const，就称之为系统的运动积分，即运动积分是粒子在轨道运动中的不变量。寻找系统的这种不变量，具有重要意义。对于一个自治哈密顿系统，哈密顿量不显含时间。此时哈密顿量本身就是一个运动积分，其积分常数E等于系统总能量。利用公式(7)和(8)，哈雅方程的解可写成时间项和广义坐标项分离的形式，表示为：

式中的W称为哈密顿特征函数。对于自由度n的自洽哈密顿系统，最多有n个两两独立且对合①如果两个运动积分的泊松括号为0，那么它们是对合的。的运动积分，因此在上式的积分常数β={β1,...,βn}和h这(n+1)个常数中，至少有一个不独立，即通过h=h(β1,...,βn)相联系，我们取βn=h。对式(9)两边微分得到两个微分关系，以及生成函数的微分关系，即式(8)；再将W代入S中，可得到关系W与新旧相空间坐标的微分关系：

现在我们以哈密顿特征函数W作为生成函数作正则变换(q,p)→(Q,P)，其中，

此时的新哈密顿量仍是K=H，并且新哈密顿函数只含新广义动量而不含新广义坐标，该正则变换是作为后面讨论运动积分坐标的基础。

更进一步，对于自治(完全)可分离的生成函数W，可写为：

则相应地可以分离出n个哈雅方程(不对i求和)，可表示为：

哈雅方程方法正可方便可分离系统的求解。

如上所述，若自由度n的自治哈密顿系统为(完全)可积的，则它存在n个互相对合的运动积分I i(i=1,2,3,···,n)(我们用I i代表已经从E,β中选好的运动积分)；此即(刘维尔)可积性定理的内容。另一方面，根据刘维尔保体积定理，以这些运动积分为某一组常数时的坐标集合，M={(q,p)|I i(q,p)=βi}，它在以I i为哈密顿函数的相流下保持不变，并且同胚于n维环面T n(若不存在共振)，因此将(完全)可积系统的轨道转化为一个n维准周期运动[63]。

关于运动积分，值得注意两点：第一，在星系动力学语境中，运动积分和运动常量是作了严格区分的两个概念；第二，关于运动积分还有更多细分的概念。对于恒星在星系引力势中的运动，给定引力场，运动常量的定义是：关于相空间坐标和时间的、沿恒星轨道不变的任意函数，即C(x(t),v(t),t)=const。运动积分也是在恒星轨道中保持常数的量(即运动常量)，但它进一步要求是不显含时间的函数，I(x(t),v(t))=const。运动积分可以分为孤立积分和非孤立积分两类[8]。孤立积分能影响恒星轨道在相空间的分布，即相空间轨迹(phase trajectory)所占据空间的维度；对于可积系统，每增加一个孤立积分，就把相空间轨迹所占据流形(即相空间的子空间)的维度降低一维。非孤立积分则对于恒星轨道在相空间的分布没有影响，因而没有实用价值。

另外，在完全可积的或部分可积的系统中，请注意对合积分与孤立积分这两个概念的区别。对于自由度为n的自治哈密顿系统，相空间为2n维，最多有2n−1个独立的(即线性无关的)运动积分，其中孤立运动积分最多可以是2n−1个；而对合是与可积性有关的概念，两两对合的运动积分不超过n个(当对合积分的个数为n时，该系统就是刘维尔当年所定义的完全可积的)。因此，若有一组运动积分是互相对合的，那么它们肯定是孤立积分；但反之不然。例如，对于球对称引力势中的恒星轨道(6维相空间)，运动积分有5个，其中孤立积分可以是4个或者5个(仅当轨道为开普勒型时)，但相互对合的积分的个数是3[8]。对于星系动力学中常见的轴对称引力势系统(例如盘星系)，一般有3个孤立运动积分，能量E、角动量的垂向(z)分量L z、以及第三个运动积分I3，但I3往往是没有解析形式的；对于三轴引力势系统，一般有3个孤立运动积分，能量E和I2,I3，但I2,I3往往是没有解析形式的。过去，前人有很多工作就是在找第2、第3运动积分，以及在研究稳态系统的恒星分布时所用的、以此为自变量的分布函数f(I)(这也是本综述文章的中心内容)。

1.2.3 作用量-角变量体系

尽管运动积分I与其广义坐标(记为θ，角变量)已经将运动化为准周期运动，但是运动积分一般不是正则坐标，即它与其广义坐标不满足哈密顿共轭，运动积分经过变换后能够满足共轭的叫做作用量J。下面我们对角变量和作用量的要求作一下简介(可见文献[8]第3.5.1节)。

如前所述，我们假定存在从I到J的可逆的变换(det(∂J/∂I)̸=0)，使作用量J满足：(1)运动积分，(2)正则坐标。于是哈密顿方程(1)可用，关于J的部分有：

于是哈密顿量只是作用量的函数，表示为：

其中，Ωi(J)是角频率。对于束缚轨道，角变量θ和运动积分所描述的就是一个n维准周期运动，即前述的同胚于n维环面T n的流形M；M上某一维的任一闭合环路γi(同胚于圆周S1)，由某一组(θi,J i)所代表，并且要求θ沿γi的环路积分为2π。

于是我们可以将相空间轨道展开成傅里叶级数：

角频率之间的比例关系丰富了体系周期类型和遍历性的复杂性(见第4章)。天体动力系统的周期特性一般来说非常强，这也是作用量-角变量体系在天体力学中经久不衰的重要原因。

其中第2个等式是由于正则变换的相体积不变性；第4个等式是转化为边界积分，并且刨除积分奇点的邻域。将奇点置于作用量-角变量坐标的中心，=0，此时等于作用量的不变量。于是，我们由正则变换的庞加莱不变量的不变性，可得：

其中，T是正则变换的矩阵，指在γi对应的相空间所在区域。也由此可见，某个作用量的值J i代表子环γi所围的面积(当对q i,p i选取合适的标度)，而θi代表转过的角度。

此外，我们可用式(15)计算角频率，用式(11)中第1行的哈密顿特征函数计算角变量的值。数值计算上见附录A。

于是，从一般的相空间坐标到作用量的关系为：相空间坐标(q,p)→运动积分I和角坐标θ→正则环面坐标(J,θ)。

1.2.4 基于守恒量的分布函数

关于单粒子分布函数f(x,v)(此时为3个自由度的，由坐标、速度所构成的6维µ相空间)，对于(准)稳态体系(指∂f/∂t=0)，由无碰撞玻尔兹曼方程可推导，分布函数可写成运动积分的函数f=f(I)，此即Jeans定理；而对于规则轨道具有非谐振频率的准稳态恒星系统的分布函数，可写成3个作用量的函数f=f(J)，此即强Jeans定理。因此人们能更加不冗余、深刻地理解恒星的分布。

早先前人所用的分布函数模型，常常是f(E,···)的函数形式，即把能量作为分布函数的参数(参见3.2节)。不幸的是，对于实际的建模工作而言，这是一个陷阱，从而阻滞了基于运动积分或作用量建模这项事业的进展[15]。对于实际星系的动力学建模，不能将E放入分布函数的参数列表中。原因很简单：实际星系是一个多组分系统，动力学模型中至少要包括暗物质晕和恒星组分。轨道的能量不是局域量，如果E出现在参数列表中，那么该轨道周围的相空间密度将以不受控制的方式变化。例如，最简单的轨道是恒星位于银河系中心的规则轨道，当添加暗物质晕组分后，中心引力势会变深，位于中心的恒星的能量会减少，而该轨道保持不变。但是，只有当模型完全组装好并且确定了其势能之后，我们才能知道该轨道的能量E0=Φ(0)。因此我们无法在建模工作一开始，就指定恒星的f(E,···)。而当分布函数只是作用量的函数时，每个轨道周围的相空间密度则不会改变[15,51]。

总之，相比于一般运动积分，作用量-角变量体系的优势是：作用量是一种特殊的运动积分，而且作为正则坐标，作用量有共轭坐标即角变量；作用量体系适用于使用微扰方法，也正适合描述周期性强的天体运动；作用量是绝热不变量，可用于缓慢变化的引力势，以及邻近势能中的轨道识别；此外如上段所述，使用以作用量为参数的分布函数来为星系(多组分系统)建模，才使得基于分布函数的建模成为可能。

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τ(τ)前面的正负符号的规定，也可能有所出入。本文采用如式(21)所示的方式。。令Fτ(λ)=4(λ+α)(λ+β)(λ+γ)ζ(λ)，Fτ(µ)及Fτ(ν)依此类推②为了书写方便，本文后面统一用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画图程序。关于算法的更多细节，请参见文献[20,24]，以及TACT程序的源码。

2.1 Stckel作弊法

关于估算一般引力势系统的作用量-角变量，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节)，其中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 v(v)可以重写为F v(π/2)+δF v(v)。

这里u0是一个不影响作用量计算的参考值，其选取是固定v而使δF u(u)最小值时的u(见文献[13])。F u(u0)和F v(π/2)的计算为：

本小节主要介绍Sanders和Binney[22]提出的算法，其是对Binney[13]提出的Stckel作弊法的三轴推广。对于一般轴对称势系统，除了Stckel作弊法，Sanders[19]还提供了对一般轴对称势的Stckel拟合的方法(下一小节介绍)。但是三轴Stckel势比轴对称的多一维参数，对其拟合所需的具体拟合模型会复杂许多，于是Sanders和Binney推广出了三轴势的Stckel作弊法。该方法借助类似于运动积分J,K的另一套运动积分(作弊法得到的新的不独立的运动积分)，以及局部近似的Stckel势，用动量和运动积分的关系式得到每一点的动量pτ，再积分得到作用量J；类似地，可求出角变量θ。本方法以当前真实的相空间点做初始条件(x0,v0)，来计算一系列动量，这样就能不需要Stckel近似势的分离势F(τ)的显式表达式；然后积分得到作用量；最小化作用量的标准差(与每个初始(x0,v0)有关)，调整焦距，从而再得到合适的作用量；对于轴对称势系统，还可做迭代以减小误差。下面介绍其具体步骤。

对于一个一般引力势，定义三轴情形的辅助函数：

当固定另外两个椭球坐标时，某个τ对应的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给出了数值调整合适焦距的方法。例如对于短轴闭环轨道(short axis closed loop orbits)，轨道形成在z=0面，焦距为y=±∆1的椭圆。在这种轨道，只有作用量Jν不为0，利用这一点，给每个能量值找出一个合适的焦距。

其步骤是，先给定一个一般势，最大√和最小能量值以及期间的能量网格点；在每个网格点E i，画一个从坐标(0,y k,0)以速度为垂直出发的半周期的轨道；记录所画轨道的相空间点并计算作用量Jν，求出使作用量的方差最小的焦距∆1。∆2的计算也可用上述方法求出。

(3)环面迭代增加精度

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)，非常影响作用量的标准差，因此通过一定的手续来估算一个更好的焦距∆。Sanders[19]采用的估计方法是，由于Stckel势的分离性，可得

由上式以及两种坐标之间的变换关系，可推得[19]：

既然假设所估算的局部势Φ就是该Stckel的势ΦS，于是就用该式来估计a2的值，其中的偏导用已知的真实Φ来数值计算梯度而得；每个a2值，需要在轨道中每个点的相应值来计算，最后做平均，这些点在下一步的轨道采样中获得。

(3)轨道采样以及运动区域的估算。进行轨道积分，即在(x,v)坐标下通过每一点的势和时间积分来算下一步的坐标(x′,v′)，即相当于数值法画轨道。Sanders和Binney的TACT程序中用守恒格式的数值方法画出轨道，并且每个轨道周期离散为300个点。

在画轨道时：1)在前几周期的一些点中，隔5个点用式(48)估算a2(c2已取1)并平均，最后得到那个最合适的扁球坐标；2)通过前后∂τ/∂t异号等条件(其中每次需要x到τ的坐标变换)来判断是否达到轨道的τ的端点，找到后就记录，得到λ+−,ν+,ν−=c2)，从而得到轨道在椭球坐标的所在区域。

假如所拟合的势Φ是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 基于环面映射的方法

本节介绍的内容包括：(1)McGill和Binney[11]提出、Binney及其合作者后来完善的，通过正则变换构造环面(作用量-角变量)来计算作用量的方法，即基于玩具模型的环面映射(torus mapping,TM)思想及其TM程序[12]；(2)基于此思想的两个扩展：迭代环面构造法(iterative torus construction,ItTC)[22,24]和轨道积分拟合生成函数法(orbit integration to generation function,O2GF)[21,24]。这类方法的特点是精度最高，计算花费大，属于收敛方法，即可以通过充分大的计算量，把作用量的误差一直不断地减小下去。

将环面(n-Torus)视为几何物体，而将解析的哈密顿量的环面变形为所求的目标哈密顿量的环面，如同微扰理论，通过寻求从玩具环面坐标到目标环面坐标的正则变换来实现[11]。即用已知的解析引力模型作为玩具模型，来逼近一般的目标引力系统，而这是在环面坐标的意义上进行的，并且通过正则变换来实现，此即环面映射的关键思想。一般所用的玩具模型为等时势或谐振子势，如式(70)和(71)，这些都是Stckel势，可解析求解。尽管环面映射法不同于前述的通过近似Stckel的分离函数Fτ(τ)的Stckel作弊法，但其需要从某种作为玩具模型的Stckel势出发，本质上相当于某种简化版的微扰法，因而仍可看作是Stckel势的应用拓展。

本节中，若无特殊说明，用(Jtoy,θtoy)代表玩具环面(即玩具势下粒子运动的环面坐标)，(J,θ)代表目标环面(即所求的一般引力势下粒子运动的环面坐标)。从玩具环面映射到目标环面，将此看作一个正则变换，选取第二类生成函数S=S(J,θtoy)，将这个生成函数对角变量傅里叶展开，其一般形式为：

其中，这里的虚单位i是为后面求导表达式的简便而加的。当沿着玩具环面在θ的任一分量增加一个周期2π，相应地，目标圆环上的像点也应完成一个周期回路，最终这要求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)求微分，得

这样就通过生成函数联系起了玩具环面和目标环面。

2.3.1 环面映射程序

Binney和McMillan[12]提供了TM程序(Torus Mapper)①TM程序的源码可下载于：https://github.com/PaulMc Millan-Astro/Torus。，该程序实现三个方面的功能，具体功能与输入输出如下：用户写出目标引力势及其一阶微分，以及给定作用量(J r,Jϕ,J z)，1)对于特定θ，可计算(x,v)；2)计算雅克比∂(x)/∂(θ)；3)给定(x,J)，计算θ。与传统的恒星轨道相比，环面的优点是：一个环面是由100个数字指定的，即S n,∂S n/∂J和一些其他参数，而不是由沿时间序列的数千个相空间坐标指定，因此可以在给定分辨率下显著缩小轨道库的大小[12]。

环面映射方法的流程见式(68)[12]，通过三步正则变换完成，变换形式为：

已知的真实势和轨道，给定目标作用量J和系数S n的试探值，由玩具坐标θtoy的网格点，根据式(66)可得到相应的玩具作用量Jtoy。在玩具势下，可直接通过解析的方法从环面坐标得到相空间坐标(xtoy,vtoy)。有了这些相空间坐标可计算相应采样点的能量，之后用Levenberg-Marquardt算法迭代以最小化哈密顿量偏差δH，得到最终的S n。类似地，可得到∂S n/∂J，再通过式(67)可得到目标角变量θ[24]。

这样，通过真实的环面坐标(目标环面坐标作为已知)，用合适的系数S n(J)和∂S n/∂J，转化到玩具环面坐标；再通过玩具势下的解析性就得到相空间坐标(xtoy,vtoy)，并由点变换得到(x,v)。点变换也是一类正则变换，在程序中仅当J z/J r>0.05且尝试在不进行点变换的情况下拟合环面失败时，再使用点变换[12]。

文献[12]中的附录图A1给出了环面映射方法的原理示意，对于一个好的玩具模型，更高阶的生成函数傅里叶系数能越来越好地拟合轨道。

2.3.2 迭代环面构造法

如上所述，TM程序[11,12]给出了一般轴对称系统的相空间坐标和环面坐标转换的功能，主要是由(J,θ)计算(x,v)。对于环面坐标的计算，除了按照上述方法，重复地构造环面，直到得到的相空间坐标穿过给定相空间坐标，TACT程序的其中一个方法对其进行了改进，即迭代环面构造法(ItTC)[22,24]。该方法是“Stckel作弊法+TM程序”的迭代算法，即用Stckel作弊法计算试探的环面坐标，然后寻找环面上与其最近的相空间点，最后补上残差并迭代[22]。

给定一个初始相空间点(x,v)，先用Stckel作弊法计算一个环面坐标(J S,θS)；然后用环面映射得到一个以此为作用量的环面，该环面的频率为Ω，在该环面上找到一个与给定点最近的相空间点(x S,v S)，即

使其最小。之后，计算所找到的这一点的环面坐标J P(x S,v S),θP(x S,v S)，得到∆J=J P−J S作为误差，假设这个误差变化很慢，于是将其作为残差补回去，得到一个更好的作用量J2=J S+∆J。

可以再次构造一个以J2为作用量的环面，然后继续从中找到一个与给定初始点最近的相空间点，……，如此往复迭代，得到更精确的作用量，这个过程是收敛的。

2.3.3 轨道积分拟合生成函数法

这里介绍Sanders和Binney[21,24]给出的轨道积分拟合生成函数法。该方法过程为：轨道时序采样(x,v)(t i)→计算玩具环面坐标(Jtoy,θtoy)(t i)→解生成函数系数矩阵。它和TM程序的区别是，该方法将环面映射思想应用于一段运动过程(即一段时间序列)，计算出在时间t上沿轨道采样的一系列相空间坐标的作用量、角变量和角频率，即构造生成函数是基于轨道积分，而不是TM所基于的某个时间点的哈密顿量偏差δH的最小二乘。

具体步骤为，给定一个真实相空间点作为初始点，首先计算一个时间N TT的轨道积分(画轨道)，并采集Nsamp个相空间采样点，其中T是具有相同能量的大圆轨道的周期；然后在每个相空间点，解析计算其玩具势中的(θtoy,Jtoy)；之后将这些玩具环面坐标点代入到式(66)和(67)来求解整个方程组的未知量(目标作用量J，生成函数傅里叶系数S n及其导数∂S n/∂J，以及目标角频率Ω和常数θ(0))。

(1)轨道积分与轨道分类

进行轨道积分。在没有转动的情况下，一般三轴势允许两类基本的非共振轨道：环轨道和盒轨道[1,21]。同一真实势，不同的初始条件可能导致不同类型的轨道，于是TACT为不同类型的轨道匹配相应的玩具势，与给定轨道类的环面相同的几何结构。

(2)玩具势的拟合

对于环轨道，玩具势为等时势：

对于盒轨道，玩具势为谐振子势：

其判定依据为轨道类型，即看L x和L z的方向变化。在用玩具引力势之前，需要在轨道上采集样点，拟合出相应玩具引力势的质量和标长参数。

(3)生成函数

在每个相空间样点t i计算相应的玩具作用量-角变量(θtoy(t i),Jtoy(t i))，并产生相应的式(66)，其中各J和S n未知。Sanders和Binney轨道积分拟合生成函数法的实现技巧为：“我们不能精确地解这些方程，因为处理的方程有无穷多个未知数。我们可以在每个方程的右边只包含有限个项，则方程右边和左边不能精确一致，因此可以采取的正确方法是：对每个方程的左右边的差作平方，然后对所有方程的差平方求和(即轨道积分)，并最小化这个总的平方和”[21]。其总平方和的表达式是：

其中，k=1,2,3代表坐标维数，N限制为有限维向量(有限精度)，|n|≤Nmax，Nmax≈6；(t i)代表样点，最外层的求和(关于指标i)是遍历所有采样点，此即“轨道积分”的含义。更多关于样点个数的选取原则请详见参考文献[21]中的2.2和2.3节。

下面最小化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作弊法，其中，Fudge v1用式(48)估算焦距，Fudge v2基于壳轨道(shell orbits)来估算焦距；Fit对应Stckel势拟合法，ItTC对应环面迭代法，O2GF对应轨道积分拟合生成函数法，AvGF为平均生成函数法；CAA和SAA分别为柱绝热近似和球绝热近似方法，下一节作简介。轨道中的每个点都独立计算一次作用量，因而作用量的数值随时间有所波动。

图2 作用量数值计算的效果对比

Sanders和Binney[24]对几种方法的在轨道类型、精度和时间等方面的性能做了更详细的对比(见其中的图2到图6)。总体上看，收敛方法比非收敛方法有更高的精度，但耗时也高1∼3个量级。其中，轨道积分拟合生成函数法精度最高，耗时较长；Stckel作弊法速度很快，是大量计算的最经济方法，但精度不如收敛性方法高，特别是对于盒轨道的误差大。另外，对于轴对称势经过环面迭代后(即ItTC法，目前TM程序只提供了轴对称的环面迭代)，Stckel作弊法的精度会提高，而耗时也变长。此外，TACT程序包中，轨道积分拟合生成函数法和Stckel作弊法还能用于三轴势。

2.5 其他一些求运动积分的方法或程序

前面两节我们主要介绍了Binney及其合作者们近年来的工作。接下来我们介绍其他一些基于Stckel势的估算运动积分或作用量的方法。

2.5.1 “近似的运动积分”解析式

若Φ在(−α,−β,−γ)连续，应有Φ(−α,−β,−γ)=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程序包以及其环面迭代功能所依赖的TM程序包，星系动力学还有一些程序可供使用，这些程序现在都可以公开获取。主要程序如下所述。

(1)Gadget，Gadget①https://wwwmpa.mpa-garching.mpg.de/galform/gadget/index.shtml以及其后续发展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 深度学习在作用量计算中的应用

在这一小节中，我们介绍最新的一个研究动态。Ibata等人[32]将人工智能(AI)中的深度学习方法应用到作用量计算上。该文用AI算法来实现在静态引力势下从粒子在相空间的“位置-速度”到“作用量-角变量”的映射，并且得到加速度场。该算法的核心原理与前述“生成函数法”类似，先从数据点中拟合玩具势；然后通过从玩具环面(θ,J)到目标环面(θ′,J′)的正则变换表达式，迭代出最合适的生成函数G；之后是正则点变换的生成函数P，最终得到所求的目标环面(θ′′,J′′)。其中两个正则变换所用网络的形式为：

该算法的程序为ACTIONFINDER，基于Python的Pytorch所写。据文献[32]介绍，经测试ACTIONFINDER得到的作用量精度要比Stckel作弊法要高，而且它还有如下优点：相比于早期的Torus Mapper对初值要求更不敏感；不需要向它提供系统内轨道的公平样本，几乎可覆盖任何感兴趣区域的样本；不需要预先知道所研究系统的哈密顿量或引力势；此外，该算法所用的深层神经网络，相比前述“生成函数法”线性分解成的傅里叶系数，能够生成更加灵活的非线性函数，并且这种灵活性可能使其更容易拟合更一般的动力系统。看来，将来经过一定的发展完善后，深度学习或其他AI技术用来计算作用量，适用性会比传统方法更广泛、性能更佳，应该有一定的发展前景，未来或成为分析更复杂体系的另一种尝试。

3 应用

在星系研究中，作用量的一个重要应用在分布函数方面，即使用以作用量作为自变量的分布函数，为星系建模[15]。一般来说，我们难以直接获得分布函数的形式。在建模以及数据拟合的实际工作中，人们通常是为不同的星系组分(如恒星盘、核球、暗物质晕等)采用某种先验的分布函数形式，然后根据数据拟合来确定其中的具体参数的取值。

大约10年前，基于分布函数的建模方法实际上很难应用(见Read[61]和Binney[15]的综述)，星系的建模工作往往是基于Jeans方程，即对分布函数所满足的无碰撞玻尔兹曼方程取矩(即所谓“矩方法”)。假定星系处于稳态(即∂f/∂t=0)，从无碰撞玻尔兹曼方程出发，分别乘上速度分量并且积分掉速度空间(即取矩)，即得到Jeans方程。对于实际的银河系建模工作，研究者还会假设轴对称(即在柱极坐标下∂f/∂ϕ=0,∂f/∂vϕ=0)，那么Jeans方程的表达式为(柱极坐标)：

其中，ν为示踪者恒星的数密度，为分布函数的0阶矩；平均速度为1阶矩；速度弥散张量为 ∫

为2阶矩。在观测上，速度弥散可由某一个位置处诸恒星的速度方差来计算。

正如Read所总结[61]，Jeans方程和矩方法的优点是：与其他方法相比非常快，可以探索较大的参数空间；并且不需要假设分布函数的形式，因为只是限制了它的矩。缺点是：必须对数据进行整理以计算速度的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Σ。

在下面的几个小节中，本文介绍近年来有代表性的建模工作。这几个工作都是基于3个运动积分(或作用量)的分布函数来为银河系建模。这些工作中，对于观测到的恒星盘成分，一般是使用基于作用量的准等温球分布函数模型，或者基于运动积分的准等温球分布函数。对于暗物质晕，这些工作使用质量分布模型(如NFW模型)而不是分布函数模型。后来，Piffl等人[51]及Binney和Piffl[49]开始尝试使用分布函数模型来直接为暗物质晕建模并拟合银河系数据；这项暗物质晕建模工作由于数据有限以及软件上的问题，目前的效果不太理想[15]。在本节的末尾，我们还介绍前人提出的针对恒星盘的另一种分布函数模型。

前人的这些工作，主要是基于RAVE和SDSS的银河系巡天数据。现在，在Gaia的天测和测光数据以及大规模光谱巡天(例如中国的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)拟合手续

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是为了在圆盘中创建中心腔，而其他两个盘则设置为零。

核球和暗物质晕的密度模型为：

因此Piffl等人[50]的质量模型具有8个参数，薄盘和厚盘分别具有3个自由参数(Σ0,R d,z d)，而暗物质晕具有2个自由参数ρ0,r0，其中薄盘和厚盘的径向标长取相同。文献[50]中的表1给出了银河系模型的固定参数取值。

(2)分布函数模型

分布函数分为盘和晕两大部分，

其中，fdisk为盘的分布函数，fhalo为晕的分布函数。Fhalo是晕与盘的质量之比。

盘的分布函数为：

其中，fthin为薄盘的分布函数，fthick为厚盘的分布函数，Fthick是厚盘与薄盘的质量之比。

其中，ν仍代表垂向频率而不是椭球坐标。但是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章给出了一种方法，通过由积分区域边界的线性组合表示积分区域中的点，来计算其各阶矩中的积分；本文附录B简介了其基本思想。有了矩的计算，文献[42]的第5章介绍并分析了分布函数中各种参数所体现的物理特征。

文献[42]最后介绍了对实际的恒星盘建模的方法：通过将分布函数组分作为基础函数线性组合。这需要在配置空间中引入一个网格，并最小化方差(如对于分量FΛ代表0阶矩，即密度时)：

其中，Λ=(α1,α2,β,δ,η,z0,s,a)为分布函数的参数集，cΛ为系数。拟合程序中不断加入满足要求特征的新组分，直至χ2的最小值不显著变化。

4 讨论和总结

前面内容聚焦于准稳态时星系的分布函数特点及建模。但是，与常见的稳态、准稳态的热力学系统相比，自引力系统有些特点本质不同，这些特点也许会限制前述方法的适用范围。另一方面，当考虑处于非稳态(远离平衡时)的星系结构和演化过程，或者考虑少体引力系统(粒子数较少的引力系统)中的具体恒星的运动，上述基于作用量的建模方法是否有效，或者说如何使得上述方法仍能应用于某些具体场合，这是非常值得探索的课题。

在本章里，我们将首先在动力系统理论的视角下(dynamical-systems theory)，研究自引力系统的动力学(见4.1节)；然后介绍和评述近年来研究者把“作用量-角变量”方法应用到完全不可积的星系现象和星系结构的进展(共振和星系棒，见4.2节)，应用于非平衡的星系过程——潮汐星流(tidal streams)的研究，包括星流的形成过程，在作用量空间、角变量空间对星流的刻画，以及把星流用来拟合主星系的引力势分布(见4.3节)；最后是全文的总结(见4.4节)。

4.1 动力系统理论视角下的研究

从统计物理角度，引力N体系统是一种长程相互作用系统。与一般的近程作用系统不同，这类系统中任何粒子之间的作用几乎都不能忽略(粒子间的作用势ψ∝r−α，且α小于位形空间维数d；假设粒子在全空间均匀分布，那么总作用势将是发散的)，并且与势能有关的宏观量不再具有统计物理所经常要求的可加性。另一方面，自引力系统由于可形成核-晕结构使其熵可能不具有最大值，因此并没有一个真正的平衡态；一般来说，相比于近程作用系统，自引力系统的弛豫过程可以更快，例如剧烈弛豫机制(violent relaxation)[55,58]。这些特点导致自引力系统与统计物理中常见的近程作用系统显著不同，往往不能用平衡态统计方法来很好地研究。

从动力学角度，作用量-角变量坐标将n自由度可积系统的2n维∑相空间变成n维环面T n，如果存在一组非零整数k1,k2,k3,···,k n，使得角频率Ω满足时，则存在共振。共振关系可有不同多种。庞加莱(Poincar)截面法就是在相空间(即坐标维+速度维)中取一个低于相空间维度的“面”，记录系统的轨道穿过这个面的点，并研究这些落点规律，该方法是降维研究运动的常见且直观的方法。以3自由度可积系统为例(如在Stckel势形式的星系中的任意一个恒星轨道)，该系统在相空间中呈现为一个三维环体(3-torus)；若系统的其中两个频率比Ω1/Ω2为有理数，则在其相应的二维环面的运动是周期的，其在环面上的运动轨迹为周期轨线，在庞加莱截面上则是一些离散的不动点；若运动是准周期的，则是布满环面的准周期轨线，落点在截面上将形成一条闭合曲线。

当可积系统受不可积微扰且成为近可积系统(近可积系统是大部分的相空间仍被规则轨道占据的受扰可积系统)，KAM定理给出了近可积系统中仍然稳定存在局部的不变环面的条件。KAM定理是20世纪中叶Kolmogorov,Arnold和Moser提出的重要定理，后续的工作放宽了定理的条件。对于截面上的不动点，当施加扰动，简单来说，经过多次运动后，双曲(不稳定)不动点将逐渐偏离，椭圆(稳定)不动点则将围绕不动点形成闭合环面，或在满足KAM定理的一定区域内(即KAM环面)稳定存在，或形成更高阶的新的两种不动点。如此迭代，将导致更复杂和不可预知的演化，从而产生混沌。此外，对于n>2自由度的近可积系统，相空间为2n维，等能面(或称能量流形)为(2n−1)维，而n维的KAM环面不再能将等能面分割成闭区域的集合，即存在阿诺德扩散现象(Arnold diffusion)。阿诺德扩散使得n≥3自由度的动力系统的行为更为复杂，在稳定的轨道之间存在着高度不稳定的轨道，这样运动轨迹就可以绕过KAM环面到其他区域，仍会出现混沌；人们对于阿诺德扩散的理解，目前仍然非常有限。大致上，我们可以说，有理环面的破坏源于运动模式的非线性共振，非线性共振进一步导致混沌运动[63]。现代的“动力系统”学科本身也在拓扑和分析上有了更严密的定义和更多的理论，我们期待着这些现代理论能在自引力系统的动力学研究中大放异彩。

4.2 应用于“不可积”的星系现象

绝大多数自引力系统(如星系)都不是可积系统；即使数值计算，由于混沌的存在，也几乎不可能精确计算某个具体天体(如某个恒星)实际的长期运动。尽管如此，对于处于准稳态的恒星集体(即星系)而言，前述章节所介绍的Binney等人的工作表明：恒星集体的统计性质，例如各种分布，是可以利用“作用量-角变量”方法来可靠地获得的。我们进一步感兴趣的是，“作用量-角变量”方法是否可以应用于共振、混沌等所谓“不可积的”的恒星运动，可以应用于哪些具体场合？最近Binney[16,17]研究了把环面映射方法应用到棒旋星系的共振轨道的情况(以银河系观测数据为实例)。他在文献[16]中，分别研究了在银河系棒的这种扰动强度下，“作用量-角变量”工具对于外林德布拉德共振、共转共振、内林德布拉德共振的适用性。他发现，对于银河系被外林德布拉德共振和共转共振所俘获的恒星轨道，环面映射法能很好地拟合出；但是更靠近棒内部的内林德布拉德共振俘获的轨道，环面映射法就显得不够好了。Binney[17]更进一步把上述方法应用到太阳邻域恒星的速度空间数据中，发现在银河系的棒模型下，环面映射法所拟合的共转共振情形与观测数据一致，而外林德布拉德共振情形不被数据所支持；也就是说，基于环面映射法和微扰理论的计算，对于棒旋星系中的共振俘获现象，可以有很好的实用性了。

4.3 应用于远离平衡的星系过程：潮汐星流

当一个星团或矮星系(称为前身天体)靠近一个大星系(称为主星系)时，前身天体各个不同位置所受的引力不同，从而发生形状和速度的变化；当所处的引力场梯度足够大时，前身天体的恒星们会被主星系巨大潮汐力逐步拉扯，或多或少的恒星脱离前身天体，常形成长条状或称纤维状的结构，此即潮汐星流(详见文献[33])。

研究者发现，对于星流的这种形成过程，以及对于星流特性的刻画，最方便的方式是基于“作用量-角变量”体系[34,35]。当前身矮星系或星团的成员恒星没有被潮汐剥离，随着前身天体在主星系引力势中运动时，在主星系参考系中，这些恒星的作用量是不守恒的。但一旦某恒星被潮汐力剥离前身天体之后，就可以忽略前身天体以及星流中的其他恒星对它的引力作用，认为只是在主星系的引力势下运动。那么，自剥离时刻起，该恒星的作用量就一直保持剥离前的数值，不再随时间变化了。下面我们沿用Eyre和Binney[36]的记号和表述方式，来简介这个方法。设前身天体(小星系或星团)位于固定背景势下的规则轨道上，其作用量为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]：

什么样的引力势分布形式(必要条件)可以形成星流等技术细节，这里不作详细介绍，读者可以参考Sanders的博士论文(第5章和第6章)[20]。特别是该论文第6章详细介绍了如何应用星流数据来拟合限定主星系的引力势分布。

最后，我们要指出的是：当前身天体为很小的星系或者星团时，星流的轨迹可以粗略地(1阶近似)看作是前身天体在主星系中的轨道；因此，过去研究者常简单假定星流为一条轨道，从而拟合主星系引力势，此即文献中所说的轨道拟合法(例如文献[38])。但是，基于上面的分析，星流的轨迹并不必然与前身天体的轨道或者任何其他轨道重合(参见文献[36]的图13)；因此轨道拟合法会导致较大的主星系引力势偏差(详见文献[20]的第6章)。

4.4 结束语

本文主要目的是介绍星系动力学建模方面的新进展，即作用量计算方法，以及近年来基于作用量(或运动积分)的分布函数建模在银河系数据上的应用，另外还简介了“作用量-角变量”方法在不可积的星系现象(共振和星系棒)和非平衡过程(例如潮汐星流的形成)的应用，以期为中国在“Gaia+LAMOST”黄金时代的银河系研究做些贡献。目的之二，厘清上述方法及其应用的理论基础和脉络。在实际或理想的星系中，唯一的一类可积系统是具备Stckel形式引力势的星系；相应地，几乎所有星系建模工作的解析性基础，都是Stckel势理论。因此，我们还在本文第1章花不少篇幅从N个粒子的经典力学(哈密顿方程和刘维尔方程)出发，依次引入平均场近似，介绍无碰撞玻尔兹曼方程；然后在经典力学的正则变换以及相空间的语境下，引入运动积分、正则环面坐标即作用量和角变量；接着就介绍唯一的严格可积的(理想的)“星系”，即介绍Stckel势理论的各种解析式。介绍这些来龙去脉，期望我国学者不仅能应用Binney等国外专家发明发展的方法来分析银河系数据，还能够在理论基础和方法上同样有所创新。

附录A 计算作用量、角频率和角变量的步骤

这里介绍Sanders和Binney的程序中数值计算作用量、角频率和角变量的步骤(文献[1]中第8节及附录；文献[22]的附录A)。程序中的这些中间量的计算以及轨道分类，来自Sanders和Binney的TACT程序[23]的/tact/aa/inc/stackel aa.h中的函数定义。

(1)作用量J

作用量的计算公式为