带电测试粒子在磁化史瓦西黑洞中的混沌运动.pdf

当史瓦西黑洞周围存在渐近均匀的外部磁场时，描述带电粒子在史瓦西黑洞附近运动的哈密顿系统会变为不可积系统.类似于这样的相对论哈密顿系统不存在有显式分析解的2 部分分离形式，给显式辛算法的构建和应用带来困难。近一年以来的系列工作提出将相对论哈密顿系统分解为具有显式分析解的2 个以上分离部分形式，成功解决了许多相对论时空构建显式辛算法的难题.最近的工作回答了哈密顿系统显式可积分离数目对长期数值积分精度有何影响、哪种显式辛算法有最佳长期数值性能这两个问题，指出哈密顿有最小可积分离数目即3 部分分裂解形式并且应用于优化的4 阶分段龙格库塔显式辛算法可取得最好精度.由此选择上述数值积分方法并利用庞加莱截面、最大李雅普诺夫指数和快速李雅普诺夫指标研究在磁化史瓦西黑洞附近运动的带电粒子轨道动力学.结果显示：针对某特定的粒子能量和角动量，较小的外部磁场很难形成混沌轨道；较大的正磁场参数容易使轨道产生混沌，并且随着磁场的增大，轨道的混沌程度也随之加强；粒子能量适当变大也可以加剧混沌程度，但负磁场参数和粒子角动量变大都会减弱混沌.关键词天体力学，黑洞，磁场，混沌，算法中图分类号：P138；文献标识码：A1引言近几年来，引力波 1-2 的数十起成功探测以及M87星系中心超大质量黑洞(M87*)照片 3-4)和银河系中心黑洞(SgrA*)照片 5-6 陆续证实了爱因斯坦广义相对论引力理论关于黑洞存在的预言.黑洞物理性质和引力检验已成为研究热点问题史瓦西黑洞是爱因斯坦场方程的一个解.由于史瓦西黑洞存在4 个独立的运动常数，即粒子能量、角动量、粒子的4 速度关系和类似的卡特常数(亦称为方位角运动)7 ，所以是可积系统尽管该系统理论上存在解析解，但分析解很难用初等函数表示，仅仅可以用椭圆积分形式表达而已.大多数可观测2022-04-08收到原稿，2 0 2 2-0 5-3 1收到修改稿*国家自然科学基金项目(119 7 3 0 2 0)资助tM到的黑洞都有等离子体吸积盘，可以在黑洞外产生磁场.这样的磁场一般比较弱小，虽然对黑洞时空没有什么影响，但对荷质比(电荷与质量之比)较大的粒子产生洛伦兹力作用对粒子运动影响不可忽略 8-12 .当史瓦西时空有外部电磁场浸入时，哈密顿雅可比方程很可能会变为不可分离变量的情形，即不存在类似的卡特常数，从而相应的哈密顿系统变为不可积系统.带电粒子在磁化黑洞周围的运动于一定条件下会产生混沌现象 13-2 0 判定试验粒子轨道的混沌性方法有庞加莱截面、最大李雅普诺夫指数 2 1 和快速李雅普诺夫指标 2 2 等.数值方法是研究混沌系统最常用的方法.对于39-1(2)64卷哈密顿系统长期数值积分方法来说，保持哈密顿相流的辛积分器 2 3 自然成为优先选择，可使能量积分误差无长期增长趋势辛算法分为显式辛算法 2 3-2 4 和隐式辛算法 2 5；显隐混合算法 2 6-3 1 也常被使用显式算法比隐式方法具有计算效率优势.相对论时空对应的哈密顿一般没有显式分析解的两部分分离形式，导致显式辛算法的构建和应用有极大困难求解这些相对论时空问题，如采用辛算法的话，过去只能是隐式 2 5 或者显隐混合 2 6-3 1 的形式.当然，结合中点置换的扩大相空间显式类辛方法 14,3 2-3 7 可以适合用于求解不能变量分离的相对论时空问题.近一年以来伍教授团队在天体物理杂志等期刊上发表了系列工作 17,2 0,3 8-4 4 ,，提出将相对论时空对应的哈密顿或时间变换的哈密顿分解为具有显式分析解的2 个以上分离部分形式，成功解决了许多相对论时空构建显式辛算法的难题.这一思路还可以应用于Yoshida高阶显式辛算法 4 5-4 6 问题是哈密顿系统显式可积分离数目是否影响长期数值积分精度还有一个问题是何种显式辛算法有最佳长期数值性能作者最近发表在天体物理杂志上的一个工作 4 4 回答了这两个问题，指出哈密顿有最小可积分离数目即3 部分分裂解形式并且应用于优化的4 阶分段龙格库塔显式辛算法(Partitioned-Runge-Kutta,PRK64)47可取得最好精度.注意这种4 阶PRK方法与同阶Yoshida算法相比，包含更多额外的时间系数和子哈密顿解的组合，故优化后的PRK算法在计算上比Yoshida算法更耗时但额外增加的时间不太多，并且拥有更小的截断误差，因而这种方法值得推荐应用.本文主要目的是沿着近期工作 4 4 的思路，采用4 阶优化的PRK显式辛方法应用于具有显式分析解的3 部分哈密顿分解形式并结合庞加莱截面、最大李雅普诺夫指数和快速李雅普诺夫指标研究在磁化史瓦西黑洞附近运动的带电粒子轨道动力学.2牛物理模型本节先给出描述带电粒子在附带渐近均匀磁场的史瓦西黑洞周围运动的动力学模型，再给出适合于求解磁化史瓦西时空哈密顿系统的辛算法.天文学报2.1史瓦西时空在球坐标(t,r,)中，t表示坐标时间，r、表示3 维空间，类似于空间直角坐标系中的c、y、么.描述一个带电荷q的试验粒子在磁化史瓦西黑洞附近运动的动力学哈密顿函数方程为：T=zgl(p-qA)(pu-qAv),H2其中，q是粒子的电荷，p是由一组标准哈密顿正则方程多一兴确定的产义动量即力表示协变广义动量，A、A,表示电磁四向量势，上下标、v表示不同坐标符号，代表球坐标(t,r,)中任意坐标符号，史瓦西黑洞度规的非零逆变分量gu为：gt=1/gt=-(grr=1/grr=1这里gt、9 r、9 o o、9 s s 是史瓦西黑洞度规的协变分量.4速度iu是坐标u关于原时的导数.两个常数协变动量分量为：pt=gtut=-(1P=gosp+qAp=r?sin262Pt、P分别表示在t与方向上的广义动量，B表示磁场强度，E是粒子的恒定能量，L是粒子的恒定角动量.(4)式中电磁场势的唯一非零协变分量采用 4 9 的形式:B将(3)-(5)式代入(1)式，磁化史瓦西哈密顿可以简化为2 自由度4 维相空间系统，即：2HP-22P1(L_Brsin?e)2个22r2 sin?04期(1)212t=-E,qB2BsinH2E22(3)L.(4)(5)(6)39-264卷Pr、p e 分别表示在r与方向上的广义动量，这里=qB.(6)式实际是经过无量纲化处理后所得到的.光速c和引力常数G单位化处理：C=G=1,无量纲处理就是以黑洞质量M和粒子质量m来进行标度变换：ttM,r rM,B B/M,E mE,pr mpr,L mML,pe mMpe,q mq,H mH.其中m为试验粒子质量，M引力源质量变为几何化单位，即：M=1.哈密顿函数(6)式除含两个常数E和L外，还有哈密顿本身总是一个守恒量.对于类时时空，这个常数是：H2当史瓦西黑洞外无渐近均匀磁场时，哈密顿函数(6)式的第4 个运动积分存在，故而是严格可积的；但当存在渐近均匀磁场时，第4 个运动积分不存在，从而导致系统不可积，进而可能会产生混沌.2.2算法实现文献 3 9 将哈密顿函数(6)式分解为4 个显式可积求解部分构建显式辛算法.作者近期的工作 4 4 指出哈密顿函数(6)式还可分解为3、5个显式可积求解部分来构建显式辛算法，发现3 部分分解方法精度最高.按照这一思路，3 个部分裂解方法如下：H=Hi+H2+H3,其中每个部分的哈密顿为：2H1=1一2（12r2 sin?01P2H2721H很明显，3 个部分哈密顿函数H1、H 2、H 3 都是可积且可以写出解析解，其解析解是时间t的显式函数.求解3 部分的算子依次定义为礼H设h为积分步长，对哈密顿H提供2 个一阶辛算子：S1=HH HR,周娜英等：带电测试粒子在磁化史瓦西黑洞中的混沌运动hS2H=S1H（2）S1H*（2）3个2 阶对称乘积可以构造4 阶Yoshida显式辛积分器 4 5;S4H=S2H(cih)S2H(czh)S2H(cih)，(15)1(7)(8)E2+2L gBr2 sin24期S1H=HH HN H.(13)其中S1H*是S1H的共轭辛算子.于是，这2 个一阶辛算子的积是一个对称组合，也就是二阶显式辛方法：(14)式中：C1=1/(22 1/3)，C2=12 C1.进一步可以建立优化的4 阶分段龙格库塔PRK64显式辛算法 4 4,4 7 ;PRKg4=S1Ha12 S1Ha1 S1Ha1o S1Ha。S1Has S1Ha,S1Hag S1Ha,S1Haa S1Ha,SIHa S1Ha1:从第(n1)步到第n步，具体离散差分方案如下：S1He2=sin on-/Hs-pfa cos On-/po.m-1,(9)fi=arctan 2(e2,e1),f2=tan(en-1-fi),gH2=fi+arctan(ei+e2)hpo,n-1+fal,(10)(11)(12)hrH2=1/(e1 cos On-1+e2 sin On-1).Hl:=pa+L-(ra)sin2 Ha/2)2/(r-H2)sin2 gHla-L-(rH2)?sin2 g2/2/rH2 _E2/rH2 2)2),(16)39-364卷S1H*p:H=po,n-1+hl-cos On-1(L-(r-1Posin 0n-1/2)2)/(rm-1 sin 0n-1)-B cos n-(L-r-1sin On-1/2)/sin On-1.天文学报sin2 gH2/21)/(rH2)3 sin2 gHa-(cos 0H2 L-(rH2)2sin?gH2/2)/sin gH2.sin 0n-1/2)/(r-1 sin?On-1)-(L-rn-1 sin?n-1/2)/n-1-E2/(rn-1-2),4期f#=arctan 2(e,e),f=tan(0n-1-f),a=f+aretan(e1?+e2)hpi+f2,r*Ha=1/(ei cos On-1+eg sin On-1).式中计算出的1,2，12 列举在表1中.所谓优化算法意味着系数ai、b,之间的自由系数在5阶(优化的4 阶显式辛算法)截断误差项中系数平方和最小化与非优化算法相比，优化算法可以减小截断误差(12)-(16)式都是求解哈密顿函数(6)式或(8)式的显式辛算法.表1PRK64算法相关系数Table 1 Correlation coefficient of PRKs4 integratorsIntegratorsCoefficient3类数值分析先评估算法优劣，再挑选好的算法并结合混沌指标来研究带电粒子轨道的动力学，特别探讨3 个动力学参数对轨道混沌的影响.3.1算算法对比取时间步长为h=1以及参数为磁场参数=8.910-4，粒子能量E=0.995和角动量L=4.6.初始条件为pr=0、=/2 一旦初始半径r给定，pe初值取正值并且由(7)式确定图1为两种4 阶算法的哈密顿误差图图1(a)为两种4 阶算法积分初始半径r=11的轨道所得哈密顿误差H=H+1/2图，其中红色的曲线表示Yoshida4阶显式辛算法S4，绿色曲线表示PRK4阶显式PRK641=12=0.0792036964311962=11=0.130311410182166Q3=Q10=0.2228614958676084=g=-0.3667132680474265=8=0.324648188689706Q6=7=0.109688477876750辛算法PRK64可以看到两条曲线都不随时间的增大而增大，在长期积分中保持平稳显式辛算法S4的哈密顿误差阶数值最后稳定在10 的-9 阶到-8 阶次之间，而优化后的显式辛算法PRK64算法的哈密顿误差阶数值稳定在10 的-13 阶到-12 阶次之间，比S4算法精度高4 个数量级.图1(b)为初始半径r=88的误差图两种4 阶算法的精度非常高大轨道初始半径比小轨道初始半径精度高的原因在于前者平均轨道周期要小于后者，但积分到=10 4 后两条曲线全部上移，原因是截断误差远小于计算机的舍入误差导致的.即使误差有长期增长趋势，但精度仍可达到10 的-13 阶次左右.还要注意到优化后的显式辛算法需要使用更多的计算时间.表2 中列出两种算法的中央处理器(Central39-464卷ProcessingUnit，CPU)时间，表明优化后的PRK方法计算时间比S4算法稍长一些，但总的来说，计算时间依旧很短，最长CPU时间不超过3 min.因此，考虑到精度和效率，PRK64算法值得在实算中使用.-6-7-8-9-1050-12-13-14-15-160-10(b)r=88-11-12HV-13014-15-1650图1两种4 阶算法的哈密顿误差图。磁场参数=8.910-4，粒子能量E=0.995,角动量L=4.6.Fig.1 Hamiltonian errors for two fourth-order integrators.Magnetic parameter =8.910-4,particle energyE=0.995,angular momentum L=4.6.表2 图1中两种算法的CPU时间Table 2 CPU times for the two integrators inFig.1IntegratorsS4PRK64r=11122r=881483.2轨道动力学图1中的两个轨道对同一算法的哈密顿误差的长期表现不同是由于两轨道有不同的动力学性质.实际上，初始半径r=11对应有序轨道，即轨道是拟周期的；而初始半径r=88对应混沌轨道，即轨周娜英等：带电测试粒子在磁化史瓦西黑洞中的混沌运动(a)r=11S4PRK.412345671gTPRKWMV124期道演化对初始条件微小改变有指数式的敏感依赖性.轨道的有序和混沌动力学属性可借助轨道的3 维空间构型和平面投影图来粗略地观察选取初始轨道半径r=25，粒子能量E=0.995和角动量L=4图2(a)-(c)画出了3 个不同磁场参数的3 维空间轨道构型，图中红色为轨道在赤道面0=元/2 上的投影在图2（a)中，=0显示试验粒子在史瓦西黑洞时空中的运动轨道限于赤道平面上，并且围绕黑洞做圆轨道周期运动在图2(b)中，=1.210-4对应3 维空间轨道不在一个平面上，其投影是一组相交环曲线在图2(c)中,=0.01对应3 维空间轨道也不在一个平面上，其投影由一些互相交错环曲线组成.图2(c)与图2 (b)相比，轨迹更杂乱、更没有规则.这两个运动轨道是混沌还是有序难以判断.轨道有序或混沌的准确判定还需要借助其他方法.3.2.1混沌指标判断带电粒子在史瓦西黑洞周围运动的轨道是否混沌有许多方法.例如庞加莱截面法、李雅普345671gt1532*24诺夫指数和快速李雅普诺夫指标等.庞加莱截面适合研究2 自由度的4 维保守系统，可用于判断系统的运动状态.若以0=元/2 为曲面,上方为0 元/2，下方为0 元/2，可以按照线性插值来求曲面上的点.当曲面上只得到几个点或者这些点构成一封闭曲线时，系统做拟周期运动；当截面上得到的点是杂乱无章随机分布在一个区域时，系统做混沌运动.李雅普诺夫指数(Lyapunov indicator)是衡量两邻近轨道随时间平均指数分离比的指标，能够反映轨道混沌的强度.最大李雅普诺夫指数判断混沌较为准确，计算方法主要有变分法和两粒子法 2 1.本文使用变分法计算：(17)T80 T13(0)1式中：S(T)和S(0)分别表示T时刻与初始时刻的切向量.入值趋于稳定到一个正值，说明有界轨道处在混沌状态；入值趋于0，说明有序轨道处于有序状态.39-564卷入是定义的一个值，用来表示轨道的混沌状态.快速李雅普诺夫指标(Fast Lyapunov indica-tor，FLI)与李雅普诺夫指数相比，可以更快地看到轨道的混沌，得到混沌指数，其定义为 50 ：FLI=lg S()I.若FLI随时间增长很快，则有界轨道混沌；若FLI随时间增长非常慢，则有界轨道有序.(a)=00=元/2(b)=1.210-49=元/2(c)=0.019=元/2图2 3 维轨道构型及其在赤道面0=元/2 上的投影Fig.2 Three-dimensional trajectories and their projectionsat the plane=/23.2.2粒子运动轨道的动力学图3 为赤道面=元/2 上的庞加莱截面.图3(a)为图2(b)中粒子轨道在赤道面=/2上的庞加莱截面图，截面上的轨道是闭合的环轨道，说明图2(b)中粒子做有序的拟周期运动.图3(b)为图2(c)中粒子轨迹的庞加莱截面图，截面上的轨道是随机离散的点分布在一个面区域上，表明粒子做混天文学报沌运动显然，相比3 维空间构型和平面投影图，庞加莱截面区分轨道有序或混沌要更清楚、准确带电试验粒子在史瓦西黑洞周围运动时，磁场破坏了第4 个运动积分的存在，导致系统不可积，是粒子运动混沌的根本原因.(18)0.20.10.0-0.1-0.200.40.30.20.10.0-0.1-0.2-0.3-0.4图3 粒子轨道在赤道面=元/2 上的庞加莱截面.图(a)：图2(b)中粒子的轨道.图(b)：图2(c)中粒子的轨道.Fig.3 Poincar cross-section of the particle orbit at theplane =/2.Panel(a):the particle trajectory is the sameas that in Fig.2(b).Panel(b):the particle trajectory is the图4 为进一步改变磁场参数值时粒子轨道演化情况.图4(a)-(c)为粒子轨道运动的庞加莱截面图，画了3 条轨道：轨道1（红色轨道r=40)，轨道2(蓝色轨道r=65)和轨道3（黑色轨道r=88)，其余轨道为绿色轨道r=11，紫色轨道r=20.固定粒子能量E=0.995与角动量L=4.6不变，但磁场参数在3 个子图中依次为=1.210-4、=5.6 10-4、=8.9 10-4.图4(a)中=1.210-4,所有轨道都在做有序的拟周期运动，因此截面上的相4期(a)=1.210-450100(b)=0.01 1020same as that in Fig.2(c).15030402005039-664卷轨道都是闭合的环轨道.图4(d)中的最大李雅普诺夫指数也同样揭示图4(a)中所有轨道都是有序的，因为当积分时间达到10 7 后，这些李雅普诺夫指数还不能达到稳定值，都有趋于0 的趋势，表明轨道的有序性.图4（g)中所有轨道的积分到最后FLI值都小于3，没有随时间指数增长的趋势，属于有序轨道的特征.图4(b)中=5.610-4，粒子轨道2 呈现弱混沌性，轨道1依旧做拟周期运动，轨道3 亦如此.图4 (e)、(h)中的李雅普诺夫指数和快速李雅普诺夫指标分别显示轨道1、3 为有序轨道，而轨道2 在图4(e)中的李雅普诺夫指数趋于一个稳定的正值,周娜英等：带电测试粒子在磁化史瓦西黑洞中的混沌运动4期表明此轨道为混沌轨道.在图4(h)中轨道2 快速李雅普诺夫指标随时间指数增长，说明轨道2 是混沌的.因此，3 种方法都判定图4(b)中轨道1、3 为有序轨道，轨道2 为混沌轨道.图4（c)中轨道1与轨道3 为混沌轨道，轨道2 形成5个有序岛屿.轨道1、3 的混沌性与轨道2 的有序性质也被图4 (f)、(i)中的最大李雅普诺夫指数和快速李雅普诺夫指标所证实.故得出结论：随着正磁场参数的增大，亦即2增大，粒子运动轨道从有序走上混沌，并且混沌程度也随之增强.这是因为带电粒子受到来自磁场的洛伦兹力表现为引力效果变大的缘故.0.20.10.0-0.1-0.20.20.10.0-0.1-0.20.20.10.0-0.1-0.2图4 轨道的庞加莱截面、最大李雅普诺夫指数、FLI随磁参数的变化.图(a)-(c)：3 条研究轨道半径(r=40,r=65,r=88)的庞加莱截面图，粒子能量E=0.995，角动量L=4.6;图(d)-(f)为图(a)-(c)的最大李雅普诺夫指数图；图(g)-(i)为图(a)-(c)的FLI.Fig.4 The variation of the Poincar cross-section,maximum Lyapunov indicator,and fast Lyapunov indicator with magneticparameters.Panels(a)-(c):Poincar cross-sections of three studied orbital radi(r=40,r=65,r=88).Particle energyE=0.995,and angular momentum L=4.6;panels(d)-(f)are the maximum Lyapunov indicator graphs of panels(a)-(c),图5中磁参数为负值时，粒子轨道的庞加莱截面随磁参数的变化考虑5条轨道，初始半径依次为r=11、r =2 0、r =4 0、r =6 5、r =8 8,颜色意义同图4.图5(a)中有两条轨道处于强混沌状态；图5(b)中也有两条轨道处于混沌状态，但属于弱混沌轨道；图5（c)中只出现一条混沌轨道，且处于弱混沌状态；图5(d)中没有混沌轨道出现，所有轨道-2(a)=1.2x10-4-33-4-5-2(b)=5.610-4-32730306090120150180210respectively;panels(g)-(i)are FLI of panels(a)-(c),respectively.(d)=1.2x10-4(e)=5.6x10-43-5(c)=8.9x10-4-2-3-4-516(g)=1.210-454320807(h)=5.6x10465432(f)=8.9x10-7(i)=8.9x1046543322341gt都是拟周期有序轨道.因此，随着负磁场参数的绝对值减小，粒子的轨道从混沌状态变为有序状态.也就是说，当3指示有界轨道的混沌性，而FLI不超过3 表明有界轨道的有序性图8(a)中初始轨道半径r=11，E小于0.9947为有序轨道，E大于0.9 9 57 时混沌开始；图中点3 所示的庞加莱截面为弱混沌轨道，点1与点2 为有序轨道.图8(b)中初始轨道半径r=88，比E=0.9928大的能量诱发混沌；图中点1所示的庞加莱周娜英等：带电测试粒子在磁化史瓦西黑洞中的混沌运动4期截面为强混沌轨道，点2 为弱混沌轨道，点3 表现为有序轨道.从图8 (a)、(b)中得出结论也能与前面图5的对应，即粒子能量E适当增大可以增强轨道混沌程度.图8(c)、(d)绘制了固定磁场参数=8.910-4与能量E=0.995，而角动量L从3.8 到5的FLI图图8（c)中初始轨道半径r=11时，粒子轨道在L4.11时为混沌的，而在大于4.11后轨道变得有序；图中点1所示的庞加莱截面为混沌轨道，点2 为有序轨道.图8(d)中初始轨道半径为r=88,轨道在大于L=6.35后开始变得有序；图中点1与点2 所示的庞加莱截面皆为弱混沌轨道.从图8（c)、(d)中得出结论也能与前面图6 的对应，即粒子角动量L变大可以减小轨道混沌程度.图8(e)、(f)绘制了固定粒子角动量L=4.6与能量E=0.995，磁场参数从0 跨度到1.6 10-3 的FLI图.诱发混沌的临界值依次为9.0 110-4、6.9 2 10-4.从图8(e)、(f)中得出结论与前面图3、4 的对应，即随着正磁场参数的增加会增强轨道混沌程度.0.20.10.0-0.1-0.20.30.20.10.0-0.1-0.2-0.3Fig.6 The variation of the Poincare cross-section with the particle energy E,particle angular momentum L=4.6.0.2(a)=8.9104E=0.990204002040600.3(b)=8.9x10-40.20.1E=0.9920.03-0.1-0.26080100(d)=1.5x10-3E=0.99380100图6 改变粒子能量E的庞加莱截面图，粒子角动量L=4.6.(c)=8.910-4E=0.9950.10.0-0.1-0.260100120020400.30.20.10.0-0.1-0.2-0.302040 60 8010012039-9-0.380140(e)=1.5x10-3E=0.992020 4060801001201401600.30.20.10.0-0.1-0.2-0.3020406080100120140(f)=1.5x10-3E=0.99564卷0.20.10.0-0.1-0.2天文学报(a)L=4.21-23-4-54期9:87(d)L=4.21654IT3321(g)L=4.21 3120.20.10.0-0.1-0.20.20.10.0-0.1-0.2图7轨道的庞加莱截面、最大李雅普诺夫指数、FLI随试验粒子角动量L的变化.磁参数=8.910-4，粒子能量E=0.995.Fig.7 The variation of the Poincar cross-section,maximum Lyapunov indicator and fast Lyapunov indicator with the angularmomentum of experimental particle L.Magnetic parameter =8.910-4 and particle energy E=0.995.90(a)r=)8070605040302010FLI=390(b)r=888070605040302010FLI=3 100.990图8 快速李雅普诺夫指标FLI与参数的依赖关系.每个FLI值取自积分时间=10 7 后的结果，有序与混沌的FLI临界值为3.图(a)、(b)中L=4.6、=8.9 10-4.图(c)、(d)中E=0.995、=8.9 10-4.图(e)、(f)中L=4.6、E=0.9 9 5.Fig.8 The dependence between fast Lyapunov indicator FLIs and parameters.Each of the FLIs are obtained after theintegration T=107,and the FLI value from order to chaos is 3.In panels(a)and(b),L=4.6 and =8.9x10-4;in panels(c)8(b)L=4.6-2-3-4-5(c)L=4.88-2-30-4-503060 90120 150 1800.9920.9940.9960.9981.000Eand(d),E=0.995 and =8.910-4;in panels(e)and(f),L=4.6 and E=0.995.(e)L=4.6(f)L=4.881260(c)r504030201068(d)r=88504030201003.876IT54321087(i)L=4.8865IT43203451gt224.04.2(h)L=4.667FLI=3FLI=34.44.6L2.53.03.54.04.55.0 5.56.01gT90(e)r=118070605040302010FLI=39080(0)r=8870605040302010FLI=304.85.00246810121416(10-)39-1064卷图9 为2 维参数空间对应的FLI分布.图9(a)是2维参数空间E和L对应的FLI分布初始轨道半径r=11时，E值取值接近0.9 9 57 左右开始发生混沌，并且混沌程度随着E值的增大而增大，增大到一定值后混沌开始减弱，呈现在图中颜色由蓝色依次过渡到红色，再由红色转换到蓝色，而混沌区域大致在E=0.9957到E=1之间，这个结果与图8(a)的一维图是一致的.还可以看出增大粒子角动量L可以减弱轨道的混沌程度，这与图8（c)的一维图，即E=0.995时显示的情况相同.图9（c)为初始轨道半径r=88时扫描参数E和L的2 维分布图情形，FLI1.000-88.0(a)r=11-80.00.99972.064.00.998-56.0-48.00.997-40.0-32.00.99624.0-16.00.995-8.000.994-0.004.04.14.24.34.44.54.64.74.84.95.0LFLI1.00088.0(c)r=88-80.00.999一0.9980.9970.9960.9950.9940.9930.992.4.0图9 2 维参数空间对应的FLI分布。积分时间均为T=107；图(a)、(c)：=8.910-4；图(b)、(d)：L=4.6.Fig.9 The corresponding FLI distribution for two-dimensional parameter spaces.The integral time is =107 for all plots;周娜英等：带电测试粒子在磁化史瓦西黑洞中的混沌运动78910111213(104)FLI1.00088.0(d)r=88-80.00.99972.00.998-64.0-56.0-48.0-40.0-32.0-24.0-16.0-8.00-0.004.24.4Lpanels(a),(c):=8.910-4;panels(b),(d):L=4.6.4期更清晰地展示了粒子角动量增大减弱了轨道的混沌程度，这与图8(b)、(d)显示的结果相同.扫描图9(b)中参数E和的2 维分布图显示在初始轨道半径r=11时，能量变大，故粒子受到的引力变大会加剧混沌程度；但能量大到接近1时，轨道有不稳定趋势，相应减弱混沌程度；而当值增大时，轨道由有序变为混沌，且FLI值逐渐变大，即混沌程度逐渐加强.图9(d)为初始轨道半径r=88时扫描参数E和的2 维分布图，可以得到与图9(b)相同的结论.图9 (b)、(d)与图8(e)、(f)扫描的一维图显示结果相同.1.000(b)r=110.9990.998E0.9970.9960.9950.99434560.997E0.9960.9950.9940.9930.9924.64.8FLI-88.080.0-72.064.0-56.048.040.032.024.0-16.0-8.000.00-72.0-64.0-56.0-48.0-40.0-32.0-24.0-16.0-8.000.005.02345678910111213(10-4)39-1164卷4总结与展望辛算法因具有保结构性质，可为天体长期定性演化研究提供可靠的数值结果.在计算效率方面，显式辛算法比隐式辛算法好，磁化史瓦西时空的哈密顿函数可以分解为具有显式可积解的3、4、5部分，由此可以进行显式辛算法的构建与应用优化的4 阶分段龙格库塔PRK显式辛算法结合哈密顿3部分分解方法能够取得最佳数值效果.利用这样的数值积分方法来研究附带渐近均匀的磁场下史瓦西黑洞周围带电粒子轨道动力学.本工作主要探讨磁场参数、粒子能量E和角动量L对粒子轨道动力学行为的影响.磁场是哈密顿系统不可积性和产生混沌的关键因素.磁场参数绝对值的变大可以使得粒子运动轨道由有序变为混沌，增强轨道混沌程度粒子能量增加亦可加强混沌程度，但粒子角动量增大反而对混沌有抑制作用.本工作所建立的算法适用于许多复杂的相对论时空的哈密顿函数或时间变换哈密顿函数，可以用来探索黑洞轨道动力学性质，还可以用于黑洞阴影数值模拟研究.致谢感谢伍歆教授对本文的帮助和建议，也感谢参与本文评审的各位专家对文章提出的诚挚建议，使得文章的质量有了显著的提高.参考文献1 Abbott B P,Abbott R,Abbott T D,et al.PhRvL,2016,116:0611022 Abbott B P,Abbott TD,Abraham S,et al.ApJ,2020,900:L133 The Event Horizon Telescope Collaboration.ApJ,2019,875:L14 The Event Horizon Telescope Collaboration.ApJ,2019,875:L65 The Event Horizon Telescope Collaboration.ApJ,2022,930:L126 The Event Horizon Telescope Collaboration.ApJ,2022,930:L137 Carter B.Physical Review Journals Archive,1968,174:15598 Ernst F J.JMP,1976,17:549 Panis R,Kolos M,Stuchlik Z.EPJC,2019,79:47910 Frolov A P,Shoom A A.PhRvD,2010,82:08403411 Karas V,Vokrouhlicky D.GReGr,1992,24:729天文学报12 Kopacek O,Karas V,Kovar J,et al.ApJ,2010,722:124013 Stuchlik Z,Kolos M.EPJC,2016,76:3214 Li D,Wu X.The European Physical Journal Plus,2019,134:9615 Vieira W M,Letelier P S.ApJ,1999,513:38316Voorhees B H.PhRvD,1970,2:211917 Hu A R,Huang G Q.The European Physical JournalPlus,2021,136:121018 Yi M,Wu X.PhyS,2020,98:08500819 Al Zahrani A M,Frolov V P,Shoom A A.PhRvD,2013,87:08404320 Zhang H X,Zhou N Y,Liu W F,et al.Univ,2021,7:48821Wu X,Huang T Y.PhLA,2003,313:7722Wu X,Huang T Y,Zhang H.PhRvD,2006,74:08300123 Ruth R D.ITNS,1983,30:266924 Zhang L N,Wu X,Liang E W.MEdRJ,2021,9:271825 Feng K,Qin M Z.Symplectic Geometric Algorithms forHamiltonian Systems.Hangzhou:Zhejiang Science andTechnology Publishing House,2009:187-21126 Zhong SY,Wu X,LiuS Q,et al.PhRvD,2010,82:12404027 Wu X,Zhong S Y.GReGr,2011,43:218528 Zhong S Y,Wu X.AcPSn,2011,60:09040229 Zhong S Y,Liu S.AcPSn,2012,61:12040130 Mei L J,Wu X,Liu F Y.EPJC,2013,73:241331 NMei L J,Ju M J,Wu X,et al.MNRAS,2013,435:224632 Liu L,Wu X,Huang G Q,et al.MNRAS,2016,459:196833