第42卷第7期 大学物理 Vol.42 No.7 2023年7月 COLLEGE PHYSICS Jul.2023 带电粒子在地磁场中的运动及 Mathematica 数值模拟 董顺成，郭芳侠 (陕西师范大学物理学与信息技术学院，陕西西安710119) 摘要：基于单粒子轨道模型和地磁场偶极子模型，考虑相对论效应，对近地球区域磁场中运动的带电粒子轨迹使用Math- ematica软件中六阶龙格—库塔算法进行数值计算和模拟，并对极光现象的产生进行了解释，同时讨论了带电粒子在地磁场中 运动的引导中心近似.结果表明：1)从地球北极方向观察，被地球磁场捕获的质子沿顺时针方向漂移，电子沿逆时针方向漂 移；2)粒子各个分运动的运动周期数值模拟结果与文献中理论值非常吻合；3)从(4R,0,0)入射的粒子投掷角小于7.38°时， 带电粒子将会与地球表面大气层碰撞而沉降，存在产生极光现象的可能.大于7.38°时，粒子将会被束缚在地磁场中，形成辐射 带；4)其他条件相同时，带电粒子投掷点距离地球越远，其漂移速度越大；投掷角越大，其漂移速度也越大；5)对于能量较低 的粒子，一阶近似下引导中心轨迹能很好地代表粒子实际运动轨迹. 关键词：地磁场；粒子运动；数值模拟 中图分类号：0441文献标识码：A 文章编号：1000-0712(2023)07-0053-08 【DOI】10.16854/j.cnki.1000-0712.220290 在近地空间中运动着能量高达GeV的宇宙射 在z轴投影上会产生较大误差.代国红[6]和钟海 线粒子，还有能量在几个keV到几百keV的带电粒 坚[7]等各自利用Matlab软件使用四阶龙格一库塔 子1，它们都会被地磁场控制，束缚在地球附近，形 算法对带电粒子在地磁场中的运动进行数值计算， 成辐射带，这会对我们的现代技术和日常生活产生 将结果可视化，但对粒子的弹射、漂移运动及投掷角 很大的影响[2] 问题的讨论过于简略，文献[7]中得出的投掷角取 20世纪50年代，Stomer[3]计算了偶极磁场中带 值范围值得商榷.此外没有研究带电粒子投掷角大 电粒子运动的轨道，在计算过程中发现地球周围存 小对漂移运动的影响. 在一个带电粒子捕获区，其辐射带被美国探索者1 本文基于单粒子轨道法和地磁场偶极子模型， 号卫星探测到，证明地球磁场可以捕获带电粒子并 使用Mathematica 软件六阶龙格一库塔算法对带电 在地球周围形成范艾伦辐射带. 粒子在地磁场中的运动进行较精确的数值计算和模 徐荣栏在《磁层粒子动力学》1中系统介绍了 拟，与四阶龙格一库塔算法相比，可以保持E/E 带电粒子在地球磁场中的运动轨道，并在此基础上 基本不变，即保证系统的能量是守恒的，在研究弹射 解释辐射带、极光等现象.但受计算机条件限制，数 运动时可以保证磁镜点对称，减小了数值计算的 值模拟并不太多.薛林首次用Mathematica软件对 误差[8] 带电粒子运动微分方程组进行数值求解，并绘制粒 探究影响带电粒子在地磁场中运动的各种因 子轨迹曲线，但在进行数值模拟时，没有考虑高能带 素.同时结合数值计算结果解释极光现象的产生条 电粒子的相对论效应，误差较大.黄朝艳基于磁层 件.引入引导中心近似，通过数值计算实例对比得出 粒子动力学理论，根据单粒子方法和磁偶极模型，详 引导中心近似的适用条件，对今后使用引导中心近 细讨论了带电粒子的运动特征以及适合不同能量的 似具有重要的参考价值. 带电粒子的计算方法.Ozturk[2通过Phyon 和Matlab 1物理模型 软件利用四阶龙格一库塔算法求解微分方程组，研 究了偶极磁场中的粒子运动，同时推导引导中心近 1.1地球磁场模型 似方程，画出其轨迹曲线，但其做出的引导中心轨迹 研究地球附近磁场时，在近地空间6R。(地球半 收稿日期：2022-06-09；修回日期：2022-09-15 作者简介：董顺成(2000—)，男，河南驻马店人，陕西师范大学物理学与信息技术学院2018级本科生. 通信作者：郭芳侠，E-mail:guofangxia@snnu.edu.cn 大学物理 (C)1994-2023 China Academic Journal Electr dittpou edxwh bntlsedu on ww.cnki.net
第 ４２ 卷第 ７ 期 大 学 物 理 Ｖｏｌ．４２ Ｎｏ．７ ２０２３ 年 ７ 月 ＣＯＬＬＥＧＥ ＰＨＹＳＩＣＳ Ｊｕｌ．２０２３ 收稿日期：２０２２－０６－０９；修回日期：２０２２－０９－１５ 作者简介：董顺成（２０００—），男，河南驻马店人，陕西师范大学物理学与信息技术学院 ２０１８ 级本科生． 通信作者：郭芳侠，Ｅｍａｉｌ：ｇｕｏｆａｎｇｘｉａ＠ ｓｎｎｕ．ｅｄｕ．ｃｎ 带电粒子在地磁场中的运动及 Ｍａｔｈｅｍａｔｉｃａ 数值模拟 董顺成，郭芳侠 （陕西师范大学 物理学与信息技术学院，陕西 西安 ７１０１１９） 摘要：基于单粒子轨道模型和地磁场偶极子模型，考虑相对论效应，对近地球区域磁场中运动的带电粒子轨迹使用 Ｍａｔｈ ｅｍａｔｉｃａ 软件中六阶龙格—库塔算法进行数值计算和模拟，并对极光现象的产生进行了解释，同时讨论了带电粒子在地磁场中 运动的引导中心近似．结果表明：１）从地球北极方向观察，被地球磁场捕获的质子沿顺时针方向漂移，电子沿逆时针方向漂 移；２）粒子各个分运动的运动周期数值模拟结果与文献中理论值非常吻合；３）从（４ Ｒｅ，０，０）入射的粒子投掷角小于７．３８°时， 带电粒子将会与地球表面大气层碰撞而沉降，存在产生极光现象的可能．大于７．３８°时，粒子将会被束缚在地磁场中，形成辐射 带；４）其他条件相同时，带电粒子投掷点距离地球越远，其漂移速度越大；投掷角越大，其漂移速度也越大；５）对于能量较低 的粒子，一阶近似下引导中心轨迹能很好地代表粒子实际运动轨迹． 关键词：地磁场；粒子运动；数值模拟 中图分类号：Ｏ ４４１ 文献标识码：Ａ 文章编号：１０００ ０７１２（２０２３）０７ ００５３ ０８ 【ＤＯＩ】１０．１６８５４ ／ ｊ．ｃｎｋｉ．１０００ ０７１２．２２０２９０ 在近地空间中运动着能量高达 ＧｅＶ 的宇宙射 线粒子，还有能量在几个 ｋｅＶ 到几百 ｋｅＶ 的带电粒 子［１］ ，它们都会被地磁场控制，束缚在地球附近，形 成辐射带，这会对我们的现代技术和日常生活产生 很大的影响［２］ ． ２０ 世纪 ５０ 年代，Ｓｔｏｍ¨ ｅｒ［３］ 计算了偶极磁场中带 电粒子运动的轨道，在计算过程中发现地球周围存 在一个带电粒子捕获区，其辐射带被美国探索者 １ 号卫星探测到，证明地球磁场可以捕获带电粒子并 在地球周围形成范艾伦辐射带． 徐荣栏在《磁层粒子动力学》［１］中系统介绍了 带电粒子在地球磁场中的运动轨道，并在此基础上 解释辐射带、极光等现象．但受计算机条件限制，数 值模拟并不太多．薛林［４］ 首次用 Ｍａｔｈｅｍａｔｉｃａ 软件对 带电粒子运动微分方程组进行数值求解，并绘制粒 子轨迹曲线，但在进行数值模拟时，没有考虑高能带 电粒子的相对论效应，误差较大．黄朝艳［５］ 基于磁层 粒子动力学理论，根据单粒子方法和磁偶极模型，详 细讨论了带电粒子的运动特征以及适合不同能量的 带电粒子的计算方法．Ｏ¨ ｚｔ ｕ¨ ｒｋ［２］ 通过 Ｐｈｙｏｎ 和 Ｍａｔｌａｂ 软件利用四阶龙格—库塔算法求解微分方程组，研 究了偶极磁场中的粒子运动，同时推导引导中心近 似方程，画出其轨迹曲线，但其做出的引导中心轨迹 在 ｚ 轴投影上会产生较大误差．代国红［６］和钟海 坚［７］ 等各自利用 Ｍａｔｌａｂ 软件使用四阶龙格—库塔 算法对带电粒子在地磁场中的运动进行数值计算， 将结果可视化，但对粒子的弹射、漂移运动及投掷角 问题的讨论过于简略，文献［７］中得出的投掷角取 值范围值得商榷．此外没有研究带电粒子投掷角大 小对漂移运动的影响． 本文基于单粒子轨道法和地磁场偶极子模型， 使用 Ｍａｔｈｅｍａｔｉｃａ 软件六阶龙格—库塔算法对带电 粒子在地磁场中的运动进行较精确的数值计算和模 拟，与四阶龙格—库塔算法相比，可以保持 Ｅｋ ／ Ｅｋ０ 基本不变，即保证系统的能量是守恒的，在研究弹射 运动时可以保证磁镜点对称，减小了数值计算的 误差［８］ ． 探究影响带电粒子在地磁场中运动的各种因 素．同时结合数值计算结果解释极光现象的产生条 件．引入引导中心近似，通过数值计算实例对比得出 引导中心近似的适用条件，对今后使用引导中心近 似具有重要的参考价值． １ 物理模型 １．１ 地球磁场模型 研究地球附近磁场时，在近地空间 ６Ｒｅ（地球半
54大学物理第42卷径R。=6.38×10°m）以内，地球附近磁场可近似为偶2带电粒子在地磁场中运动的轨迹极磁场9)，在笛卡尔坐标系中，以磁极方向为z轴，2.1质子在地磁场中的运动地球磁场的磁感应强度的表达式为B.R.以太阳风中能量为10MeV的质子为例其电B(x,y,z)=[3xze+量g=1.60218×10-19C，m=1.67262×10-27kg.设粒子r起始位置为（4R00），投挪角为30°（入射速度与3yze,+(22 -x*-y')e.](1)磁场方向夹角），利用Mathematica求解微分方程组地磁场关于z轴具有轴对称性，x0z平面内地磁场磁（3），运行60s后，得到粒子运动轨迹如图2所示，感线分布图如图1所示Y=/R7oWRx/R0x/R(a)0—60s质子运动轨迹图1地磁场在x0z平面磁感线分布图1.2粒子运动模型单粒子轨道模型是把等离子体当做独立的粒子系统，忽略粒子间的相互作用，同时不考虑带电粒子的运动对外界电磁场的影响【0],来讨论粒子在电磁场中的运动1.3粒子运动方程在近地空间中，假设能量为E、带电有量为q、质量为m的粒子在电场强度为E和磁感应强度为B的场中运动，考虑相对论效应，粒子的运动服从牛x/R顿-洛伦兹方程，即d(ymv)(b)0—60质子轨迹（顶视）q(E+uxB)(2)dt图2质子在地磁场中运动轨迹其中为相对论因子从地球北极方向看，质子沿顺时针方向飘移.粒在近地空间辐射带中，通常典型的辐射带粒子子的运动可分为3种运动的合成：能量高于105eV，此时地磁场占主导，因此可以不考1）回旋运动：带电粒子有垂直磁感线方向的速虑地球附近电场的影响[2,5]，可取E=0.代入式（2)度，产生围绕磁感线的圆周运动；得到分量式：2）弹射运动：偶极场类似于磁镜装置，带电粒ymx=q(yB,-zB,),子会在南北方向（磁镜点处）产生周期性振荡运动；(3)ymj=q(zB,-xB,),3）漂移运动：地磁场在地球径向存在磁场梯ym2=q(xB,-yB)度，带电粒子会发生沿地球纬度方向的漂移运动[]利用Mathematica软件Plot指令绘制出一定时间内质子的x(）、(t）和z()图像，如图3所示大海靓理(C)1994-2023China Academic JoumalElectrhttp/dxwAbit/:eduhonwww.cnki.net
５４ 大 学 物 理 第 ４２ 卷 径Ｒｅ ＝ ６．３８×１０６ ｍ）以内，地球附近磁场可近似为偶 极磁场［９］ ，在笛卡尔坐标系中，以磁极方向为 ｚ 轴， 地球磁场的磁感应强度的表达式为 Ｂ（ｘ，ｙ，ｚ）＝ －Ｂ０Ｒ３ ｅ ｒ ５ ［３ｘｚｅｘ ＋ ３ｙｚｅｙ ＋（２ｚ ２ －ｘ２ －ｙ２ ）ｅｚ］ （１） 地磁场关于 ｚ 轴具有轴对称性，ｘＯｚ 平面内地磁场磁 感线分布图如图 １ 所示． 2 2 1 1 0 0 1 2 1 2 x/Re y/Re 图 １ 地磁场在 ｘＯｚ 平面磁感线分布图 １．２ 粒子运动模型 单粒子轨道模型是把等离子体当做独立的粒子 系统，忽略粒子间的相互作用，同时不考虑带电粒子 的运动对外界电磁场的影响［１０］ ，来讨论粒子在电磁 场中的运动． １．３ 粒子运动方程 在近地空间中，假设能量为 Ｅｋ 、带电有量为 ｑ、 质量为 ｍ 的粒子在电场强度为 Ｅ 和磁感应强度为 Ｂ 的场中运动，考虑相对论效应，粒子的运动服从牛 顿－洛伦兹方程，即 ｄ（γｍｖ） ｄｔ ＝ ｑ（Ｅ＋ｖ×Ｂ） （２） 其中 γ 为相对论因子． 在近地空间辐射带中，通常典型的辐射带粒子 能量高于１０５ ｅＶ，此时地磁场占主导，因此可以不考 虑地球附近电场的影响［２，５］ ，可取 Ｅ ＝ ０．代入式（２） 得到分量式： γｍ ｘ¨ ＝ ｑ（ｙ ·Ｂｚ － ｚ ·Ｂｙ）， γｍ ｙ¨ ＝ ｑ（ｚ ·Ｂｘ －ｘ ·Ｂｚ）， γｍ ｚ¨ ＝ ｑ ｘ·Ｂｙ －ｙ ·Ｂｘ r i p （３） ２ 带电粒子在地磁场中运动的轨迹 ２．１ 质子在地磁场中的运动 以太阳风中能量为 １０ ＭｅＶ 的质子为例，其电 量 ｑ ＝ １．６０２ １８×１０ － １９Ｃ，ｍ ＝ １．６７２ ６２×１０ － ２７ ｋｇ．设粒子 起始位置为（４Ｒｅ，０，０），投掷角为 ３０°（入射速度与 磁场方向夹角），利用 Ｍａｔｈｅｍａｔｉｃａ 求解微分方程组 （３），运行 ６０ ｓ 后，得到粒子运动轨迹如图 ２ 所示． Copy ropy iopy g h t h g 5 4 3 g h t h g 3 4 55 4 3 g h t h g 3 4 5 ０—６０ ｓ 质子运动轨迹 x/Re z/Re y/Re 2 0 1 1 2 5 5 4 3 2 1 0 1 2 3 4 5 4 3 2 1 0 1 2 3 4 5 ０—６０ ｓ 质子轨迹（顶视） 图 ２ 质子在地磁场中运动轨迹 从地球北极方向看，质子沿顺时针方向飘移．粒 子的运动可分为 ３ 种运动的合成： １）回旋运动：带电粒子有垂直磁感线方向的速 度，产生围绕磁感线的圆周运动； ２）弹射运动：偶极场类似于磁镜装置，带电粒 子会在南北方向（磁镜点处）产生周期性振荡运动； ３）漂移运动：地磁场在地球径向存在磁场梯 度，带电粒子会发生沿地球纬度方向的漂移运动［６］ ． 利用 Ｍａｔｈｅｍａｔｉｃａ 软件 Ｐｌｏｔ 指令绘制出一定时 间内质子的 ｘ（ｔ）、ｙ（ｔ）和 ｚ（ｔ）图像，如图 ３ 所示．
第7期55董顺成，等：带电粒子在地磁场中的运动及Mathematica数值模拟2×1071×10号-1x107z/R-2x10/WR4020600t/sx(t)2(0)(0)x/R.图30—78s质子x(）（）和z(）曲线(a)0—25s电子运动轨迹曲线图3显示x（）、y（）和z（）曲线都具有周期性，x（）和y(）曲线周期对应于粒子漂移运动周期z（）曲线周期对应于弹射运动周期，弹射运动周期远远小于漂移周期根据PankajK.Soni等["]推导的粒子周期理论公式，粒子弹射运动周期：Ro)cT, = 0.117[1-0.463 5(sin αe)]（4）R.J粒子漂移运动周期为2mqBR1[(5)(sinT4=R.m?其中α为带电粒子从赤道面入射时的投掷角，0°≤x/R.α≤180°,R。为带电粒子人射处距离地球球心的距(b）0—100s电子轨迹（顶视）离.代人数据得：T=2.34056s，T。=77.96360s.计图4电子在地磁场中运动轨迹算图3中z（）曲线与横轴的第三个交点坐标（2.34145，0），即质子弹射运动周期为2.34145s.1x10/y（0）曲线与横轴t的第3个交点坐标为（77.963580），即质子漂移运动周期为77.96358s.可见数值5×106计算的结果与粒子周期理论值非常吻合u/z2.2电子在地磁场中的运动设电子在地磁场中运动的初始条件与2.1中质5X1子相同.起始位置（4R。，0，0），投掷角30°.以太阳风中能量为10MeV的电子为例，电有量g=-1.60218×10-19c，m=9.10956×10-kg.电子在地01234t/s磁场中的运动轨迹如图4所示图50—5电子（t）曲线从北半球方向观察，电子漂移运动沿逆时针方向，与质子的漂移运动方向刚好相反.电子荷质比非常大（约是质子的1000倍），其回旋运动半径很小，3带电粒子在地磁场中运动分析因此图4中几乎看不到回旋运动.再观察5s内电子由于电子回旋运动不明显，并且轨迹曲线过于z(t)图像，如图5.紧密，因此下面主要以质子为例进行分析讨论可知电子在z方向的弹射运动周期非常小，因此图4中电子轨迹曲线非常紧密大学如理(C)1994-2023China Academic JoumalElectrhttp/dxwAbt/:eduh@nyww.cnki.net
第 ７ 期 董顺成，等：带电粒子在地磁场中的运动及 Ｍａｔｈｅｍａｔｉｃａ 数值模拟 ５５ C C oC pC yC oriCg iriCg iriCg oriCg t/s L/m x(t) y(t) z(t) 图 ３ ０—７８ ｓ 质子 ｘ（ｔ），ｙ（ｔ）和 ｚ（ｔ）曲线 图 ３ 显示 ｘ（ｔ）、ｙ（ｔ）和 ｚ（ｔ）曲线都具有周期 性，ｘ（ｔ）和 ｙ（ｔ）曲线周期对应于粒子漂移运动周期． ｚ（ｔ）曲线周期对应于弹射运动周期，弹射运动周期 远远小于漂移周期． 根据 Ｐａｎｋａｊ Ｋ． Ｓｏｎｉ 等［８］推导的粒子周期理论 公式，粒子弹射运动周期： τｂ ＝ ０．１１７ Ｒ０ Ｒｅ ( ) ｃ ｖ １－０．４６３ ５（ｓｉｎ αｅｑ ） ３ ４ [ ] （４） 粒子漂移运动周期为 τ ｄ ＝ ２πｑＢ０Ｒ３ ｅ ｍｖ２ １ Ｒ０ １－ １ ３ （ｓｉｎ αｅｑ ）０．６２ [ ] （５） 其中 αｅｑ 为带电粒子从赤道面入射时的投掷角，０°≤ αｅｑ≤１８０°，Ｒ０ 为带电粒子入射处距离地球球心的距 离．代入数据得：τｂ ＝ ２．３４０ ５６ ｓ，τ ｄ ＝ ７７．９６３ ６０ ｓ． 计 算图 ３ 中 ｚ（ｔ）曲线与横轴 ｔ 的第三个交点坐标 （２ ３４１ ４５，０），即质子弹射运动周期为 ２．３４１ ４５ ｓ． ｙ（ｔ）曲线与横轴 ｔ 的第 ３ 个交点坐标为（７７．９６３ ５８， ０），即质子漂移运动周期为 ７７．９６３ ５８ ｓ． 可见数值 计算的结果与粒子周期理论值非常吻合． ２．２ 电子在地磁场中的运动 设电子在地磁场中运动的初始条件与 ２．１ 中质 子相同．起始位置（４Ｒｅ，０，０），投掷角 ３０°．以太阳风 中 能 量 为 １０ ＭｅＶ 的 电 子 为 例，电 有 量 ｑ ＝ －１．６０２ １８×１０ － １９Ｃ，ｍ ＝ ９． １０９ ５６ × １０ － ３１ ｋｇ．电子在地 磁场中的运动轨迹如图 ４ 所示． 从北半球方向观察，电子漂移运动沿逆时针方 向，与质子的漂移运动方向刚好相反．电子荷质比非 常大（约是质子的 １ ０００ 倍），其回旋运动半径很小， 因此图 ４ 中几乎看不到回旋运动．再观察 ５ ｓ 内电子 ｚ（ｔ）图像，如图 ５． 可知电子在 ｚ 方向的弹射运动周期非常小，因 此图 ４ 中电子轨迹曲线非常紧密． Copy ropy iopy g h t h g 5 4 3 g h t h g 3 4 55 4 3 g h t h g 3 4 5 ０—２５ ｓ 电子运动轨迹曲线 x/Re z/Re y/Re 2 0 1 1 2 5 5 4 3 2 1 0 1 2 3 4 5 4 3 2 1 0 1 2 3 4 5 ０—１００ ｓ 电子轨迹（顶视） 图 ４ 电子在地磁场中运动轨迹 0 0 1 2 3 4 5 1×107 5×106 5×106 1×107 t/s z/m 图 ５ ０—５ ｓ 电子 ｚ（ｔ）曲线 ３ 带电粒子在地磁场中运动分析 由于电子回旋运动不明显，并且轨迹曲线过于 紧密，因此下面主要以质子为例进行分析讨论．
56大学物理第42卷3.1带电粒子回旋运动南北两极处磁感应强度最大，形成一个中间弱两边强的对称性磁场位形[]，类似于磁镜装置.带电粒带电粒子在磁场中围绕磁感线作回旋运动的拉子从磁场弱处向磁场强处螺旋前进时，会受到一个ymu莫尔半径为(6)p9B与运动方向相反的洛伦兹力F=g（u×B）,随着磁Q=9B场不断增强，粒子沿磁感线方向的速度不断减小，最(7)拉莫尔频率为ym终反射，反射的位置被称为“磁镜点”.设磁镜点处的粒子在地磁场中运动存在一个特殊情况，当粒磁感应强度为B，赤道面处最弱磁场强度B，由子在赤道面上初始速度垂直于磁感线时，粒子将y'm?守恒[2,10],即于磁矩儿直在赤道面上作回旋运动和漂移运动，没有南北两2B极间的弹射运动.以太阳风中能量为30MeV的质子'm (u sin a..)mm?(8)u=为例，其初始位置是（5R。，0,0），投掷角为90°，即入2Bmin2Bmax射速度方向与磁场方向垂直.模拟的运动轨迹以及可得出赤道投掷角各方向的运动如图6和图7所示Baie(9)α=arcsinNB.由此可见，磁镜点的位置只与赤道投角有关，如果磁镜点的位置比较低，在120km左右的电离E层，带电粒子会跟中性的大气或者等离子体发生碰撞，在碰撞过程中会发生光辐射[12]，也就是美丽的0极光现象下面讨论赤道投掷角符合什么条件时粒子会沉降在大气层中，可能产生极光.根据文献【13]，地磁场在球坐标系中可表示为2B.Rcos0B,=x/R,r图60—15s质子运动轨迹(10)BaR'sineB.34×107B.=0其中6为径与≥轴的夹角，设北半球纬度Φ=90°-2×1076，则任意一点地磁感应强度为4coso+sin@I0(11)B=/B,+B,+B, =B,RT6设球坐标系下磁镜点坐标（r，m，0），则由式2×107(11)可得4x1074.cos'o.+sin'o.051015Bm=B,R(12)t/sro2(0)x(t)y(t)图70—100s质子x(t),y(t)和z()曲线设赤道面粒子投射处坐标r。，，则1结合图6和图7，质子一直在赤道面上的运动B..=B.R.(13)V（z（)=0），图6直观形象地展示出质子的回旋运动将式（12)和式（13)代入式（9)，可解得和漂移运动(csc'a.-r, sin'a.3.2南北两极间的弹射运动(14)9=arccos地球偶极磁场模型中，赤道处磁感应强度最小3r大海姚理(C)1994-2023China AcademicJourmalElectrhttpykdxwAbinteegluhnyww.cnki.net
５６ 大 学 物 理 第 ４２ 卷 ３．１ 带电粒子回旋运动 带电粒子在磁场中围绕磁感线作回旋运动的拉 莫尔半径为 ρ ＝ γｍｖ⊥ ｑＢ （６） 拉莫尔频率为 Ω ＝ ｑＢ γｍ （７） 粒子在地磁场中运动存在一个特殊情况，当粒 子在赤道面上初始速度垂直于磁感线时，粒子将一 直在赤道面上作回旋运动和漂移运动，没有南北两 极间的弹射运动．以太阳风中能量为 ３０ ＭｅＶ 的质子 为例，其初始位置是（５Ｒｅ，０，０），投掷角为 ９０°，即入 射速度方向与磁场方向垂直．模拟的运动轨迹以及 各方向的运动如图 ６ 和图 ７ 所示． Copy ropy iopy 5 2 1 0 1 2 6 5 6 4 3 2 1 0 1 2 3 4 5 6 4 3 2 1 0 1 2 3 4 5 6 图 ６ ０—１５ ｓ 质子运动轨迹 0 0 5 10 15 4×107 2×107 2×107 4×107 t/s L/m x(t) y(t) z(t) 图 ７ ０—１００ ｓ 质子 ｘ（ｔ），ｙ（ｔ）和 ｚ（ｔ）曲线 结合图 ６ 和图 ７，质子一直在赤道面上的运动 （ｚ（ｔ）＝ ０），图 ６ 直观形象地展示出质子的回旋运动 和漂移运动． ３．２ 南北两极间的弹射运动 地球偶极磁场模型中，赤道处磁感应强度最小， 南北两极处磁感应强度最大，形成一个中间弱两边 强的对称性磁场位形［１１］ ，类似于磁镜装置．带电粒 子从磁场弱处向磁场强处螺旋前进时，会受到一个 与运动方向相反的洛伦兹力 Ｆ ＝ ｑ（ｖ⊥ ×Ｂ），随着磁 场不断增强，粒子沿磁感线方向的速度不断减小，最 终反射，反射的位置被称为“磁镜点”．设磁镜点处的 磁感应强度为 Ｂ ｍａｘ，赤道面处最弱磁场强度 Ｂ ｍｉｎ ，由 于磁矩 μ ＝ γ２ ｍｖ２ ⊥ ２Ｂ 守恒［２，１０］ ，即 μ ＝ γ２ ｍ （ｖ ｓｉｎ αｅｑ ）２ ２Ｂ ｍｉｎ ＝ γ２ ｍｖ２ ２Ｂ ｍａｘ （８） 可得出赤道投掷角 αｅｑ ＝ ａｒｃｓｉｎ Ｂ ｍｉｎ 槡Ｂ ｍａｘ （９） 由此可见，磁镜点的位置只与赤道投掷角有关． 如果磁镜点的位置比较低，在 １２０ ｋｍ 左右的电离 Ｅ 层，带电粒子会跟中性的大气或者等离子体发生碰 撞，在碰撞过程中会发生光辐射［１２］ ，也就是美丽的 极光现象． 下面讨论赤道投掷角符合什么条件时粒子会沉 降在大气层中，可能产生极光．根据文献［１３］ ，地磁场 在球坐标系中可表示为 Ｂｒ ＝ －２Ｂ０Ｒ３ ｅ ｃｏｓ θ ｒ ３ ， Ｂθ ＝ －Ｂ０Ｒ３ ｅ ｓｉｎ θ ｒ ３ ， Ｂφ ＝ ０ （１０） 其中 θ 为矢径与 ｚ 轴的夹角，设北半球纬度 ＝ ９０°－ θ，则任意一点地磁感应强度为 Ｂ ＝ Ｂ２ ｒ ＋Ｂθ ２ ＋Ｂ２ 槡 φ ＝ Ｂ０Ｒ３ ｅ ４ ｃｏｓ２ θ＋ｓｉｎ２ θ ｒ 槡 ６ （１１） 设球坐标系下磁镜点坐标（ｒｍ ，θｍ ，０），则由式 （１１）可得 Ｂ ｍａｘ ＝ Ｂ０Ｒ３ ｅ ４ ｃｏｓ２ θｍ ＋ｓｉｎ２ θｍ ｒ ６ 槡 ｍ （１２） 设赤道面粒子投射处坐标 ｒ０ ，π ２ ，０ r i ，则 Ｂ ｍｉｎ ＝ Ｂ０Ｒ３ ｅ １ ｒ ６ 槡０ （１３） 将式（１２）和式（１３）代入式（９），可解得 θｍ ＝ ａｒｃｃｏｓ ｃｓｃ２ αｅｑ ｒ ６ ｍ －ｒ ６ ０ ｓｉｎ４ 槡 αｅｑ 槡３ ｒ０ ３ （１４）
第7期57董顺成，等：带电粒子在地磁场中的运动及Mathematica数值模拟设中.为磁镜点对应的北半球的纬度，中，=90°-9.当磁镜点位于距离地面120km的大气层时，带电粒子会进人大气层，r.=R+120km.以投挪0点坐标4R。，0为例，绘制质子的中~α曲线，2如图8所示-3--R.X/R.80(a)投挪角=30°6040N0-30VR图8投挪点坐标为4R。0时中~α曲线X/R.(b)投角am=5.30由图8可知，当投掷角满足5.21°<α<7.38°存在中，此时磁镜点正好位于距离地球表面120km的大气层面上，粒子将沉降到距离地球表面120km的大气层中；投掷角小于5.21°时，不存在对dRe应的磁镜点纬度坐标中.值，此时磁镜点位于距离地球表面120km的大气层内，带电粒子尚未到达磁镜点就已沉降.因此，只要投挪角α.<7.38°，有可能WRx/Re产生极光(c)投掷角α.=0°随着α的减小，中.迅速增大，投掷角α=图9从赤道面投掷时，一个弹射周期内质子运动轨迹5.50°时对应的中.已经达到60°，意味着投掷角α<5.50°的带电粒子，将沉降到60°以上高纬度地区大3.3漂移运动气层中，因此在高纬度地区更容易出现神奇的极光带电粒子在地磁场中漂移运动包括梯度漂移和现象.以能量10MeV的质子为例，投掷角分别为曲率漂移[14].由于磁场大小的不均匀性，产生梯度30°5.2°和0°，各自一个弹射周期内轨迹如图9漂移，对应的梯度漂移速度公式[7]为图9（a）所示的磁镜点距离地球很远，粒子会1 mRxBm直束缚在地球附近，形成范艾伦辐射带：图9（b）中(15)(B×VB)B"2B磁镜点已经非常接近地球，并且位于地球的高纬度由于场方向的不均匀性，粒子沿磁场方向弯曲地区，图9（c）中投掷角为0°，带电粒子几乎直接射前行会受到离心力作用，产生曲率漂移，对应的曲率入地球，实际上这种情况并不会发生，当带电粒子非漂移速度公式[15]为常接近地球时，尚未到达磁镜点就已经和大气层发m, Rx×B_mi(BxV B)生碰撞而沉降，(16)UReqB2RqB综上所述，带电粒子在一个弹射周期内是否会式中R，为磁感线曲率半径，方向由曲率中心指向沉降入大气层中甚至能否产生极光现象，与带电粒外.由式（15）和（16），总的梯度一曲率漂移速度子投掷角有关，投掷角大时，粒子会被束缚在地磁场，投挪角比较小时，带电粒子会沉降，存在产生极Dm((BxB)(17)U=UVE+URED.+qB3光的可能，并且粒子更容易在高纬度地区沉降入大2根据式（17），漂移方向与粒子的正负有关，前面气层.大学如理(C)1994-2023ChinaAcademicJourmal Electrhttpy/dxwAbt/eecluh@nwww.cnki.net
第 ７ 期 董顺成，等：带电粒子在地磁场中的运动及 Ｍａｔｈｅｍａｔｉｃａ 数值模拟 ５７ 设 ｍ 为磁镜点对应的北半球的纬度， ｍ ＝ ９０°－θｍ ． 当磁镜点位于距离地面 １２０ ｋｍ 的大气层 时，带电粒子会进入大气层，ｒｍ ＝ Ｒｅ ＋１２０ ｋｍ．以投掷 点坐标 ４Ｒｅ，π ２ ，０ ( ) 为例，绘制质子的 ｍ ～ αｅｑ 曲线， 如图 ８ 所示． Co po yo ro y i p g C h t( c) 图 ８ 投掷点坐标为 ４Ｒｅ，π ２ ，０ ( ) 时 ｍ ～ αｅｑ曲线 由图 ８ 可知，当投掷角满足 ５ ．２１° ＜αｅｑ ＜ ７ ．３８°， 存在 ｍ ，此 时 磁 镜 点 正 好 位 于 距 离 地 球 表 面 １２０ ｋｍ 的大气层面上，粒子将沉降到距离地球表面 １２０ ｋｍ 的大气层中；投掷角小于 ５ ．２１°时，不存在对 应的磁镜点纬度坐标 ｍ 值，此时磁镜点位于距离 地球表面 １２０ ｋｍ 的大气层内，带电粒子尚未到达磁 镜点就已沉降．因此，只要投掷角 αｅｑ ＜ ７ ．３８°，有可能 产生极光． 随着 αｅｑ 的减小， ｍ 迅速增大，投掷角 αｅｑ ＝ ５．５０°时对应的 ｍ 已经达到 ６０°，意味着投掷角αｅｑ ＜ ５．５０°的带电粒子，将沉降到 ６０°以上高纬度地区大 气层中，因此在高纬度地区更容易出现神奇的极光 现象．以能量 １０ ＭｅＶ 的质子为例，投掷角分别为 ３０°、５．２°和 ０°，各自一个弹射周期内轨迹如图 ９． 图 ９（ａ）所示的磁镜点距离地球很远，粒子会一 直束缚在地球附近，形成范艾伦辐射带；图 ９（ｂ）中 磁镜点已经非常接近地球，并且位于地球的高纬度 地区，图 ９（ｃ）中投掷角为 ０°，带电粒子几乎直接射 入地球，实际上这种情况并不会发生，当带电粒子非 常接近地球时，尚未到达磁镜点就已经和大气层发 生碰撞而沉降． 综上所述，带电粒子在一个弹射周期内是否会 沉降入大气层中甚至能否产生极光现象，与带电粒 子投掷角有关，投掷角大时，粒子会被束缚在地磁 场，投掷角比较小时，带电粒子会沉降，存在产生极 光的可能，并且粒子更容易在高纬度地区沉降入大 气层． 2 0 - 0 2 5 4 3 2 0 - 0 2 3 - 4 5 x/Re z/Re y/Re 投掷角 αｅｑ ＝ ３０° 2 1 0 1 2 0 x/Re z/Re y/Re 5 4 3 2 1 0 1 2 3 4 5 投掷角 αｅｑ ＝ ５ ．３° 2 1 0 1 2 0 x/Re z/Re y/Re 5 4 3 2 1 0 1 2 3 4 5 投掷角 αｅｑ ＝ ０° 图 ９ 从赤道面投掷时，一个弹射周期内质子运动轨迹 ３．３ 漂移运动 带电粒子在地磁场中漂移运动包括梯度漂移和 曲率漂移［１４］ ． 由于磁场大小的不均匀性，产生梯度 漂移，对应的梯度漂移速度公式［７］ 为 ｖ Ｂ ＝ １ ２ ｍｖ２ ⊥ ｑＢ２ Ｒｃ ×Ｂ Ｒ２ ｃ ＝ ｍｖ２ ⊥ ２ｑＢ３ （Ｂ× Ｂ） （１５） 由于场方向的不均匀性，粒子沿磁场方向弯曲 前行会受到离心力作用，产生曲率漂移，对应的曲率 漂移速度公式［１５］ 为 ｖＲｃ ＝ ｍｖ２ ‖ ｑＢ２ Ｒｃ ×Ｂ Ｒ２ ｃ ＝ ｍｖ２ ‖ ｑＢ３ （Ｂ× Ｂ） （１６） 式中 Ｒｃ 为磁感线曲率半径，方向由曲率中心指向 外．由式（１５）和（１６），总的梯度—曲率漂移速度 ｖ ＝ ｖ Ｂ ＋ｖＲｃ ＝ ｍ ｑＢ３ ｖ ２ ‖ ＋ ｖ ２ ⊥ ( ) ２ （Ｂ× Ｂ） （１７） 根据式（１７），漂移方向与粒子的正负有关，前面
58大学物理第42卷的分析也印证这一点，从北半球来看，正电荷沿顺时可以发现，其他条件相同时，粒子投掷角越大，其对应的漂移速度也越大针方向漂移，负电荷沿逆时针方向飘移.漂移速度的大小与多重因素有关，比如粒子沿磁感线的速度（与4引导中心投掷角有关），磁场的大小（与人射位置有关）等在单粒子方法研究中，粒子的运动往往非常复3.3.1入射位置对带电粒子漂移速度大小的影响杂，研究者引人了引导中心近似概念，引导中心是粒子做回旋运动的几何中心[]，如图12所示.设质子分别从赤道面（3R0.0）、（4R0.0）和引导中心（5R.0.0）位置上入射，动能为10MeV，投挪角30°粒子保持不变，运行Mathematica软件30s，模拟结果如图10.A6F54H323R1o。4R.图12带电粒子引导中心R示意图-15R.-2R为粒子引导中心位置矢量，为粒子的瞬时位置，-3p为粒子回旋半径，它们之间的关系为-4r=R+p(18)06-5-4-3-2-1012345Z/R.徐荣栏等叫指出，当粒子围绕引导中心一周x/R时，如果磁场的空间与时间T满足图10不同投射位置对应的质子运动轨迹（顶视图）rl(p B) /<IB/.(19)dB</BI1由图10可以看出，在其他条件相同的情况下，dt粒子入射位置距离球心越远，其漂移速度越大引导中心的概念仍然成立.此条件被称为小扰动条件，当粒子运动满足这个条件时，磁场可以利用3.3.2投掷角对带电粒子漂移速度大小的影响泰勒级数展开近似为设粒子投掷角分别为20°40°和60°以10MeVB(r)-B(R)+(p. V)B(R)(20)动能从（4R，0,0）投掷，运行模拟软件30s，得到的Walt[13]等人得出引导中心方程：轨迹曲线如图11.mR.b=-qp(B)·b2d,-1或(VB),(21)3Hdt2 B2B其中b=为磁场方向单位失量1200B400Oztirk[2]在此基础上进一步计算得出相对论条600件下引导中心方程组：-2dR_ymuxVB+u.bdt2qB2(22)de-b.B-10z/R.dtmx/R通过引人引导中心近似，研究带电粒子在地磁图11不同投掷角对应的质子运动轨迹（顶视图）大海城理(C)1994-2023China AcademicJournalElectrhttpy/kdxwAbit/ee@luh@nyww.cnki.net
５８ 大 学 物 理 第 ４２ 卷 的分析也印证这一点，从北半球来看，正电荷沿顺时 针方向漂移，负电荷沿逆时针方向飘移．漂移速度的 大小与多重因素有关，比如粒子沿磁感线的速度（与 投掷角有关），磁场的大小（与入射位置有关）等． ３．３．１ 入射位置对带电粒子漂移速度大小的影响 设质子分别从赤道面（３Ｒｅ，０，０）、（４Ｒｅ，０，０）和 （５Ｒｅ，０，０）位置上入射，动能为 １０ ＭｅＶ，投掷角 ３０° 保持不变，运行 Ｍａｔｈｅｍａｔｉｃａ 软件 ３０ ｓ，模拟结果如 图 １０． Copy ropy iopy 0 6 5 5 6 4 3 2 1 0 1 2 3 4 5 6 4 3 2 1 0 1 2 3 4 5 6 3py 4py 5py 图 １０ 不同投射位置对应的质子运动轨迹（顶视图） 由图 １０ 可以看出，在其他条件相同的情况下， 粒子入射位置距离球心越远，其漂移速度越大． ３．３．２ 投掷角对带电粒子漂移速度大小的影响 设粒子投掷角分别为 ２０°、４０°和 ６０°，以 １０ ＭｅＶ 动能从（４Ｒｅ，０，０）投掷，运行模拟软件 ３０ ｓ，得到的 轨迹曲线如图 １１． x/Re z/Re 5 5 4 0 3 2 1 0 1 2 3 4 5 4 3 2 1 0 1 2 3 4 5 y/Re 20° 40° 60° 图 １１ 不同投掷角对应的质子运动轨迹（顶视图） 可以发现，其他条件相同时，粒子投掷角越大， 其对应的漂移速度也越大． ４ 引导中心 在单粒子方法研究中，粒子的运动往往非常复 杂，研究者引入了引导中心近似概念，引导中心是粒 子做回旋运动的几何中心［２］ ，如图 １２ 所示． R r B O ρ 图 １２ 带电粒子引导中心 Ｒ 示意图 Ｒ 为粒子引导中心位置矢量，ｒ 为粒子的瞬时位置， ρ 为粒子回旋半径，它们之间的关系为 ｒ ＝ Ｒ＋ρ （１８） 徐荣栏等［１］指出，当粒子围绕引导中心一周 时，如果磁场的空间与时间 Ｔ 满足 r i ρ· Ｂ ＜＜ Ｂ ， Ｔ ｄＢ ｄｔ p ＜＜ Ｂ （１９） 引导中心的概念仍然成立．此条件被称为小扰 动条件，当粒子运动满足这个条件时，磁场可以利用 泰勒级数展开近似为 Ｂ（ｒ）≈Ｂ（Ｒ）＋（ρ· ）Ｂ（Ｒ） （２０） Ｗａｌｔ［１３］ 等人得出引导中心方程： ｍＲ ¨ ·ｂ ＝ －ｑρ ２ Ω ２ （ Ｂ）·ｂ 或 ｄｖ‖ ｄｔ ＝ － １ ２ ｖ ２ ⊥ Ｂ （ Ｂ）‖ （２１） 其中 ｂ ＝ Ｂ Ｂ 为磁场方向单位矢量． Ｏ ¨ ｚｔ ｕ¨ ｒｋ［２］ 在此基础上进一步计算得出相对论条 件下引导中心方程组： ｄＲ ｄｔ ＝ γｍｖ２ ２ｑＢ２ １＋ ｖ ２ ‖ ｖ r i ２ ｂ× Ｂ＋ｖ‖ ｂ， ｄｖ‖ ｄｔ ＝ － μ γ２ ｍ ｂ· Ｂ （２２） 通过引入引导中心近似，研究带电粒子在地磁
59第7期董顺成，等：带电粒子在地磁场中的运动及Mathematica数值模拟场中运动时，可以不用考虑小尺度回旋运动，从而考MROT虑更大时空尺度的弹射运动和漂移运动以太阳风中能量为10MeV的质子为例，粒子起始位置为（4R，0,0），投掷角30°，60s内引导中心引导中心轨连轨迹曲线如图13所示质子运动轨迹R引导中心轨迹质子运动轨迹图151MeV质子0~500引导中心和实际运动轨迹JR.1×1075x106图1310MeV质子0~60s引导中心及实际运动轨迹I5x106根据图13.在粒子运动的前15s.引导中心能很好的表征质子实际运动，之后引导中心曲线与粒X101002003004000500子运动曲线逐渐偏离，趋势是引导中心曲线越来越1/s往z轴正方向偏移.图14则显示出引导中心在z轴引导中心Z0)质子运动2(0)正方向的偏离情况图16引导中心Z(）)和质子运动z(t)）曲线1.5x10图16为引导中心Z(）和质子z()运动曲线，虽然1.0×1在300s之后Z()略高于z(），但与图15相比误差已5.0×10I引导中心20经很小，几乎可以忽略不计进一步对能量更小的质子运动(0)0.5MeV粒子进行数值模拟，二者轨迹几乎完全重合综上所述，文献2中推导的引导中心方程适208合能量较低的带电粒子，当粒子能量过高以至于带电粒子围绕磁感线的回旋半径过大时，引导中心方图1410MeV质子引导中心Z(）和质子运动z(）曲线程不适用，下一步我们将尝试推导二阶近似下的引导中心方程，从而实现适用性更好的引导中心近似两条曲线振荡频率完全同步，可以用来表征漂移运动.但引导中心的Z（）越来越往上偏离，此时5结论用引导中心运动来表征质子在南北两极间的弹射运本文基于地磁场磁偶极子模型和单粒子模型，动不再合适.原因是当质子动能过大时，根据式（6），考虑相对论效应，研究了近地区域（6R以内）中高其回旋半径也过大，螺旋特征较为明显，而式（20）能带电粒子在近地球区的运动问题，利用泰勒级数展开中仅使用了回旋半径中的一阶项，因Mathematica软件使用精确度更高的六阶龙格一库此引导中心与实际粒子运动有很大差异，塔算法对粒子的动力学方程组进行数值求解并模下面选择动能较小的粒子，以太阳风中能量为拟，具体结论如下：1MeV的质子为例，其他条件保持不变，观察二者轨1）从地球北极方向观察，被地球捕获的质子沿迹.运行500s结果如图15所示，顺时针方向漂移，电子则沿逆时针方向漂移.其他条由图15，当质子动能比较小时，质子螺旋运动件相同时，带电粒子投掷点距离地球越远，其漂移速特征并不明显，引导中心曲线与实际运动轨道高度度越大.投掷角越大，其漂移速度也越大；吻合，此时引导中心较适合表征质子的实际运动。大海好理(C)1994-2023ChinaAcademicJoumal Electrdhttpy//dxwAbt/eeduh@nvww.cnkinct
第 ７ 期 董顺成，等：带电粒子在地磁场中的运动及 Ｍａｔｈｅｍａｔｉｃａ 数值模拟 ５９ 场中运动时，可以不用考虑小尺度回旋运动，从而考 虑更大时空尺度的弹射运动和漂移运动． 以太阳风中能量为 １０ ＭｅＶ 的质子为例，粒子 起始位置为（４Ｒｅ，０，０），投掷角３０°，６０ ｓ 内引导中心 轨迹曲线如图 １３ 所示． Copy ropy i g h g i 5 4 3 i g h g i 3 4 5 yopy 5 4 3 i g h g i 3 4 5 图 １３ １０ ＭｅＶ 质子 ０ ～ ６０ ｓ 引导中心及实际运动轨迹 根据图 １３，在粒子运动的前 １５ ｓ，引导中心能 很好的表征质子实际运动，之后引导中心曲线与粒 子运动曲线逐渐偏离，趋势是引导中心曲线越来越 往 ｚ 轴正方向偏移．图 １４ 则显示出引导中心在 ｚ 轴 正方向的偏离情况． 0 0 20 40 60 5.0×106 5.0×106 1.0×107 1.5×107 1.0×107 t/s L/m Z(t) z(t) 图 １４ １０ ＭｅＶ 质子引导中心 Ｚ（ｔ）和质子运动 ｚ（ｔ）曲线 两条曲线振荡频率完全同步，可以用来表征漂 移运动．但引导中心的 Ｚ（ｔ）越来越往上偏离，此时 用引导中心运动来表征质子在南北两极间的弹射运 动不再合适．原因是当质子动能过大时，根据式（６）， 其回旋半径也过大，螺旋特征较为明显，而式（２０） 泰勒级数展开中仅使用了回旋半径中的一阶项，因 此引导中心与实际粒子运动有很大差异． 下面选择动能较小的粒子，以太阳风中能量为 １ ＭｅＶ 的质子为例，其他条件保持不变，观察二者轨 迹．运行 ５００ ｓ 结果如图 １５ 所示． 由图 １５，当质子动能比较小时，质子螺旋运动 特征并不明显，引导中心曲线与实际运动轨道高度 吻合，此时引导中心较适合表征质子的实际运动． x/Re z/Re 2 1 0 1 2 5 4 3 2 1 0 1 2 3 4 5 y/Re 5 4 3 2 1 0 1 2 3 4 5 图 １５ １ ＭｅＶ 质子 ０ ～ ５００ ｓ 引导中心和实际运动轨迹 0 0 100 200 300 400 500 5×106 5×106 1×107 1×107 t/s L/m Z(t) z(t) 图 １６ 引导中心 Ｚ（ｔ）和质子运动 ｚ（ｔ）曲线 图 １６ 为引导中心 Ｚ（ｔ）和质子ｚ（ｔ）运动曲线，虽然 在 ３００ ｓ 之后 Ｚ（ｔ）略高于 ｚ（ｔ），但与图 １５ 相比误差已 经很小，几乎可以忽略不计．进一步对能量更小的 ０．５ ＭｅＶ 粒子进行数值模拟，二者轨迹几乎完全重合． 综上所述，文献［２］中推导的引导中心方程适 合能量较低的带电粒子，当粒子能量过高以至于带 电粒子围绕磁感线的回旋半径过大时，引导中心方 程不适用，下一步我们将尝试推导二阶近似下的引 导中心方程，从而实现适用性更好的引导中心近似． ５ 结论 本文基于地磁场磁偶极子模型和单粒子模型， 考虑相对论效应，研究了近地区域（６Ｒｅ 以内）中高 能 带 电 粒 子 在 近 地 球 区 的 运 动 问 题，利 用 Ｍａｔｈｅｍａｔｉｃａ 软件使用精确度更高的六阶龙格—库 塔算法对粒子的动力学方程组进行数值求解并模 拟，具体结论如下： １）从地球北极方向观察，被地球捕获的质子沿 顺时针方向漂移，电子则沿逆时针方向漂移．其他条 件相同时，带电粒子投掷点距离地球越远，其漂移速 度越大．投掷角越大，其漂移速度也越大；
60第42卷大学物理运动［J].南昌大学学报（理科版），2016，40（01）：2）粒子各个分运动的运动周期数值模拟结果35-38.与文献中周期理论值非常吻合；[7]钟海坚，陈宗华，赵炳炎.基于MATLAB的地磁场中带3）从（4R。，0，0）人射的粒子投挪角小于7.38°电粒子运动模拟分析[J].大学物理实验，2021，34时，将会与地球表面大气层碰撞而沉降，可能会产生(01) :83-86极光现象.大于7.38°时，带电粒子将会被地磁场捕[8]Pankaj K. Soni, Bharati Kakad, Amar Kakad. Simulation获，形成对地球具有保护作用的辐射带；study of motion of charged particles trapped in Earth's4）文献[2]中推导的引导中心方程适合能量magnetosphere[J].Advances in Space Research, 202],较低的带电粒子.对于能量较高的粒子，可能需要考67(2):749-761.虑二阶近似。[9】刘振兴.太空物理学【M].哈尔滨：哈尔滨工业大学，2005.参考文献：[10] 马腾才，胡希伟，陈银华.等离子体物理原理［M].合肥：中国科学技术出版社，2012.[]]徐荣栏，李磊.磁层粒子动力学M1.北京：科学出版[1]刘爱红，余守宪.轴对称缓变磁场中的等离子体与磁社，2005.约束原理[J].物理与工程，2002（05）：14-17+23.M.Kaan Ozt irk.Trajectories of charged particles trapped[2] [12]周苏.极区粒子沉降能量变化及其形成机制研究in earth's magnetic field[J]. American journal of[D].中国科学技术大学，2018.physics,2012,80(05):420-428.[13] Martin Walt. Introduction to geomagnetically trapped ra-[3] Carl stomer.The polar aurora[M]. London: Oxford Uni-diation ［ MJ.New York: Cambridge UniversityversityPress,1955:213.Press; 1994.[4] 薛林.带电粒子在地磁场中的运动[J].大学物理，【14]李定，陈银华，马锦秀，等.等离子体物理学【M].北2006,(08):56-58京：高等教育出版社，2006.[5] 黄朝艳，韩建伟.带电粒子在近地球区运动的计算方法[15] 陈耀.等离子体物理学基础【M].北京：科学出版J].空间科学学报，2007，（05）：367-373.社,2019.[6] 代国红，黄伟军，肖伟，等.带电粒子在地磁场中的漂移Motion of charged particles in earth's magneticfield and Mathematica numerical simulationDONG Shun-cheng, GUO Fang-xia(School of Physics and Information Technology,Shaanxi Normal University, Xi'an, Shaanxi 710119, China)Abstract:Based on the single particle orbit model and the geomagnetic field dipole model, and considering therelativistic effect, the sixth-order Runge-Kutta algorithm in the Mathematica software is used to numerically calcu-late and simulate the trajectories of charged particles moving in the magnetic field in the near-Earth region. Thegeneration of auroras is explained, and the approximation of the guiding center of the charged particles moving inthe geomagnetic field is discussed. The results show that: (i) Observed from the earth's north pole, the protonscaptured by the earth's magnetic field drift in the clockwise direction, and the electrons in the anti-clockwise direc-tion;(ii)The numerical simulation results of the motion period of each particle motion are in good agreement withthe theoretical values in the literature;(ili) When the throwing angle of the incident particle from (4 R,O,0) isless than 7.38°, the charged particle will collide with the atmosphere on the earth's surface and sink, and there isthe possibility of producing auroras. When it is greater than 7.38°, the particles will be bound in the geomagneticfield, forming a radiation belt; (iv) When other conditions are the same, the farther the charged particle throwingpoint is from the earth, the greater the drift velocity; the greater the throwing angle, the greater the Its drift speed isalso greater; (v) For particles with low energy,the guided center trajectory can well represent the actual trajecto-ry of particles in the first-order approximation.Key words:earth's magnetic field;the particle motion; numerical simulation大海好理(C)1994-2023ChinaAcademicJourmal Electrhttpy//dxwAbit/ee@uh@nvww.cnki.net
６０ 大 学 物 理 第 ４２ 卷 ２）粒子各个分运动的运动周期数值模拟结果 与文献中周期理论值非常吻合； ３）从（４ Ｒｅ，０，０）入射的粒子投掷角小于７．３８° 时，将会与地球表面大气层碰撞而沉降，可能会产生 极光现象． 大于７．３８°时，带电粒子将会被地磁场捕 获，形成对地球具有保护作用的辐射带； ４）文献［２］中推导的引导中心方程适合能量 较低的带电粒子．对于能量较高的粒子，可能需要考 虑二阶近似． 参考文献： ［１］ 徐荣栏，李磊．磁层粒子动力学［Ｍ］．北京：科学出版 社，２００５． ［２］ Ｍ．Ｋａａｎ Ｏ¨ ｚｔ ｕ¨ ｒｋ． Ｔｒａｊｅｃｔｏｒｉｅｓ ｏｆ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅｓ ｔｒａｐｐｅｄ ｉｎ ｅａｒｔｈ ’ｓ ｍａｇｎｅｔｉｃ ｆｉｅｌｄ ［Ｊ］． Ａｍｅｒｉｃａｎ ｊｏｕｒｎａｌ ｏｆ ｐｈｙｓｉｃｓ，２０１２，８０（０５）：４２０４２８． ［３］ Ｃａｒｌ ｓｔｏｍｅｒ ¨ ． Ｔｈｅ ｐｏｌａｒ ａｕｒｏｒａ［Ｍ］． Ｌｏｎｄｏｎ：Ｏｘｆｏｒｄ Ｕｎｉ ｖｅｒｓｉｔｙ Ｐｒｅｓｓ，１９５５：２１３． ［４］ 薛林．带电粒子在地磁场中的运动 ［Ｊ］．大学物理， ２００６，（０８）：５６５８． ［５］ 黄朝艳，韩建伟．带电粒子在近地球区运动的计算方法 ［Ｊ］．空间科学学报，２００７，（０５）：３６７３７３． ［６］ 代国红，黄伟军，肖伟，等．带电粒子在地磁场中的漂移 运动 ［Ｊ］．南昌大学学报 （理科版），２０１６，４０ （０１）： ３５３８． ［７］ 钟海坚，陈宗华，赵炳炎．基于 ＭＡＴＬＡＢ 的地磁场中带 电粒子运动模拟分析 ［Ｊ］．大学物理实验，２０２１，３４ （０１）：８３８６． ［８］ Ｐａｎｋａｊ Ｋ． Ｓｏｎｉ，Ｂｈａｒａｔｉ Ｋａｋａｄ，Ａｍａｒ Ｋａｋａｄ． Ｓｉｍｕｌａｔｉｏｎ ｓｔｕｄｙ ｏｆ ｍｏｔｉｏｎ ｏｆ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅｓ ｔｒａｐｐｅｄ ｉｎ Ｅａｒｔｈ’ｓ ｍａｇｎｅｔｏｓｐｈｅｒｅ［Ｊ］． Ａｄｖａｎｃｅｓ ｉｎ Ｓｐａｃｅ Ｒｅｓｅａｒｃｈ，２０２１， ６７（２）：７４９７６１． ［９］ 刘振兴．太空物理学［Ｍ］．哈尔滨：哈尔滨工业大学， ２００５． ［１０］ 马腾才，胡希伟，陈银华． 等离子体物理原理［Ｍ］． 合 肥：中国科学技术出版社，２０１２． ［１１］ 刘爱红，佘守宪．轴对称缓变磁场中的等离子体与磁 约束原理［Ｊ］．物理与工程，２００２（０５）：１４１７＋２３． ［１２］ 周苏．极区粒子沉降能量变化及其形成机制研究 ［Ｄ］．中国科学技术大学，２０１８． ［１３］ Ｍａｒｔｉｎ Ｗａｌｔ． Ｉｎｔｒｏｄｕｃｔｉｏｎ ｔｏ ｇｅｏｍａｇｎｅｔｉｃａｌｌｙ ｔｒａｐｐｅｄ ｒａ ｄｉａｔｉｏｎ ［Ｍ ］． Ｎｅｗ Ｙｏｒｋ： Ｃａｍｂｒｉｄｇｅ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，１９９４． ［１４］ 李定，陈银华，马锦秀，等．等离子体物理学［Ｍ］．北 京：高等教育出版社，２００６． ［１５］ 陈耀．等离子体物理学基础 ［Ｍ ］．北京：科学出版 社，２０１９． Ｍｏｔｉｏｎ ｏｆ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅｓ ｉｎ ｅａｒｔｈ’ｓ ｍａｇｎｅｔｉｃ ｆｉｅｌｄ ａｎｄ Ｍａｔｈｅｍａｔｉｃａ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ＤＯＮＧ Ｓｈｕｎｃｈｅｎｇ，ＧＵＯ Ｆａｎｇｘｉａ （Ｓｃｈｏｏｌ ｏｆ Ｐｈｙｓｉｃｓ ａｎｄ Ｉｎｆｏｒｍａｔｉｏｎ Ｔｅｃｈｎｏｌｏｇｙ，Ｓｈａａｎｘｉ Ｎｏｒｍａｌ Ｕｎｉｖｅｒｓｉｔｙ，Ｘｉ’ａｎ，Ｓｈａａｎｘｉ ７１０１１９，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｂａｓｅｄ ｏｎ ｔｈｅ ｓｉｎｇｌｅ ｐａｒｔｉｃｌｅ ｏｒｂｉｔ ｍｏｄｅｌ ａｎｄ ｔｈｅ ｇｅｏｍａｇｎｅｔｉｃ ｆｉｅｌｄ ｄｉｐｏｌｅ ｍｏｄｅｌ，ａｎｄ ｃｏｎｓｉｄｅｒｉｎｇ ｔｈｅ ｒｅｌａｔｉｖｉｓｔｉｃ ｅｆｆｅｃｔ，ｔｈｅ ｓｉｘｔｈ－ｏｒｄｅｒ ＲｕｎｇｅＫｕｔｔａ ａｌｇｏｒｉｔｈｍ ｉｎ ｔｈｅ Ｍａｔｈｅｍａｔｉｃａ ｓｏｆｔｗａｒｅ ｉｓ ｕｓｅｄ ｔｏ ｎｕｍｅｒｉｃａｌｌｙ ｃａｌｃｕ ｌａｔｅ ａｎｄ ｓｉｍｕｌａｔｅ ｔｈｅ ｔｒａｊｅｃｔｏｒｉｅｓ ｏｆ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅｓ ｍｏｖｉｎｇ ｉｎ ｔｈｅ ｍａｇｎｅｔｉｃ ｆｉｅｌｄ ｉｎ ｔｈｅ ｎｅａｒＥａｒｔｈ ｒｅｇｉｏｎ． Ｔｈｅ ｇｅｎｅｒａｔｉｏｎ ｏｆ ａｕｒｏｒａｓ ｉｓ ｅｘｐｌａｉｎｅｄ，ａｎｄ ｔｈｅ ａｐｐｒｏｘｉｍａｔｉｏｎ ｏｆ ｔｈｅ ｇｕｉｄｉｎｇ ｃｅｎｔｅｒ ｏｆ ｔｈｅ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅｓ ｍｏｖｉｎｇ ｉｎ ｔｈｅ ｇｅｏｍａｇｎｅｔｉｃ ｆｉｅｌｄ ｉｓ ｄｉｓｃｕｓｓｅｄ． Ｔｈｅ ｒｅｓｕｌｔｓ ｓｈｏｗ ｔｈａｔ：（ⅰ）Ｏｂｓｅｒｖｅｄ ｆｒｏｍ ｔｈｅ ｅａｒｔｈ’ｓ ｎｏｒｔｈ ｐｏｌｅ，ｔｈｅ ｐｒｏｔｏｎｓ ｃａｐｔｕｒｅｄ ｂｙ ｔｈｅ ｅａｒｔｈ’ｓ ｍａｇｎｅｔｉｃ ｆｉｅｌｄ ｄｒｉｆｔ ｉｎ ｔｈｅ ｃｌｏｃｋｗｉｓｅ ｄｉｒｅｃｔｉｏｎ，ａｎｄ ｔｈｅ ｅｌｅｃｔｒｏｎｓ ｉｎ ｔｈｅ ａｎｔｉｃｌｏｃｋｗｉｓｅ ｄｉｒｅｃ ｔｉｏｎ；（ⅱ）Ｔｈｅ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｒｅｓｕｌｔｓ ｏｆ ｔｈｅ ｍｏｔｉｏｎ ｐｅｒｉｏｄ ｏｆ ｅａｃｈ ｐａｒｔｉｃｌｅ ｍｏｔｉｏｎ ａｒｅ ｉｎ ｇｏｏｄ ａｇｒｅｅｍｅｎｔ ｗｉｔｈ ｔｈｅ ｔｈｅｏｒｅｔｉｃａｌ ｖａｌｕｅｓ ｉｎ ｔｈｅ ｌｉｔｅｒａｔｕｒｅ；（ⅲ）Ｗｈｅｎ ｔｈｅ ｔｈｒｏｗｉｎｇ ａｎｇｌｅ ｏｆ ｔｈｅ ｉｎｃｉｄｅｎｔ ｐａｒｔｉｃｌｅ ｆｒｏｍ （４ Ｒｅ，０，０）ｉｓ ｌｅｓｓ ｔｈａｎ ７．３８°，ｔｈｅ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅ ｗｉｌｌ ｃｏｌｌｉｄｅ ｗｉｔｈ ｔｈｅ ａｔｍｏｓｐｈｅｒｅ ｏｎ ｔｈｅ ｅａｒｔｈ’ｓ ｓｕｒｆａｃｅ ａｎｄ ｓｉｎｋ，ａｎｄ ｔｈｅｒｅ ｉｓ ｔｈｅ ｐｏｓｓｉｂｉｌｉｔｙ ｏｆ ｐｒｏｄｕｃｉｎｇ ａｕｒｏｒａｓ． Ｗｈｅｎ ｉｔ ｉｓ ｇｒｅａｔｅｒ ｔｈａｎ ７．３８°，ｔｈｅ ｐａｒｔｉｃｌｅｓ ｗｉｌｌ ｂｅ ｂｏｕｎｄ ｉｎ ｔｈｅ ｇｅｏｍａｇｎｅｔｉｃ ｆｉｅｌｄ，ｆｏｒｍｉｎｇ ａ ｒａｄｉａｔｉｏｎ ｂｅｌｔ；（ⅳ）Ｗｈｅｎ ｏｔｈｅｒ ｃｏｎｄｉｔｉｏｎｓ ａｒｅ ｔｈｅ ｓａｍｅ，ｔｈｅ ｆａｒｔｈｅｒ ｔｈｅ ｃｈａｒｇｅｄ ｐａｒｔｉｃｌｅ ｔｈｒｏｗｉｎｇ ｐｏｉｎｔ ｉｓ ｆｒｏｍ ｔｈｅ ｅａｒｔｈ，ｔｈｅ ｇｒｅａｔｅｒ ｔｈｅ ｄｒｉｆｔ ｖｅｌｏｃｉｔｙ；ｔｈｅ ｇｒｅａｔｅｒ ｔｈｅ ｔｈｒｏｗｉｎｇ ａｎｇｌｅ，ｔｈｅ ｇｒｅａｔｅｒ ｔｈｅ Ｉｔｓ ｄｒｉｆｔ ｓｐｅｅｄ ｉｓ ａｌｓｏ ｇｒｅａｔｅｒ；（ⅴ）Ｆｏｒ ｐａｒｔｉｃｌｅｓ ｗｉｔｈ ｌｏｗ ｅｎｅｒｇｙ，ｔｈｅ ｇｕｉｄｅｄ ｃｅｎｔｅｒ ｔｒａｊｅｃｔｏｒｙ ｃａｎ ｗｅｌｌ ｒｅｐｒｅｓｅｎｔ ｔｈｅ ａｃｔｕａｌ ｔｒａｊｅｃｔｏ ｒｙ ｏｆ ｐａｒｔｉｃｌｅｓ ｉｎ ｔｈｅ ｆｉｒｓｔ－ｏｒｄｅｒ ａｐｐｒｏｘｉｍａｔｉｏｎ． Ｋｅｙ ｗｏｒｄｓ：ｅａｒｔｈ’ｓ ｍａｇｎｅｔｉｃ ｆｉｅｌｄ；ｔｈｅ ｐａｒｔｉｃｌｅ ｍｏｔｉｏｎ；ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ