地磁平靜期間磁層高能粒子非垂直地磁截止剛度研究

楚偉, 秦剛, 黃建平, 許嵩, 澤仁志瑪, 申旭輝

1 應急管理部國家自然災害防治研究院， 北京 100085 2 哈爾濱工業(yè)大學(深圳)， 深圳 518055

0 引言

太陽高能粒子是由太陽耀斑和日冕物質拋射爆發(fā)時導致行星際和磁層空間粒子環(huán)境增強的事件(Le and Zhang, 2017; Le et al., 2017; Zhao et al., 2018).如果太陽質子事件中的質子能量達到相對論的能量，這些相對論能量的質子可以穿越磁層，并與地球大氣發(fā)生核反應，從而使得地面儀器探測到相對論質子與地球大氣發(fā)生核反應的次級成分，這種現(xiàn)象稱為地面水平增強事件，簡稱Ground Level Enhancements(GLE)事件.這種事件主要出現(xiàn)在太陽活動高年(Le and Liu, 2020; Zhao and Le,2020).太陽高能粒子事件是偶發(fā)的，而銀河宇宙線是恒定存在的，只是其強度受太陽活動的調制.銀河宇宙線和太陽高能粒子進入磁層時會導致磁層高能粒子環(huán)境更惡劣，將對磁層空間的高能粒子環(huán)境產生重要的影響.高能質子能夠穿透電子元件防護層導致單粒子反轉，對宇航員的健康構成威脅，同時也會降低太陽能電池的壽命，影響極區(qū)短波通訊等，嚴重情況下將造成電子元器件損傷引發(fā)衛(wèi)星故障，甚至將衛(wèi)星徹底摧毀(Baker, 1998, 2001, 2002).因此研究行星際和磁層高能粒子耦合機制是日地空間物理的重點問題.

從上述公式中可以發(fā)現(xiàn)，當磁緯度為90°時，截止剛度達到最小.但上述公式中沒有給出截止剛度與粒子投擲角的相關關系，同時后續(xù)研究人員更多地將研究內容集中在垂直方向上，即天頂方向.Smart和Shea(2005)指出，在計算真實的截止剛度時，由于行星際高能粒子的各項同性，需要考慮來自不同方向粒子的截止剛度(Masarik and Reedy, 1995).同時行星際高能粒子的投擲角作為一個重要參量(可以研究粒子的加速以及絕熱聚焦效應)，對研究行星際和磁層高能粒子耦合具有至關重要的作用，因此針對行星際高能粒子的傳輸研究都需要考慮投擲角的擴散狀況.

Shea等(1965)使用數(shù)值模擬方法研究了地球表面超過300個點的垂直截止剛度，發(fā)現(xiàn)在南非、南大西洋和加那利群島附近，其截止剛度和地面觀測值具有高達15%的偏差.Daniel和Stephens(1966)研究了印度的HYDERABAD臺站附件的高能粒子截止剛度的方向依賴性.Bland和Cioni(1968)研究了有限方向(方位角方向15°間隔和極角方向10°間隔)的非垂直方向截止剛度，并利用兩個不同觀測位置(Trapani和Aire-sur-Adour)氣球測量試驗發(fā)現(xiàn)，測量和計算得到的截止剛度在上述兩個不同的觀測位置的均方差分別為4.8%和4.6%.Fanselow和Stone(1972)研究表明衛(wèi)星實際觀測結果比理論計算結果得到的垂直方向截止緯度范圍通常要低3°～5°，也就意味著垂直方向的截止剛度應該不是該點的最小地磁截止剛度.Cooke等(1991)針對變化磁場引起的截止剛度變化進行了研究，并針對計算過程中只考慮粒子剛度而忽略粒子方向的做法進行了指正.Bhattacharyya和Mitra(1997)研究了截止剛度隨著變化磁場的變化關系.Smart和Shea兩人針對截止剛度問題做了一系列的工作(Shea et al. 1965，1968； Smart，1999; Smart et al.,1999, 2000; Shea and Smart， 2001; Smart and Shea，2003a，b，2005)，針對地面、衛(wèi)星觀測和數(shù)值模擬存在的偏差展開了論述，并得到在進行衛(wèi)星觀測和數(shù)值模擬結果精確匹配過程中需要考慮不同方向粒子入射狀況的結論.Derome和Buénerd(2001)發(fā)現(xiàn)衛(wèi)星觀測到的剛度小于模擬獲得的截止剛度.Shimazu等(2006)研究了行星際磁場方向對截止緯度的影響.Dorman等(2008)研究了非垂直方向的截止剛度，發(fā)現(xiàn)當剛度為7～8 GV時，垂直方向和非垂直方向的截止剛度相差可能達到1 GV；Dmitriev等(2010)通過POES衛(wèi)星的觀測數(shù)據(jù)，確定了高能粒子截止緯度的橢圓形表達式.國內研究者通過數(shù)值模擬(甄杰和楚偉，2013； Chu and Qin, 2016)和觀測研究(朱邦耀，1979；樂貴明等，2003)發(fā)現(xiàn)地磁活動對截止剛度的影響.

由于使用解析表達式的垂直方向高能粒子截止剛度不能真實反映實際狀態(tài)的磁層高能粒子截止剛度，為了盡可能準確反映真實磁層的地磁截止剛度，需要研究真實磁場狀態(tài)下不同方向(投擲角)的截止剛度.

通過測試粒子的方式對不同能量(剛度)的粒子進行倒向追蹤計算，得到給定點的截止剛度.由于地球磁場的結構，使用此種方法需要計算包括回旋運動在內的高能粒子運動，同時需要追蹤大量的粒子運動，因此往往需要耗費大量的計算資源和計算時間，因此本文只給出地磁平靜期間的典型高能粒子截止剛度狀態(tài).

1 數(shù)值模擬方法

1.1 高能粒子剛度

高能粒子剛度作為描述高能粒子能量大小的物理量，其定義如下：R=mvc/Ze，其中R表示粒子剛度，m為高能粒子質量，v為粒子速度，c為光速，Z為電荷數(shù)，e為元電荷電量.剛度在數(shù)值上等于單位單核的能量，單位為V.截止剛度表征到達給定點的最小剛度.地磁粒子截止剛度可以描述地磁場對高能粒子的屏蔽作用.Le等(2006)推導了粒子剛度和能量的關系式，依據(jù)該關系式，我們可計算不同剛度粒子的能量.圖1為剛度與質子能量的關系圖，從圖上可以發(fā)現(xiàn)，兩者不是簡單的線性關系.

圖1 質子剛度Rigidity與能量Energy的關系圖Fig.1 The relationship between rigidity and energy for proton

Weygand和Raeder(2005)研究表明，使用目前通用的背景場模型都可以產生比較可靠的截止剛度模式，其中使用T96模型(Tsyganenko，1995； Tsyganenko and Stern，1996)作為外源場計算得到的截止剛度對應的緯度與衛(wèi)星觀測結果平均偏差為4.0°±1.4°，而使用T89模型(Tsyganenko，1989)作為外源場計算得到的數(shù)值為3.9°±2.4°，兩者的偏差并不明顯.由于本文研究平靜期間的截止剛度狀況，不采用T05(Tsyganenko and Sitnov，2005)模式，同時為了計算快捷方便，采用T89模式(Tsyganenko，1989)作為外源場.使用單粒子運動理論，結合地球磁層的內源背景場IGRF12模型以及T89外源場模型作為背景磁場，倒向追蹤粒子運動狀態(tài)，得到不同投擲角的粒子截止剛度，從中選擇最小剛度作為給定點的截止剛度.在計算過程中，只考慮背景磁場作用，忽略電場影響，因此運動過程中粒子的能量/剛度不會發(fā)生改變.同時由于是計算行星際粒子進入磁層，不考慮回旋運動和彈跳運動的投擲角損失錐造成的影響.

1.2 模擬條件的設定

模擬過程中采用四階龍格庫塔積分，初始高度設定為450 km高度.內邊界設定為1個地球半徑，外邊界由TS模式模型給出的磁層頂為邊界.同時為了保證計算效率最大化，排除捕獲粒子長時間做彈跳運動的影響，模擬過程中最大運動時長設定為1小時，最大運動距離設定為1000個地球半徑，最大運動步數(shù)為106，這樣既能確保計算結果的準確性，又能有效提高計算效率.

積分步長的設定將在一定程度上決定計算速度以及計算準確性.對于積分步長的設定目前主要存在兩種方法.第一種方法是使用龍格庫塔積分程序rkqs(Levin，1998)，設定誤差范圍內自適應的調節(jié)積分時間步長；第二種方法為使用粒子回旋周期的一定比例(通常情況下設定為粒子局地回旋周期的1/100，詳見Smart et al., 2000；Kress et al., 2010).在本文模擬過程中，在準確描述粒子運動前提下，盡量提高粒子積分效率.我們采用第二種方法，模擬過程中時間步長設定為粒子回旋周期的1/100.由于緯度對磁場的影響遠遠大于經度，因此模擬過程中緯度采用1°的間隔，經度采用10°的間隔.同時鑒于目前近地空間探測器的投擲角分辨率普遍大于5°，因此投擲角采用5°的間隔，粒子初始垂直平面方向運動速度進行均勻分配(Weygand and Raeder，2005).

我們在此使用的內部場模型是國際地磁參考場(IGRF)模型，該模型是地球主要磁場及其長期變化的標準數(shù)學描述.IGRF模型針對特定年份進行了標準化處理，反映了當時可用的最精確值.為了計算內源場和外源場貢獻，并進行計算過程中的坐標轉換以及磁層頂外邊界確認，我們使用Tsyganenko提供的軟件包GEOPACK-2008.以上所有有關磁場和輸入的參數(shù)均可從http:∥geo.phys.spbu.ru/～tsyganenko/modeling.html網站下載.

為了進一步提高計算效率，背景磁場使用插值進行計算，首先計算每一個格點的背景磁場，然后通過線性插值方法得到粒子運動位置的磁場強度.

2 數(shù)值模擬結果

圖2 計算過程中的截止剛度的半影區(qū)示意圖允許(白色)和禁止(黑色)在西經40°、北緯30°垂直方向，距離地面450 km高度處，剛度計算間隔為1 MeV的最大截止剛度、最小截止剛度和有效截止剛度示意圖.Fig.2 Schematic diagram of penumbra area of cut-off rigidityThe maximum, minimum and effective cut-off rigidity at intervals of 1 MeV for the calculation of the rigidity at a height of 450 km from the ground in the vertical direction of 40 degrees west longitude, 30 degrees north latitude for allowed (white) and prohibited (black).

2.1 不同投擲角截止剛度的對比

為了使研究結果具有普遍性，此研究中我們選取地磁平靜期間進行模擬.

圖3為截止剛度的全球分布圖像，左上圖為垂直方向的截止剛度分布圖像，右上圖為非垂直方向對應的截止剛度圖像，左下圖為全向的最小截止剛度圖像(包括不同方向的粒子截止剛度)，右下圖為垂直方向截止剛度和最小截止剛度的差值全球分布圖像.對比三個全球截止剛度的分布圖像，可以發(fā)現(xiàn)，通常所說的垂直方向(也就是天頂方向)一般意義上不是最小有效截止剛度對應的方向，尤其是在中緯度地區(qū)以及西經0～90°范圍內的差值最大，差值比例可能高達60%至70%.因此如果使用垂直方向的截止剛度作為地磁有效截止剛度，將會出現(xiàn)普遍高估的情形.從全球平均意義來看，垂直方向截止剛度將平均高估13%左右，最大偏差將達到70%，由此可見研究不同投擲角的高能粒子進入磁層的狀態(tài)對研究地磁截止剛度具有重要意義.從圖4可以發(fā)現(xiàn)，全方向計算得到的最小截止剛度和垂直方向計算得到的最小截止剛度偏差隨著緯度呈現(xiàn)蝴蝶狀分布形態(tài)，這也意味著中緯度地區(qū)的偏差比值達到最大.

2.2 最小截止剛度對應的投擲角的全球統(tǒng)計

圖5為最小有效截止剛度對應的投擲角的全球分布圖像(為了研究方便，我們使用185°表示垂直方向角度).從圖中可以發(fā)現(xiàn)，行星際高能粒子到達磁層的投擲角全球分布沒有明顯經緯度依賴性.

高能粒子的投擲角將影響高能粒子由行星際進入磁層，導致行星際高能粒子在磁層空間的分布投擲角依賴性明顯.從圖6我們可以發(fā)現(xiàn)，高能粒子由行星際空間進入磁層的最強有力的方向為沿著磁力線方向，也就是投擲角為0°或者180°時，此時高能粒子平行方向的能量巨大，可以直接沿著磁力線到達地球.同時由于研究的高能粒子能量非常高，天頂方向也是高能粒子進入磁層的一個重要通道.其他的投擲角方向分布符合中心點為90°的正態(tài)分布.

圖3 高能粒子截止剛度全球分布圖像(a) 垂直方向的截止剛度分布圖像； (b) 非垂直方向對應的截止剛度圖像； (c) 全向的最小截止剛度圖像(包括不同方向的粒子截止剛度)； (d) 垂直方向截止剛度和最小截止剛度的差值全球分布圖像.Fig.3 Global distribution image of cut-off rigidity(a) is vertical cut-off rigidity distribution image; (b) is non-vertical cut-off rigidity; (c) is omnidirectional minimum cut-off rigidity (including particle cut-off rigidity in different directions), and (d) is vertical direction Global distribution of the difference between the cut-off rigidity and the minimum cut-off rigidity.

圖4 截止剛度差值隨著經緯度變化圖像Fig.4 Cutoff rigidity difference changes with latitude and longitude

圖5 最小有效截止剛度對應的投擲角全球分布圖像(為了研究方便，我們使用185°表示垂直方向角度)Fig.5 The global distribution of the pitch angle corresponding to the minimum effective cut-off rigidity. (For convenience in the study, we use 185 degrees to represent the vertical direction)

圖6 最小截止剛度對應的投擲角占比柱狀圖(為了方便問題研究，我們采用185°表示垂直方向的角度，因此圖中出現(xiàn)185°的情形)Fig.6 The histogram of the pitch angle corresponding to the minimum cut-off rigidity (To facilitate the study of the problem, we use 185 degrees to represent the angle in the vertical direction, so the situation of 185 degrees appears in the figure)

一般假設銀河宇宙線為各項同性分布，因此由銀河宇宙線造成的進入磁層的高能粒子，其投擲角分布形態(tài)應呈現(xiàn)正態(tài)分布，此研究對磁層中高能粒子的投擲角分布具有重要作用.磁層中高能粒子的投擲角分布主要有三種可能的形態(tài)分布，各向同性分布、蝴蝶型分布以及薄餅型分布.此處獲得截止剛度對應的投擲角分布圖像與薄餅形態(tài)最接近，但是在兩端呈現(xiàn)上翹的形態(tài).

3 小結

高能粒子，尤其是行星際高能粒子進入磁層空間，將對磁層空間的粒子輻射環(huán)境產生影響.研究高能質子的截止剛度，尤其是中高緯度地區(qū)的高能粒子截止剛度，對研究行星際高能粒子在磁層空間形成的分布具有重要意義.本文使用單粒子數(shù)值模擬方法研究了近地空間高能粒子截止剛度與粒子投擲角的相關關系，研究發(fā)現(xiàn)：

(1)天頂方向或者垂直方向的截止剛度通常不是最小地磁截止剛度.

(2)最小地磁截止剛度對應的投擲角方向最大為沿著磁場方向，即0°方向；其次為天頂方向，也就是通常所說的垂直方向；然后為180°方向，即與磁場反方向.這種高能粒子進入磁層的投擲角分布狀態(tài)符合理論預期，因為沿著磁力線方向是粒子最容易進入磁層的方式，因此對應的截止剛度比較小.

(3)全球范圍截止剛度對應的投擲角分布符合兩端上翹的正態(tài)分布形態(tài)，不考慮兩端最大占比，中心位置出現(xiàn)在90°附近.

(4)通過地磁平靜期間的數(shù)值模擬發(fā)現(xiàn)，使用垂直方向的截止剛度對比最小截止剛度平均將高估13.17%，最大可能高估70%.

(5)不同經緯度高能粒子的截止剛度與投擲角不存在明顯關系.

本文使用了數(shù)值模擬方法進行高能粒子的運動模擬，由于使用的背景磁場和磁層頂邊界模型的不盡精確，得到的結果可能不能完全反映到達近地空間最小的高能粒子剛度，但能反映截止剛度的總體趨勢.為了得到更普適的截止剛度模式，在將來工作中，我們將利用數(shù)值模擬結合衛(wèi)星觀測數(shù)據(jù)，進一步修正不同經緯度的高能粒子截止剛度，以期做到盡可能精確可靠，為研究近地空間高能粒子環(huán)境變化提供精確的背景模型(Wang et al. 2020).隨著我國第一個電磁監(jiān)測試驗衛(wèi)星的成功發(fā)射，星上搭載的高能粒子探測器將對近地空間高能粒子展開軌道高度全球范圍、高投擲角分辨率的高能粒子探測數(shù)據(jù)，將有利于對此問題進行更加精確的研究.Rodger(2006)使用地基觀測數(shù)據(jù)和理論預期進行了對比，發(fā)現(xiàn)在磁擾動期間(Kp指數(shù)小于5的情況下)觀測和理論預期符合得很好，但是當Kp指數(shù)大于7，理論預期要比實際觀測結果大很多.我們后續(xù)將針對擾動條件下不同方向的截止剛度展開研究，屆時將對比平靜期間與擾動期間的變化磁場導致的截止剛度變化狀況.

致謝本研究得到了國家重點研發(fā)專項-地球物理探測衛(wèi)星數(shù)據(jù)分析處理技術與地震預測應用研究(2018YFC1503502-05)課題、地震活動期間ZH-1電離層掩星數(shù)據(jù)異常研究(ZDJ2019-03)課題以及ISSI-BJ(2019IT-33)課題的支持，同時得到了組內各位同事的鼎力幫助，在此感謝同事們的辛勤付出，使得本文得以完成.感謝IGRF提供的內源場源代碼以及Tsyganenko提供的外源場模型源代碼.同時在此對審稿老師提出的寶貴建議表示感謝.