ABSTRACT

Title of Dissertation: PROTON ENERGIZATION DURING
MAGNETIC RECONNECTION
IN MACROSCALE SYSTEMS

Zhiyu Yin
Doctor of Philosophy, 2024

Dissertation Directed by: Professor James F. Drake
Department of Physics

Magnetic reconnection is a widespread process in plasma physics that is crucial for the

rapid release of magnetic energy and is believed to be a key factor in generating non-thermal

particles in space and various astrophysical systems. In this dissertation, a set of equations are

developed that extend the macroscale magnetic reconnection simulation model kglobal to include

particle ions. The extension from earlier versions of kglobal, which included only particle elec-

trons, requires the inclusion of the inertia of particle ions in the fluid momentum equation. The

new equations will facilitate the exploration of the simultaneous non-thermal energization of ions

and electrons during magnetic reconnection in macroscale systems. Numerical tests of the prop-

agation of Alfvén waves and the growth of firehose modes in a plasma with anisotropic electron

and ion pressure are presented to benchmark the new model. The results of simulations of mag-

netic reconnection accompanied by electron and proton heating and energization in a macroscale

system are presented. Both species form extended powerlaw distributions that extend nearly three



decades in energy. The primary drive mechanism for the production of these nonthermal parti-

cles is Fermi reflection within evolving and coalescing magnetic flux ropes. While the powerlaw

indices of the two species are comparable, the protons overall gain more energy than electrons

and their power law extends to higher energy. The power laws roll into a hot thermal distribution

at low energy with the transition energy occurring at lower energy for electrons compared with

protons. A strong guide field diminishes the production of non-thermal particles by reducing the

Fermi drive mechanism. In solar flares, proton power laws should extend down to 10’s of keV,

far below the energies that can be directly probed via gamma-ray emission. Thus, protons should

carry much more of the released magnetic energy than expected from direct observations.

In Encounter 14 (E14), the Parker Solar Probe encountered a reconnection event in the he-

liospheric current sheet (HCS) that revealed strong ion energization with power law distributions

of protons extending to 500keV. Because the energetic particles were streaming sunward from an

x-line that was anti-sunward of PSP, the reconnection source of the energetic ions was unambigu-

ous. Using upstream parameters based on the data observed by PSP, we simulate the dynamics

of reconnection applying kglobal and analyze the resulting spectra of energetic electrons and

protons. Power law distributions extending nearly three decades in energy develop with proton

energies extending to 500keV, consistent with observations. The significance of these results for

particle energization in the HCS will be discussed.



RECONNECTION IS ALL YOU NEED

by

Zhiyu Yin

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2024

Advisory Committee:
Professor James F. Drake, Chair/Advisor
Dr. Michael M. Swisdak, Co-Advisor
Professor Thomas M. Antonsen, Jr.
Professor Alexander Philippov
Professor Christopher Reynolds (Dean’s Representative)



© Copyright by
Zhiyu Yin

2024



Preface

The journey to completing this dissertation has been both challenging and rewarding. My

interest in magnetic reconnection was sparked during a cosmic ray seminar in my early graduate

studies, where a talk by Jim led me to explore the acceleration mechanisms of plasma particles

further. This initial curiosity evolved into a deep commitment to understanding the energization

process of plasma particles during macro reconnection.

Throughout this research, I encountered numerous challenges, from debugging code to en-

suring numerical simulation stability, each serving as an opportunity for growth and learning. The

support and guidance of my advisors and colleagues have been invaluable, and their contributions

have greatly enhanced the quality of this work.

This dissertation is organized into several chapters. Chapter 1 provides an overview of

magnetic reconnection, motivations for building the new model, brief reviews of the electron ver-

sion of the new model, and a detailed outline of the thesis. Chapter 2 delves into the theoretical

background and testing of the new simulation model kglobal, while Chapter 3 provides the sim-

ulation results of the model. Chapter 4 gives the applications of kglobal. Chapter 5 offers the

conclusions of this research and discusses the future potential applications of the model.

I hope that this work will contribute to the field of plasma and heliophysics and inspire

future research in these areas.

ii



Zhiyu Yin

July 16, 2024

iii



Dedication

To my parents and my love, Ping Ren.

iv



Acknowledgments

It hath been a long and arduous quest,

At last, it finds its destined place of rest.

My heart doth overflow with joy and zest,

Emotions deep within, I do attest.

With gratitude, my soul is richly blessed,

For trials and triumphs that we have addressed.

Too many are the words I wish to test,

For this endeavor hath been our very best.

First of all, I want to extend my most sincere thanks to my research supervisors, Prof. Jim

Drake and Dr. Marc Swisdak. Both of you have guided me from confusion to clarity, helping

me grow from immaturity to a more seasoned researcher. You provided numerous opportunities

I had never dreamed of, broadening my horizons in the realm of real scientists. Your belief in me

as a scientist and physicist has been instrumental. Without your support, I couldn’t have come

this close to achieving my childhood dream. Thank you, Jim and Marc. This colorful journey we

shared will stay with me throughout my life.

I am deeply grateful to the members of my thesis and dissertation committee: Prof. Tom

Antonsen, Prof. Sasha Philippov, and Prof. Chris Reynolds. Your guidance, assistance, and

v



advice have been invaluable, greatly enhancing my understanding of physics and supporting my

future career.

I also want to express my gratitude to my collaborators, Prof. Mihir Desai and Dr. Tai

Phan, for providing me with the opportunity to apply my contributions to real-world scenarios.

Your collaboration and support allowed me to see the practical impact of my research, bridging

the gap between theoretical studies and real-world applications.

I would also like to thank my group members: Dr. Meriem Alaoui, Anna Fitzmaurice,

Haotian Da, Hanqing Ma, Dr. Dominic Payne, and Qianyi Ruan. Joining the group during the

Covid-19 pandemic was challenging, but your communication, advice, and support made all the

difference. Thank you for your academic companionship and constant encouragement.

My heartfelt thanks go to the alumni of my group: Prof. Michael Shay, Prof. Paul Cas-

sak, Prof. Yi-Hsin Liu, Dr. Kevin Schoeffler, Dr. Tak Chu Li, Dr. Joel Dahlin, Dr. Gareth

Roberg-Clark, and Dr. Harry Arnold. Your guidance, support, and shared experiences have been

invaluable. It is an honor to be part of this distinguished family.

My deepest thanks also go to Prof. Peter Shawhan, my academic advisor, and Prof. Eun-

Suk Seo, my former research supervisor. Both of you have been instrumental in helping me

navigate the challenging early days of my graduate studies. Your support and guidance eased my

transition and made me more resilient. I also want to thank the members of my former groups,

Dr. Leandro de Paula, Dr. Kichun Kim, Dr. Jayoung Wu, and Dr. Hongguang Zhang, for helping

me overcome those difficulties and for your unwavering encouragement.

I could not have endured the challenges of this Ph.D. journey without the unwavering sup-

port of my friends. My heartfelt thanks go to Dr. Mingshu Zhao for your steadfast support and

the countless joyful moments we shared. Dr. Yihang Wang, your kindness and sincerity were

vi



a constant source of strength, helping me navigate the darkest times. Jiateng Huang and Ziqi

Qin, your comforting presence and the relaxation you provided allowed me to find peace under

immense pressure. Shuyang Wang, your guidance and assistance in both academic and personal

matters were crucial to my progress. Wance Wang, our stimulating discussions, both academic

and otherwise, fostered my intellectual growth. Dr. Ziyue Zou, your readiness to assist in times

of difficulty was a great reassurance. Thanks to Haoying Dai for your daily interesting conversa-

tions and for listening and reasoning through my tough times. Thanks to Shuzhe Zeng for your

academic help, which has given me much to think about. I am deeply grateful to Kaixin Huang

for every encouraging conversation we had; your support provided the motivation I needed dur-

ing challenging times. Dr. Xuchen You and Dr. Xiaoyun Huang, your consistent support and

practical assistance were instrumental in overcoming obstacles. Dr. Yalun Yu and Xiangyu Li,

your care and the opportunities you shared were greatly appreciated, significantly enhancing my

experience. Dr. Haonan Xiong and Dr. Tong Wei, your companionship and thoughtful conversa-

tions alleviated my loneliness, bringing much-needed joy and perspective. Dr. Chongchong He,

Dr. Kungang Li, Dr. Yixu Wang, and Dr. Gong Cheng, your intriguing ideas and perspectives

were always a source of enjoyment and inspiration. Hua Huang and Xi Wang, thank you for your

amazing and heartwarming companionship as roommates.

I want to express my gratitude to Prof. Yi Wang, Dr. Hongyu Wang, Prof. Jiali Huang, Jie

Zhai and Dr. Jing Ao for the insightful conversations we had and the valuable information we

shared.

During my graduate studies in computer science at the Georgia Institute of Technology,

I had the pleasure of gaining incredible experiences working with Yun Zhang, Yisong Chen,

Xiaoyue Chen, Azhar Masani, Daniel Lesser, and Jae Hyeon Lee.

vii



I would also like to extend my heartfelt gratitude to Prof. Stéphane Coutu, my under-

graduate research supervisor. You provided exceptional opportunities to work in the lab and

introduced me to the fascinating field of physics and the universe. Before your guidance, my

dream of becoming a physicist was merely a lofty aspiration. Under your mentorship, I gained

not only fundamental knowledge and techniques but also developed the right attitude essential

for a physicist. Thank you for your unwavering support and invaluable lessons.

I feel fortunate to have attended Pennsylvania State University for my undergraduate stud-

ies, where I learned to build my life as an adult and consider my future seriously. During this

time, I extend my heartfelt thanks to Zangnan Yu for your companionship wherever you are. I

also express my gratitude to Shixin Zhao and Jiaming Zhang for our gatherings and for looking

after my family. Thank you, Terry Tang, for your career advice, and Runyu Mao, for our relaxing

conversations and your support. I appreciate Wei Cheng for your heartfelt company and Shuhan

Tian for your initial companionship when we came to the United States. Thanks to Jiayi Wang

and Hui Wang for all your help, communication, and warm invitations to spend time with your

family. Thank you, Boris Huang and Mike Wu, for your support and companionship during my

difficult time.

星垂平野阔，月涌大江流。作为北大附中的校友，每当回想我的高中，我还是难以

平复自己的内心。令公桃李满天下，何用堂前更种花。感谢张济达老师，为我指引物理

学的道路；感谢金钟鸣老师，既是良师也是益友；感谢王莹老师，在暗黑时期给予我唯

一的光明。这些恩情，铭刻在心，永不忘怀。孤帆远影碧空尽，唯见长江天际流。谢谢

你，给我长久陪伴的杨天骐；谢谢你，给我温柔交流的蔡宇坤，谢谢你；给我情绪抚慰

的郑峒；谢谢你，给我无尽帮助的王天雄；谢谢你，给我有趣灵魂的郝逸洲，谢谢你，

给我中肯鼓励的杨一博；谢谢你们，付震、黄无忌，谢谢你们，袁兆杰、原宇彤，谢谢

viii



你们，秦靖、李思恒，谢谢你，李家鑫，谢谢你，尹义平。

云山苍苍，江水泱泱，先生之风，山高水长。作为北京十五中分校的学生，也作为

曾经十五中初中的实习教师，无论是陶然亭畔的湖光山色，还是盆儿胡同的四季更替，

点滴细节，皆铭记于心。自蒙半夜传衣后，不羡王祥得佩刀感谢郭玉敏老师的深信不

疑，郭莎玲老师的关怀备至，吴晓光老师的雪中送炭，彭硕老师的鼎力相助，王会敏老

师的慧眼识珠，卢宇老师的仗义执言，范晓芸老师的细致入微，以及池晓波老师的教诲

如山。这些点点滴滴，皆为我成长的珍贵财富。还有高军老师，肖颖老师，于冬霞老

师，尽管不是你们的学生，但作为同事，也感谢您们的信任与照顾。还有我的挚友们，

荷笠带斜阳，青山独归远。把感谢送给肖迪，随着海风洒在星海的海滩；把感谢送给吴

慧超，随着夏夜留在丰台的车站；把感谢送给胡井暄，随着红魔呐喊在工体的国安；把

感谢送给王可，随着冬日驻在永恒的泰山；把感谢送给王梦妍，随着浓雾抚慰着天路的

草原；把感谢送给张帆，随着广西回到记忆的云南；把感谢送给李楠，随着单车畅游在

喧闹的二环；把感谢送给郜汝闯，随着哨声回响在费城的激战；把感谢送给周昊，送给

刘翔，送给董毅松，送给每一个你们，我的朋友们、同学们，还有我的学生们。

浮云游子意，落日故人情。我也十分有幸，不仅在学校中遇到了良师益友，还在这

一路上还结识了许多朋友，并得到了无数的帮助和支持。我想感谢陈晓航、谢小妍、何

雯，想感谢王丹，想感谢郭鹏，吴振媛，也想感谢揭妤。感谢不只学校中的田野，还有

你们这世界给予我的关切。

这里也想怀念我家的皮皮（狗），自从你离开后，我常常想起你，心里充满感激，

也常在梦中见到你。天堂可好？你一定也在想我吧。

最后，也要感谢我所有的家人，陈星姨妈、钱大行姨夫，陈震夏哥哥、夏蓉容姨

妈、陈健舅舅，陈越舅舅，韩超哥哥、郭碧薇姐姐、邸明明大姑、尹玫二姑、尹丹小

姑、韩晓风姑父、郭是忠姑父，还有阿姨马淑英。空床卧听南窗雨，谁复挑灯夜补衣，

ix



特别要感谢我的姥爷陈昌谦，姥姥陈嬗忱，奶奶邸金俊，爷爷尹旭。因为在落日时我曾

看到看到你们抹出粉色的天空，也在睡梦里探索了你们创造的树洞。

搴帷拜母河梁去，白发愁看泪眼枯。惨惨柴门风雪夜，此时有子不如无。感谢我的

父亲尹怡青和母亲陈晓千，让我把这一切献给你们。人生三十载，来美十一年，我收到

你们的所给予不只是生命，也有去探索世界的勇气、责任和努力。你们的一句“照顾好自

己，养活好自己”，是送给我成为一个无忧无虑孩子的祝福，而我也在努力成为一个能超

出你们期待的人。

君游东山东复东，安得奋飞逐西风。愿我如星君如月，夜夜流光相皎洁。最后的最

后，把我的感谢留给任凭。曾经的信件，代表我过去的思念；如今的文献，也写下和你

的诗篇；未来的百年，都会是美丽的晴天。感谢你这一路的鼓励、鞭策和陪伴，还有见

证着我的成长，还有还有我们那未来未知却有趣的每一天、每一刻、每一点。

敬谢诸般万象，亦谢此生千事。无论吾身何人，前途何往，皆因这一切，铸吾今时

之身，亦必成未来之才。千磨万击还坚劲，任尔东西南北风。在此留下吾之姓名，尹志

宇是也。

x



Table of Contents

Preface iv

Dedication iv

Acknowledgements v

Table of Contents xi

List of Tables xiii

List of Figures xiv

List of Abbreviations xvi

Chapter 1: Introduction 1
1.1 History of Magnetic Reconnection . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Historical Theoretical Development . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Multi X-line Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2 Particle Energy Gain in Magnetic Reconnection . . . . . . . . . . . . . . . . . . 10
1.2.1 Fermi Acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.2 Parallel Electric Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.2.3 Gradient-B Drift and Betatron Acceleration . . . . . . . . . . . . . . . . 24
1.2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

1.3 Limitations of Common Plasma Models . . . . . . . . . . . . . . . . . . . . . . 26
1.4 The Original kglobal Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.5 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Chapter 2: The Upgraded Computational Model for Particle Energization during Macroscale
Reconnection 34

2.1 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.1.1 kglobal with Particle Electrons . . . . . . . . . . . . . . . . . . . . . . 34
2.1.2 kglobal with particle electrons and ions . . . . . . . . . . . . . . . . . . 38

2.2 Alfvén Wave Propagation and the Growth of Firehose Modes in an Anisotropic
Plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Chapter 3: Simultaneous Proton and Electron Energization during Macroscale Mag-
netic Reconnection 49

xi



3.1 Simulation Model Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Chapter 4: Applications 73
4.1 Solar Flare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2 Heliospheric Current Sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.2.1 Simulations of Magnetic Reconnection . . . . . . . . . . . . . . . . . . 78

Chapter 5: Conclusion 81
5.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.2 Future Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.2.1 Earth’s Magnetotail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.2.2 Coronal Holes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.2.3 Anomalous Cosmic Rays . . . . . . . . . . . . . . . . . . . . . . . . . . 90

5.3 Future Developments and Enhancements of the kglobal Model . . . . . . . . . . 93

Appendix A: Theoretical Derivations 96
A.1 Derivation of the Fluid Ion Momentum Equation . . . . . . . . . . . . . . . . . 96
A.2 Energy conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Appendix B: Spectral Fitting 106
B.1 Kappa Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Appendix C: HCS Related Topics 109
C.1 PSP Instrumentation and Current Sheet Coordinate System . . . . . . . . . . . . 109
C.2 Details of HCS Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

Bibliography 116

xii



List of Tables

3.1 Parameters for Simulation Domains . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2 Density Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.3 Relative Energy and Spectra Parameters . . . . . . . . . . . . . . . . . . . . . . 72

5.1 Upstream temperature of electrons and protons in magnetotail . . . . . . . . . . 87

xiii



List of Figures

1.1 Sweet-Parker Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Petschek Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Multi-Island Magnetic Reconnection Model . . . . . . . . . . . . . . . . . . . . 11
1.4 Fermi Reflection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5 Contracting Magnetic Island . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6 Particle Energy Gain from the Merging of Two Islands . . . . . . . . . . . . . . 23
1.7 Fermi Acceleration vs Electric Fields . . . . . . . . . . . . . . . . . . . . . . . . 25
1.8 Spectra from Original kglobal . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.1 Alfvén Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.2 Firehose Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.1 Evolution of Particle Proton Density . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2 Energy and Firehose Parameter at Late Time . . . . . . . . . . . . . . . . . . . . 56
3.3 Time Evolution Spectra of Energetic Protons . . . . . . . . . . . . . . . . . . . 58
3.4 Spectra vs Different Simulation Domain Sizes . . . . . . . . . . . . . . . . . . . 60
3.5 Spectra vs Fraction of Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.6 Comparison of Initial and Late time Spectra between Protons and Electrons . . . 64
3.7 Spectra Fitting and Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.8 Linear-Linear Plot of Initial and Late-Time Spectra for Protons and Electrons . . 69

4.1 PSP Crossing Data and 2D image of Simulation . . . . . . . . . . . . . . . . . . 77
4.2 Comparison of Spectra from HCS Observations vs Simulations . . . . . . . . . . 79

5.1 Comparison of Initial and Late-Time Spectra between Protons and Electrons for
Magnetotail Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.2 Time Evolution Spectra with Pickup Ions in ACR Simulations . . . . . . . . . . 94

C.1 Time development of flux ropes in a simulation of the E14 HCS reconnection event113
C.2 The firehose parameter and the out-of-plane magnetic field Bz late in time from

the simulation in Fig C.1. The color bar indicates the values of α = 1− 4π(p∥ −
p⊥)/B

2 and Bz. Negative α values identify regions that are firehose unstable,
where the plasma has undergone significant energization. In these areas, the en-
ergetic particles influence the dynamics of flux ropes, suppressing reconnection
and the merging of islands. Within the flux ropes, Bz increases significantly
compared to its upstream value due to plasma compression in the reconnecting
current layer. This enhancement in Bz reduces the magnetic shear during flux
rope mergers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

xiv



C.3 Differential flux of electrons and protons at early and late times from the simula-
tion of Fig. C.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

xv



List of Abbreviations

ACR Anomalous Cosmic Rays

DSA Diffusive Shock Acceleration

E14 Encounter 14

GCR Galactic Cosmic Rays

HCS Heliospheric Current Sheet

KGLOBAL Kinetic Global

MHD Magnetohydrodynamic
MMS Magnetospheric Multiscale Mission

NASA National Aeronautics and Space Administration
NSF National Science Foundation

PIC Particle-In-Cell
PSP Parker Solar Probe

RHESSI Reuven Ramaty High Energy Solar Spectroscopic Imager

SCS Small-scale Current Sheet

xvi



Chapter 1: Introduction

1.1 History of Magnetic Reconnection

Plasmas, the fourth state of matter, constitute approximately 99.999% of all observable

matter, highlighting the importance of their study in understanding the fundamental nature of the

universe [1, 2]. Central to plasma physics is the process of magnetic reconnection, which enables

the transformation of magnetic energy into plasma energy by altering the magnetic topology

[3, 4]. Magnetic reconnection allows for the rapid reconfiguration of magnetic fields, facilitating

the release of magnetic energy into kinetic and thermal energy, as well as the acceleration of

charged particles into non-thermal distributions [5, 6, 7, 8], particularly power-law distributions

that extend over many orders of magnitude in energy [5, 7, 9, 8, 10]. This process is crucial

in the dynamics of a wide range of astrophysical systems, including the solar corona, planetary

magnetospheres, magnetic sectors of the solar wind, pulsar wind nebulae, and black holes [11,

6, 5, 12, 13, 7, 8, 14, 15, 16, 17]. In laboratory settings, magnetic reconnection is identified as

the primary cause of the sawtooth crash in fusion devices that leads to disruptions in confinement

[18, 19, 20, 21].

The significance of magnetic reconnection in both astrophysical and laboratory contexts

emerged from early observations of solar flares, which are sudden bursts of brightness above the

sun’s surface [22, 23, 24]. The first recorded solar flare, known as the ‘Carrington Flare,’ was

1



documented in 1859 by R. C. Carrington and R. Hodgson [25, 26, 27]. This event disrupted

telegraph systems and produced intense aurorae. Later, Giovanelli noted that solar flares are

associated with magnetic null points, regions where the magnetic field reverses direction [28].

Building on these observations, James Dungey proposed a model in which oppositely oriented

magnetic field lines separate and reconnect due to diffusion, forming a thin current sheet, and

applied reconnection to Earth’s magnetosphere [29], which laid the foundation for the modern

understanding of magnetic reconnection.

1.1.1 Historical Theoretical Development

1.1.1.1 The Sweet-Parker Model

The Sweet-Parker model is a foundational framework in the theory of magnetic reconnec-

tion, particularly within the context of astrophysical plasmas such as those found in the solar

atmosphere [30, 31]. Developed in the mid-20th century by P.A. Sweet and E.N. Parker, this

model provides a steady-state description of the reconnection process, where magnetic field lines

of opposite polarity are brought together, reconfigure, and release stored magnetic energy.

In the Sweet-Parker model (shown in Fig. 1.1), reconnection occurs in a two-dimensional

geometry where oppositely directed magnetic fields flow towards each other, forming a thin,

elongated current sheet at the boundary between the opposing magnetic fields. The magnetic

fields then reconnect across this current sheet. The resulting bent field lines expand downstream,

converting magnetic energy into kinetic and thermal energy, as well as accelerating particles.

Magnetic reconnection is self-driven since the expanding, reconnected field lines lead to a pres-

sure reduction at the magnetic x-line, which pulls in the upstream plasma and magnetic field to

2



sustain the process.

The Sweet-Parker model operates under several key assumptions that simplify the complex

dynamics of magnetic reconnection. First, the model assumes a steady-state scenario (shown in

Fig. 1.1), where the inflows and outflows of plasma (with velocities Vin and Vout, respectively)

and magnetic fields (Bin) are balanced, meaning the amount of material entering the current

sheet is equal to the amount leaving it. Additionally, the model presumes an elongated current

sheet structure (L ≫ δ), where both the kinetic energy of the inflowing plasma and the magnetic

energy of the outflow are neglected. Another important assumption is that resistivity (η) is only

considered within the current sheet, ignoring any effects outside this region. The model also

disregards 3D effects, focusing only on 2D dynamics.

First, the mass flux into the current sheet must equal the mass flux out, which can be

expressed as:

LVin ∼ δVout (1.1)

Similarly, the magnetic energy flux entering the current sheet is balanced by the kinetic energy

flux exiting it. This conservation of energy can be expressed as:

LVin

(
B2

in

4π

)
∼ δVout

(
ρV 2

out

2

)
(1.2)

By combining these equations, it is possible to show that the outflow velocity, Vout, scales with

the upstream Alfvén speed VA, defined as:

3



2L

2𝛿

𝐵𝑖𝑛

𝑉𝑖𝑛
𝑉𝑜𝑢𝑡

𝑉𝑖𝑛

𝐵𝑖𝑛

Figure 1.1: Cartoon diagram of the Sweet-Parker Model.

4



Vout ∼ VA ≡ Bin√
4πρ

(1.3)

Additionally, due to the 2D and steady-state assumptions, the out-of-plane electric field is

uniform inside and outside the layer, which gives:

VinBin

c
∼ ηJ (1.4)

Using Ampère’s law, the current density J is obtained as:

J ∼ c

4π

Bin

δ
(1.5)

The inflow velocity is then determined by balancing the resistive diffusion of magnetic flux with

the convective inflow, yielding:

Vin ∼
Dη

δ
(1.6)

where Dη ≡ ηc2

4π
is the resistive diffusivity.

The Sweet-Parker model provides a scaling law for the reconnection rate, which is propor-

tional to the inflow velocity and inversely proportional to the square root of the Lundquist number

S. The dimensionless reconnection rate r is given by:

r ∼ Vin

VA

∼ 1

S1/2
(1.7)

where the Lundquist number, defined as:

5



S ≡ LVA

Dη

(1.8)

measures the ratio of the resistive diffusion time scale to the Alfvén wave crossing time scale.

Typically, S is very large in astrophysical plasmas, often between 109 and 1020, leading to the

prediction of slow reconnection rates in the Sweet-Parker model. This slow reconnection rate is

inadequate to explain the rapid energy release observed in phenomena like solar flares. Thus,

while the Sweet-Parker model successfully demonstrates that reconnection occurs faster than by

diffusion alone, it requires additional mechanisms such as turbulence or 3D effects to account for

the fast energy dissipation observed in explosive astrophysical events [32, 33].

1.1.1.2 The Petschek Model and Advances in Collisionless Reconnection

The Petschek model, proposed by H.E. Petschek, advanced the theory of magnetic recon-

nection by suggesting that the presence of slow-mode shocks between the inflow and outflow

regions (shown in Fig.1.2) can substantially increase the rate of magnetic energy dissipation by

shortening the length of the current sheet (L ≫ L∗), which was taken to be macroscopic in the

Sweet-Parker model [34]. These shocks facilitate the conversion of magnetic energy into kinetic

and thermal energy. This model departs from the earlier Sweet-Parker framework by allowing

for faster reconnection with reconnection rate r ∝ 1
lnS

, which could potentially explain the rapid

energy release observed in various astrophysical phenomena.

6



Despite the theoretical appeal of the Petschek model, however, numerical simulations have

shown that this configuration only emerges as a solution of the MHD equations under certain

conditions. Specifically, the formation of slow-mode shocks and the resulting fast reconnection

rate are sustained only when there is a localized enhancement in the resistivity within the current

sheet [35]. In the absence of such localized resistivity, the system tends to revert to the slower

Sweet-Parker configuration, where the reconnection rate is limited by the resistive diffusion of

the magnetic field [36]. On the other hand, in the high Lundquist limit, magnetic reconnection

in large systems becomes highly turbulent [36, 37, 38]. Instead of a single, smooth reconnection

region, the current sheet can fragment and break up into multiple localized structures called

magnetic flux ropes. These flux ropes are small-scale magnetic loops or “islands” (also known as

plasmoids) of magnetic field lines that are twisted around each other. This fragmentation happens

due to instabilities such as the tearing mode instability, which allows the current sheet to break

into smaller regions of intense current [36, 39, 40, 41, 42, 43]. The presence of multiple flux ropes

suggests a more complex reconnection process, with many small regions where reconnection is

happening simultaneously rather than a single large reconnection outflow. This model will be

presented in the next subsection.

1.1.2 Multi X-line Model

Unlike the Sweet-Parker and Petschek models, the fragmentation of current sheets and the

development of plasmoids seems to be a common feature of systems undergoing reconnection

[36, 37]. Plasmoids have been detected both in solar flares [44, 45] and the Earth’s magnetotail

[46, 47, 48, 49, 50, 51]. Numerous simulations, ranging from fluid models [36, 52] to fully ki-

7



2𝐿

2𝐿∗

Figure 1.2: Cartoon diagram of the Petschek Model.

8



netic ones [39], have reported the occurrence of plasmoid formation. Plasmoids have become

an important component in theoretical models of magnetic reconnection. For instance, they can

increase the rate of magnetic energy release by reducing the effective length of the Sweet-Parker

current sheet and/or activating anomalous resistivity effects associated with small-scale plasma

dynamics. The scenario involving multiple plasmoids has been proposed to explain the acceler-

ation of electrons during reconnection events [53]. Recently, a magnetic reconnection plasmoid

model has been developed for flares associated with Sagittarius A*, the supermassive black hole

at the center of our Galaxy [54, 55].

This process results in what is known as multi-island reconnection, where multiple mag-

netic islands or structures emerge during reconnection. A schematic representation of multi-

island magnetic reconnection is provided in Fig. 1.3. The size of these islands ranges from

electron-scale to above the ion-scale. Neighboring islands with oppositely directed magnetic

fields give rise to localized reconnection events. Such small-scale reconnections can be initiated

by turbulence caused by the twisting and braiding of magnetic flux tubes. In large-scale recon-

nection events, such as those associated with solar flares, plasmoids created by current-driven

instabilities like the tearing mode instability occupy the magnetic field reversal region. In multi-

island magnetic reconnection, beyond the acceleration facilitated by the reconnection-induced

electric field, the contraction and merging of the magnetic islands also contribute to particle ac-

celeration [39, 56, 57, 58, 59, 60], as will be discussed further in Sec. 1.2.1.

9



1.2 Particle Energy Gain in Magnetic Reconnection

To investigate the various factors influencing electron energy evolution, we begin by con-

sidering the guiding-center approximation as described by Northrop [62]. For simplicity, we

specifically treat electrons here, although the extension to protons/ions is straightforward. The

change in energy, ϵ, of an electron in the guiding-center frame can be expressed as:

dϵ

dt
=

(
µ

γ

)
∂B

∂t
− e(v∥b+ vc + vg) · E (1.9)

Here, b = B/|B| is the unit vector along the magnetic field, µ = meγ
2v2⊥/2B represents the

magnetic moment, γ is the relativistic Lorentz factor, v∥ = v · b denotes the parallel component

of the velocity, and vc and vg are the curvature drift and gradient-B drift velocities, respectively,

given by:

vc =
v2∥b

Ωce

× κ (1.10)

vg =
v2⊥b

2Ωce

× ∇B

B
(1.11)

In these equations, Ωce = eB/γmec is the electron cyclotron frequency, and the magnetic curva-

ture is κ = b · ∇b.

Equation (1.9) can be reformulated in a more insightful manner:

dϵ

dt
= γmev

2
∥(uE · κ) + qE∥v∥ +

µ

γ

(
∂B

∂t
+ uE · ∇B

)
(1.12)

10



Figure 1.3: Left Panel: (a) A diagram illustrating volume-filling magnetic islands generated by
tearing instability in the magnetic field reversal region. (b) A PIC simulation depicting island for-
mation during magnetic reconnection. Here, jez represents the out-of-plane current density[53].
Right Panel: (c) A schematic illustration of neighboring magnetic islands with locally oppositely
directed magnetic fields (shown by arrows), undergoing magnetic reconnection[59]. The entire
figure is reprinted from the reference[61].

11



This expression highlights three distinct mechanisms by which particles can gain energy:

1. Fermi Acceleration: The first term, γmev
2
∥(uE · κ), describes Fermi acceleration, where

particles gain energy by reflecting off contracting magnetic field lines. The term uE · κ

indicates that this process is related to the curvature of the outflow magnetic fields. In this

scenario (shown in Fig. 1.4), as a magnetic field line contracts away from an X-line with

an Alfvénic outflow velocity CA, a charged particle moving along the field line experiences

a head-on collision with the moving field. For instance, if the curved field line moves to

the left with the Alfvén velocity and a charged particle travels rightward with a parallel

velocity v∥ = −v0, the particle reflects elastically from the field line. After the reflection,

the particle’s velocity increases to v∥ = v0 + 2CA, resulting in a potentially significant

energy gain.

2. Acceleration by Parallel Electric Fields (E∥): The second term, qE∥v∥, represents the

acceleration of particles due to electric fields that are aligned with the magnetic field. There

are two types of electric fields associated with this term: the kinetic-scale electric field, a

highly localized inductive parallel electric field that exists only at small scales because

highly mobile electrons tend to short-circuit these fields, limiting their range [63]; and the

macroscale electric field, generated by parallel electron pressure gradients, which plays a

significant role in particle heating[64, 65].

3. Betatron and Gradient-B Drift Acceleration: The third term, (µ/γ) (∂B/∂t+ uE · ∇B),

encompasses both betatron acceleration and gradient-B drift acceleration, both of which

are associated with the conservation of the magnetic moment µ.

• Betatron Acceleration: The component (µ/γ)∂B/∂t corresponds to betatron ac-

12



celeration. This process occurs as the magnetic field strength B changes over time,

with the particle’s perpendicular velocity adjusting to maintain the conservation of µ.

During magnetic reconnection, the perpendicular magnetic field decreases, and as a

result, the perpendicular energy ϵ⊥, becomes smaller.

• Gradient-B Drift Acceleration: The term (µ/γ)uE ·∇B represents gradient-B drift

acceleration. This mechanism takes place when a particle moves across magnetic field

gradients, causing B to vary as a function of position. Due to the conservation of µ,

the particle gains or loses energy as it drifts through regions with different magnetic

field strengths via the E×B drift.

13



These mechanisms accelerate particles by affecting different components of their energy.

The first and second mechanisms primarily enhance the parallel component of the energy. Par-

allel electric fields tend to generate beams of particles aligned with the magnetic field, which

may eventually thermalize due to various instabilities. Fermi acceleration, which is directly tied

to parallel energy, increases the parallel temperature T∥. In contrast, the betatron acceleration

and gradient-B drift mechanisms primarily increase the perpendicular energy in conjunction with

a rising magnetic field strength. Significant energy gain via betatron acceleration as well as

gradient-B drift typically requires a substantial increase in the magnetic field, which does not oc-

cur during magnetic reconnection. Thus, betatron acceleration and gradient-B drift are typically

energy sinks during reconnetion.

1.2.1 Fermi Acceleration

What is now known as Fermi acceleration was first proposed by Enrico Fermi in 1949 to

explain the origin of cosmic rays [66]. In his original formulation, Fermi proposed that charged

particles could gain energy through interactions with randomly moving ‘magnetic irregularities’

in the interstellar medium. When a particle collides head-on with a moving magnetic irregularity,

it gains energy, whereas a tail-on collision results in a loss of energy. Fermi argued that head-on

collisions are statistically more likely than tail-on collisions, leading to a net average acceleration

of the particles. The average gain in energy per collision can be expressed as:

δϵ ∝
(
V

c

)2

ϵ (1.13)

where ϵ is the energy of the particle and V is the velocity of the magnetic irregularity. Fermi

14



also considered that energetic particles have a finite lifetime, during which they can lose en-

ergy through collisions with ambient particles. The evolution of the distribution f ≡ f(ϵ, t) of

energetic particles can be described by the following equation:

∂f

∂t
+

∂

∂ϵ
(ϵ̇f) = − f

T
. (1.14)

Here, T represents the characteristic lifetime of the energetic particles, and ϵ̇ is the rate of energy

gain which, for interactions with magnetic irregularities, is given by:

ϵ̇ = ϵ

(
V

c

)2
1

τ
, (1.15)

where τ is the average time between interactions. Assuming a steady state (∂/∂t = 0), the above

equation simplifies to:

∂

∂ϵ
(αϵf) = − f

T
, (1.16)

where α = (V/c)2/τ . The solution to Eq. 1.16 takes the form of a power law:

f(ϵ) ∝ 1

ϵ1+(αT )−1 , (1.17)

which is the functional form required to explain the observed cosmic ray spectrum.

Fermi’s original concept was later adapted to describe acceleration at astrophysical shocks,

leading to the development of what is now known as Diffusive Shock Acceleration (DSA) [67,

68, 69, 70]. A shock propagating through a plasma compresses and slows the upstream plasma as

it crosses the shock. Particles undergo acceleration as they reflect back and forth across the shock

15



[71, 72]. This process is a specific case of Fermi acceleration, often referred to as first-order

Fermi acceleration since the rate of energy gain of a particle is proportional to the shock velocity.

It is the leading model for the acceleration of galactic and intergalactic cosmic rays [73].

In addition to the specific case of first-order Fermi acceleration at shocks, Fermi’s gen-

eral concept has broader applications in particle acceleration processes. The second-order Fermi

acceleration mechanism occurs in more turbulent environments, where particles gain energy

through repeated interactions with randomly moving magnetic irregularities [74, 75]. Unlike

the first-order process, where energy gain is proportional to the velocity of the shock, the second-

order process is proportional to the square of the velocity (V/c)2. This makes second-order Fermi

acceleration less efficient but still significant in environments dominated by turbulence.

Fermi acceleration, whether first-order or second-order, fundamentally relies on the con-

cept of particles gaining energy through interactions with moving magnetic structures. These in-

teractions can be modeled in various astrophysical settings, such as supernova remnants [76, 77],

the solar wind [78], and the interstellar medium [79, 80], where they contribute to the acceleration

of particles to high energies, forming the power-law energy spectra observed in cosmic rays.

1.2.1.1 Fermi Acceleration in Reconnection

During magnetic reconnection, magnetic islands or flux ropes are formed as elongated

magnetic structures that contract at both ends with Alfvén speed, CA (See Fig. 1.5). Particles

trapped within these islands can reflect multiple times between the ends, gaining energy with

each reflection. If the particle’s velocity is significantly greater than the Alfvén speed, which is

typically the case for electrons, the energy gain per reflection can be substantial.

16



During magnetic reconnection with adiabatic invariance, the conservation of parallel action

can be applied to determine the particle energy gain in an elongated magnetic island. The parallel

action is given by:

J = v||l (1.18)

where l is the circumference of the ellipse, whose expression was given by Ramanujan[81]:

l = π[3(L+ ω)−
√

(3L+ ω)(L+ 3ω)] (1.19)

In Fig. 1.5, the circumference of the ellipse of the island for L ≫ w is:

l ≃ (3−
√
3)πL ≃ 4L (1.20)

And thus, with the contracting rate L̇ = CA, the derivative of l is as follows,

l̇ ≃ 4L̇ = −4CA (1.21)

On the other hand, the time derivative of J is:

dJ

dt
= 0 = v̇||l + v||l̇ (1.22)

which leads to the change in the parallel velocity:

v̇∥ = −v||
l̇

l
= −v∥

(−4CA)

4L
= v||

CA

L
(1.23)

17



As a result in Eq.1.23, particle energy gain in a contracting island follows a first-order

Fermi process.

18



𝐶!

𝑣" + 2𝐶!

𝑣"

Figure 1.4: Cartoon diagram of a charged particle reflecting from an outflow magnetic field line at
the Alfvén speed. The particle velocity increases by 2CA and particle’s final velocity is v0+2CA.

19



Another perspective on the energy increase can be understood through the merging of two

magnetic islands. As illustrated in Fig. 1.6, two islands, each with an initial radius ri, merge to

form a larger island with a final radius rf . Initially, the path length is approximately li ≃ 2(2πri),

and after the merger, the final path length becomes lf = 2πrf .

During this process, the total area enclosed by the islands remains conserved [82], leading

to the equation:

2(πr2i ) = πr2f (1.24)

from which we obtain:

rf =
√
2ri (1.25)

As a result, the ratio of the final to initial path lengths is given by:

lf
li

=
rf
2ri

=
1√
2
. (1.26)

By conserving the parallel action, we can then deduce the velocity scaling:

vf =
√
2vi (1.27)

Thus, the energy (∝ v2) doubles during the merging of the two islands.

20



2𝐿

2𝜔
𝐶𝐴𝐶𝐴

Figure 1.5: Cartoon diagram of a magnetic island contracting at the Alfvén speed.

21



1.2.2 Parallel Electric Fields

A seemingly straightforward mechanism for particle acceleration in reconnection regions

involves parallel electric fields. Two types of parallel electric fields will be explored here, large-

scale parallel fields associated with the gradient of electron pressure, and kinetic-scale parallel

fields. The role of these fields in particle acceleration has been the subject of numerous studies

[8, 83, 64].

In a previous study exploring the strong guide field regime [84], the kinetic-scale paral-

lel electric field (E∥) drove essentially all of the electron energy gains, yet it failed to generate

an energetic non-thermal component, and in particular did not produce a power-law distribu-

tion. Fig.1.7 shows the cumulative energy gain of electrons resulting from different accelera-

tion mechanisms—including parallel electric fields (E∥), Fermi reflection, and betatron acceler-

ation—across varying guide field strengths (bg), highlighting that E∥ plays an important role in

driving electron energy gains in strong guide field regimes. In the scenario with a weak guide

field (as shown in Fig.1.7(a)), the dominant contributions come from Fermi reflection (positive)

and betatron acceleration (negative), which partially offset each other, resulting in the bulk of

the electron energy gain, while the effect of parallel electric fields is negligible. In the simula-

tion with bg = 4 (as shown in Fig.1.7(c)), all terms, including parallel electric fields (E∥), have

insignificant contributions to the electron energy gain. In the intermediate case with bg = 1 (as

shown in Fig.1.7(b)), both Fermi reflection and E∥ play significant roles. It is important to note

that the total electron energy gain is approximately 50% higher in the bg = 0.2 case compared

to the bg = 4.0 case. In summary, Fig.1.7 reveals how the contribution of E∥ compares to other

mechanisms, showing that the influence of E∥ only remains dominant for higher bg (under higher

22



𝑟𝑖 𝑟𝑖

𝑟𝑓

Figure 1.6: Cartoon diagram of the merger of two islands with identical radii and magnetic field
strengths.

23



guide fields, the reconnection produce almost no energetic electrons [84, 85]), with implications

for understanding the efficiency of electron energization during magnetic reconnection. In a later

study [86] it was found that for low guide fields the electric field acts as an energy sink instead of

a driver for highly energetic electrons and primarily contributes to heating the lower energy range.

(Note: In this dissertation, energization refers to accelerating particles to non-thermal energies,

often leading to the formation of a power-law distribution, while heating refers to accelerating

particles thermally, resulting in a Maxwellian distribution.)

On the other hand, large-scale parallel electric fields, extending beyond localized boundary

layers, enhance the process described in Fig.1.4 by confining electrons within the reconnection

exhaust, allowing them to experience multiple Fermi accelerations[87, 64, 65, 88]. Furthermore,

the study [85] shows that the heating by electric fields remains nearly constant for electrons under

different guide field strengths at late times.

1.2.3 Gradient-B Drift and Betatron Acceleration

The betatron mechanism and gradient-B drift are two fundamental processes that play a

significant role in the acceleration of charged particles in both laboratory[89, 90] and astrophysi-

cal plasmas[91]. The magnetic moment of a particle µ is an adiabatic invariant. This invariance

holds as long as the temporal scales involved are much longer than the particle’s cyclotron period.

Consequently, changes in the magnitude of the magnetic field directly influence the evolution of

the particle’s perpendicular energy ϵ⊥. Betatron acceleration occurs when the magnetic field

strength B increases over time, causing the particle to adjust its perpendicular velocity to main-

tain a constant magnetic moment.

24



Figure 1.7: (a)–(c) Cumulative electron energy gain from different acceleration mechanisms for
various guide field strengths (bg). The red, blue, and magenta curves represent contributions from
Fermi reflection, parallel electric fields (E∥), and betatron acceleration, respectively. The black
dashed line shows the total energy gain. Panels (a), (b), and (c) correspond to bg = 0.2, bg = 1.0,
and bg = 4.0, respectively. (d) Normalized total electron energy gain as a function of guide field
strength (bg). Symbols represent different contributions: squares for total energy, diamonds for
Fermi reflection, and stars for E∥. The red line indicates a fit of the form ∼ [1 + 4b2g]

−1. Reprint
from [84].

25



During magnetic reconnection, the release of magnetic energy generally leads to a reduc-

tion in the magnetic field strength. Given that betatron acceleration and gradient-B drift rely on

an increase in B to effectively impart energy to particles, these mechanisms are not expected

to play a significant role in the context of the simulations presented in this dissertation[92, 86].

Therefore, the impact of betatron acceleration and gradient-B drift is considered negligible in the

modeled scenarios.

1.2.4 Conclusion

Based on the discussions above, the dominant drive mechanism of particle energy gain

during magnetic reconnection is Fermi reflection in growing and merging magnetic flux ropes

[39, 56, 92, 93, 86, 94]. Magnetic reconnection produces bent magnetic field lines that form an

Alfvénic exhaust, carry energy away from the x-line, and transfer energy from the field to the

surrounding plasma [31, 95]. Particles gain energy through repeated reflections in contracting

magnetic field lines, thereby leading to the observed power-law tails in particle energy distribu-

tions [39, 92, 86]. The rate of energy gain from Fermi reflection is proportional to a particle’s

energy and is therefore greatest for the most energetic particles.

1.3 Limitations of Common Plasma Models

Both particle-in-cell (PIC) and magnetohydrodynamic (MHD) [36, 96, 97] models have

limitations that compromise their utility for studying non-thermal particle energization.

PIC is a widely used computational technique for simulating plasma dynamics in both lab-

oratory and astrophysical contexts. This method is particularly powerful for studying systems

26



where the behavior of individual charged particles and their interactions with electromagnetic

fields are crucial, such as in magnetic reconnection, shock waves, and wave-particle interactions.

PIC simulations correctly model kinetic physics by resolving kinetic scales [98, 99, 100, 101, 86,

94]. “Macro-particles” that represent a collection of real particles move through the simulation

domain. The electric and magnetic fields are evaluated at the gridpoints surrounding the cells

containing macro-particles and in turn accelerate the particles via Lorentz forces. This model is

ideal for simulations at kinetic scales, but it has limitations at macroscales due to the computa-

tional expense, which prevents the simulation of systems significantly larger than kinetic scales.

The spatial scales associated with solar flares can exceed 104 km, while the Debye length—a

measure of the scale of kinetic effects—can be on the order of centimeters, representing a scale

separation of 1010. Hence, computational costs drastically limit the ability of PIC simulations to

simulate the extensive range of scales occurring in magnetic reconnection events. A consequence

of this inadequate separation of scales is the demagnetization of energetic electrons and ions in

PIC simulations of reconnection, which in turn may inhibit particle energy gain and prevent the

formation of the extended power laws that are observed in nature [86, 94].

MHD combines the principles of fluid dynamics and electromagnetism to describe the be-

havior of electrically conducting fluids such as plasmas. The MHD model is particularly useful

for studying large-scale plasma dynamics in astrophysical and laboratory settings, where the

fluid approximation is valid and kinetic effects are less dominant. Hence, in MHD, small-scale

motions are ignored. While capable of representing macroscopic events, MHD simulations fre-

quently assume the plasma is thermal (i.e., represented by a single temperature) and hence by

definition cannot produce non-thermal particles. The electric and magnetic fields from MHD

models have been used to explore non-thermal energization of test particles. However, these test

27



particles do not interact self-consistently with the fields so energy is not conserved, which limits

the test particle model’s fidelity [102].

1.4 The Original kglobal Model

To address these challenges, a new formalism was developed that permitted the first self-

consistent simulations of electron energization during magnetic reconnection in a macroscale

system [103, 88, 85]. The model was named kglobal, standing for ”kinetic global” and originally

did not consider ion energization. To clearly distinguish it from the version that does consider

ions (and serves as the main focus of this dissertation), the version that explores electron acceler-

ation exclusively is referred to as the original version. The theoretical background and details of

the original version of the simulation model will be presented in Sec 2.1.1. The original kglobal

model combines aspects of both the PIC and MHD descriptions to address electron acceleration

in large-scale astrophysical reconnection [103, 88]. Particle electrons are distributed on the MHD

grid, as in PIC models, but are evolved with MHD fields using the guiding-center approximation,

where all kinetic scales, including the Larmor radii of particles, are ordered out of the equations.

Notably, kglobal includes the effects of Fermi reflection within evolving magnetic flux ropes and

the pressure anisotropy that feeds back on the MHD fields to reduce the tension force driving

magnetic reconnection. Magnetic reconnection simulations with kglobal produced, for the first

time, extended power-law spectra of non-thermal electrons that aligned with existing observa-

tional data [85]. This notable result can be seen in Fig. 1.8, which presents a log-log plot of the

differential number density versus the normalized energy. The results show the temporal evolu-

tion of the electron distribution in a simulation with a modst guide field, Bg/B0 = 0.25, where

28



an initial Maxwellian spectrum develops an energetic tail that hardens over time, eventually sta-

bilizing into a power law. The figure also demonstrates that smaller guide fields lead to harder

spectra that extend to higher energies. Additionally, the late-time distributions of energetic elec-

trons indicate that while the power-law indices vary weakly with the domain size, the maximum

energy of the electrons increases with larger system sizes.

29



The implications of the original kglobal simulations were noteworthy. However, observa-

tions show that magnetic reconnection can simultaneously produce non-thermal electrons and

non-thermal protons [5, 104]. Hence, developing a self-consistent simulation model is vital to

elucidate the mechanisms of ion and electron energization during magnetic reconnection and the

resulting partitioning of energy between the two species. In Chapter 2, we extend the kglobal

equations to include particle ions, which requires the inclusion of the inertia of the particle ions

in the bulk fluid ion momentum equation. The resulting equations will facilitate the exploration of

the self-consistent energization of both ions and electrons during magnetic reconnection. Based

on these equations, the current version also includes particle protons and produces non-thermal

spectra for both species. The proton distribution features a non-thermal tail that extends over

more than three decades in energy.

1.5 Outline of Thesis

In Chap. 2, we will revisit the original kglobal model in detail, show additions made in the

upgraded model and present tests of the upgraded model that verify its capabilities. The model

consists of an MHD backbone with macroparticle protons and electrons distributed on the MHD

grid. The guiding center equations with the MHD electric and magnetic fields describe the mo-

tion and energy gain of the macroparticle protons and electrons. These particles feedback onto

the MHD fluids and fields and conserve energy. We show kglobal correctly simulates Alfvén

waves in the presence of a pressure anisotropy, which is crucial for modeling Fermi reflection

since pressure anisotropy can reduce magnetic field tension and the associated Fermi drive mech-

anism. Specifically, we demonstrate that the model reproduces the correct growth rate of the

30



Figure 1.8: A log-log plot of the differential number density versus the normalized energy. Panel
(a) shows the temporal evolution of the electron distribution in a simulation with Bg/B0 = 0.25.
The initial Maxwellian spectrum develops an energetic tail that hardens over time and eventually
stabilizes into a power law. Insert in Panel a shows the late time spectra for different guide fields.
Smaller guide fields result in harder spectra that extend to higher energies. Panel b shows the late-
time distributions of energetic electrons in a simulation with Bg/B0 = 0.25 for various values
of Sν (effective system size) reveal that the power-law indices of the spectra vary weakly with
the domain size. The notable difference is that as the size of the system increases, the maximum
energy of the protons also increases. Reprinted with permission from [85]

31



linear firehose instability, which onsets when the local field-line tension goes to zero, using the

same techniques as used in the original electron model [103].

In Chap. 3, we present results from 2D magnetic reconnection using the upgraded kglobal

model. Reconnection begins from particle noise and produces multiple islands (or flux ropes

since they consist of an axial magnetic field wrapped by the inplane magnetic field) that merge,

forming larger islands and accelerating electrons and protons. The pressure anisotropy from the

plasma macroparticles becomes large enough to eliminate the magnetic tension force along out-

flow exhausts and within islands when particle acceleration is strong. Both species form extended

powerlaw distributions that extend nearly three decades in energy. The primary drive mechanism

for the production of these nonthermal particles is Fermi reflection within evolving and coalesc-

ing magnetic flux ropes. While the powerlaw indices of the two species are comparable, the

protons overall gain more energy than electrons and their power law extends to higher energy.

The power laws roll into a hot thermal distribution at low energy with the transition energy oc-

curring at lower energy for electrons compared with protons. A strong guide field diminishes the

production of non-thermal particles by reducing the Fermi drive mechanism.

In Chap. 4, we explore the complex processes of charged-particle acceleration during mag-

netic reconnection, focusing on solar flares and the heliospheric current sheet (HCS). The chapter

begins with a discussion of solar flares, highlighting the role of magnetic reconnection in trans-

forming magnetic energy into particle energy and the formation of non-thermal particle popu-

lations. Observational data from instruments like RHESSI are discussed to provide context for

the particle acceleration mechanisms in solar flares. The chapter then shifts focus to the HCS,

the largest current sheet in the solar system, analyzing data from the Parker Solar Probe (PSP)

during Encounter 14. The detection of sunward-directed reconnection exhausts and the observa-

32



tion of proton energy spectra are discussed. Using the kglobal model, simulations are performed

to investigate the underlying physical mechanisms driving extreme proton energization observed

in the HCS. The results demonstrate that merging magnetic flux ropes during reconnection can

accelerate protons to energies far exceeding the available magnetic energy per particle. This

chapter combines observational insights with advanced simulations to deepen our understanding

of particle acceleration in large-scale reconnection events.

In Chap. 5, we present conclusions, discuss the results in the context of previous efforts

to explore reconnection-driven plasma particle acceleration, and suggest future extensions and

potential applications of the kglobal model.

33



Chapter 2: The Upgraded Computational Model for Particle Energization dur-

ing Macroscale Reconnection

The material in this chapter has been adapted from [105] with permission from the authors.

2.1 Model Equations

This chapter presents the upgraded kglobal model. To ensure a comprehensive understand-

ing, I will first review the governing equations of the original kglobal model before explaining

the enhancements made in the updated version. The upgraded model has been tested and found

to accurately represent the propagation of an Alfvén wave in the presence of ion and electron

pressure anisotropy, as well as the linear firehose instability. These benchmarks are crucial be-

cause bent and reconnected magnetic field lines propagate at the Alfvén speed, and reconnection

is moderated by the reduction of field-line tension as a plasma nears firehose marginal stability.

2.1.1 kglobal with Particle Electrons

The original model includes three distinct plasma species: the ion fluid (number density ni),

the electron fluid (nef ), and particle electrons (nep) [103, 88]. No conversion between species is

allowed. The ratio of the particle to electron fluid density is a free parameter. However, electron

heating and energization were shown to be insensitive to this ratio, as discussed in reference [85].

34



The particle electrons are governed by the guiding-center equations [62]. Their perpendicular

velocity is determined by conservation of the magnetic moment:

µep =
p2ep,⊥
2B

= const. (2.1)

where pep is the momentum and B is the magnitude of the local magnetic field. Their perpen-

dicular motion is given by the usual MHD E×B drift while they stream parallel to the ambient

magnetic field at their parallel velocity. The variation in the parallel velocity reflects contributions

from Fermi reflection, magnetic mirroring, and the large-scale parallel electric field

dpep,∥
dt

= pep,∥vE · κ− µe

γe
b · ∇B − eE∥ (2.2)

with vE the E × B drift velocity, b = B/B a unit vector in the direction of the magnetic

field, κ = b · ∇b the field curvature, and γe the relativistic Lorentz factor. The guiding-center

electrons move through the MHD grid and, crucially, do not need to have a thermal distribution.

Their gyrotropic pressure tensor feeds back on the MHD momentum equation.

The ion fluid forms an MHD-like backbone that satisfies the usual continuity equation

∂ni

∂t
= −∇ · (nivi) (2.3)

and a modified momentum equation

mi
∂(nivi)

∂t
= −∇ · Ti − (∇ · Tef )⊥ +

J×B

c
− (∇ · Tep)⊥ + eniE∥b (2.4)

35



where the ion stress tensor with inertial contributions is

Ti = PiI+minivivi (2.5)

with Pi the ion scalar pressure and I the unit tensor. The corresponding electron fluid tensor is

Tef = PefI+menefγefv
2
ef,∥bb (2.6)

with Pef the electron fluid scalar pressure. Finally, Tep is the particle electron gyrotropic stress

tensor, including inertial contributions,

Tep = Tep∥bb+ Pep⊥(I− bb) (2.7)

where Tep∥ is the component of the stress tensor along the magnetic field B, and Pep⊥ is the usual

perpendicular pressure,

Pep⊥ =

∫
dpe

p2e⊥
2meγe

f (2.8)

where in the frame drifting with vE = cE × B/B2 there are no perpendicular flows so f =

f(x, pe∥, pe⊥, t). Tep∥ includes the mean parallel drifts of the particle electrons,

Tep∥ =

∫
dpe

p2e∥
meγe

f (2.9)

with pe∥ being the electron particle parallel momentum.

Note that in the ion momentum equation the electron fluid and particle inertia are neglected

36



compared to the ion inertia on the left side of Eq. (2.4), while the electron fluid and particle

pressures are retained on the right side of this equation. This result is based on an ordering in

which the ion velocity scales as the Alfvén speed and the electron particle and fluid thermal

velocities and parallel streaming velocities are of the order of the electron Alfvén speed, CAe =√
mi/meCA, where the Alfvén speed is defined as CA =

√
B2/(4πmini). In this ordering the

electron momentum scales as meCAe =
√

me/mimiCA and so is much smaller than that of the

ion fluid. In contrast, the electron pressure is meC
2
Ae = miC

2
A, which is the same as the ion

pressure. Thus, the pressure terms on the right side of Eq. (2.4) must be retained even though

the electron inertia in this equation is neglected. Finally, the ions follow an adiabatic equation of

state

d

dt

(
Pi

nγ
i

)
= 0 (2.10)

where the adiabatic index γ is taken to be 5/3.

The number density of the electron fluid is given by

nef = ni − nep (2.11)

while the parallel motion follows from the requirement of vanishing parallel current

nefvef,∥ = nivi,∥ − nepvep,∥ (2.12)

The zero current condition follows from the constraint on the current from Ampère’s law

Vd

CA

=
J

neCA

∼ Bc

4πneCAL
∼ di

L
≪ 1 (2.13)

37



where Vd is the drift speed associated with the current, di is the ion inertial scale length and L is

the characteristic magnetic scale length. In the kglobal model all kinetic scales are ordered out,

which yields the zero current condition in Eq. (2.12). The fluid electrons also follow an adiabatic

equation of state

d

dt

(
Pef

nγ
ef

)
= 0 (2.14)

The magnetic field is advanced via Faraday’s law,

∂B

∂t
= −c∇× E⊥ (2.15)

with

E⊥ = −1

c
vi ×B (2.16)

Finally, E∥ is produced by magnetic-field-aligned gradients in the electron pressure and is given

by [88]

E∥ = − 1

e(nep + nef )
(b · ∇ · Tep + b · ∇Tef ) (2.17)

Thus, Equations (2.1)-(2.17) constitute the complete set of equations for the kglobal model with

particle electrons.

2.1.2 kglobal with particle electrons and ions

In the upgraded version of kglobal, particle ions are introduced as a new species and, as is

the case for the particle electrons, are treated in the guiding-center limit. Similar to the original

model, the populations of the four plasma species do not change after initialization. The per-

38



formance differences arising from various values of the free parameter—the ratio between the

number density of particle ions (nip) and the total number density of ions (nif + nip)—will be

discussed in Section 3.2. The basic ordering of the electrons and the equations describing their

dynamics are basically unchanged from the original model. Minor modifications to the charge

neutrality and zero current conditions are discussed below, but the details of the electron dynam-

ics that were discussed earlier are not repeated here.

Due to the small mass of the electrons, their inertial terms in the fluid ion momentum

equation were neglected in the original kglobal equations. That simplification is not possible

for the particle ions in the upgraded model, which complicates the derivation and the final form.

However, the perpendicular motion of the fluid and particle ions is still governed by the same

E×B drift. As a consequence, the generic form of the ion momentum equation, can be written

as

∂

∂t
mi(nifvif + nipvip) = S (2.18)

with vif and nif the fluid ion velocity and density, vip and nip the particle ion velocity and

density, and S the usual momentum driver, including convective inertial terms. The differences

in the inertia between the fluid and particle ions arises solely from the parallel motion so the total

momentum equation takes the form

∂

∂t
minivif +

∂

∂t
minip(vip∥ − vif∥)b = S (2.19)

where ni = nif + nip. The inertial terms parallel to B are easily evaluated, which leaves an

equation for the fluid ion momentum that includes the inertia of the particle ions. The details of

39



the evaluation of the parallel inertial terms are found in the Appendix.

Most of the other equations in Section 2.1.1 remain unchanged, aside from minor notation

adjustments (e.g., ni → nif + nip). Specifically, Equations (2.1) and (2.2) governing the motion

of the particle electrons are unchanged. Equations (2.3) and (2.10) become equations for the ion

fluid density,

∂nif

∂t
= −∇ · (nifvif ) (2.20)

and ion fluid pressure,

d

dt

(
Pif

nγ
if

)
= 0 (2.21)

The guiding center equations for the particle ions are essentially identical to those of the electrons

in Equations (2.1) and (2.2),

µip =
p2ip,⊥
2B

= const. (2.22)

and

dpip,∥
dt

= pip,∥vE · κ− µib · ∇B + eE∥ (2.23)

The ions are taken to be non-relativistic. Including the new species also requires changes in the

expressions for the electron fluid density

nef = nif + nip − nep (2.24)

and parallel electron fluid velocity

nefvef,∥ = nifvif,∥ + nipvip,∥ − nepvep,∥ (2.25)

40



The adiabatic equation of state for the electron fluid, Equation (2.14), remains the same, as do

equations (2.15), and (2.16) for B and E⊥. Equation (2.17) for E∥ remains the same because

only the electrons contribute.

In both the original and upgraded model, electron inertia can be discarded compared with

that of the ions in the ion momentum equation, as discussed in Sec. 2.1.1. The distinction between

the inertia of collective and individual particle electrons, along with their movements in relation to

the ambient magnetic field and their specific interactions with the MHD fluid, is comprehensively

discussed in [103]. However, as shown in equation (2.18), the inertia of the particle ions cannot

be neglected. The derivation of the new fluid ion momentum equation is detailed in Appendix

A.1, with the final result being

mi
∂ (nivif )

∂t
= −∇·Tif − (∇·Tef )⊥+J×B/c−∇·Tip− (∇·Tep)⊥+eniE∥b−Ii∥ (2.26)

where Ii∥ is the difference in the parallel inertia between the fluid and particle ions, given by

Ii∥ =
∂

∂t
minip(vip∥ − vif∥)b = −bmivif∥

(
nip

nif

∇ · nifvif∥b−∇ · nipvip∥b− nifvi,⊥ · ∇nip

nif

)
+

nip

nif

bb · ∇ · Tif − bb · ∇ · Tip +minip

vip∥ − vif∥
B

(I− bb) · ∂B
∂t

(2.27)

where ni = nip+nif , and the ion fluid stress tensor, Tif , is defined in Equation (2.5). The particle

ion pressure tensor is given by

Tip = Pip∥bb+ Pip⊥(I− bb) +minipvipvip (2.28)

with the parallel and perpendicular pressures calculated in the frame of the bulk motion of the

41



particle ions,

Pip⊥ =

∫
dpi

p2i⊥
2mi

f (2.29)

and

Pip∥ =

∫
dpi

p2i∥
mi

f (2.30)

with pip the ion momentum. Note that in the inertia term on the left side of Equation (2.26), ni

is the total ion density and not just the fluid ion density nif . In the limit nip → 0, Ii∥ → 0 and

Tip → 0 and equation (2.26) reduces to equation (2.4).

The ratios nip/nif and nep/nef are typically chosen at t = 0 to be identical throughout

the domain, but are generally free to evolve in space and time during simulations. The electron

spectra reported in previous kglobal reconnection simulations were insensitive to nep/nef [85],

but the sensitivity of reconnection simulations with both particle electrons and ions to these ratios

will also need to be explored.

To summarize, the kglobal model with particle electrons and ions consists of Eqs. (2.1)-

(2.2), (2.6)-(2.9), (2.14)-(2.17), and (2.20)-(2.30). Because of the complexity of Eqs. (2.26) and

(2.27), an alternative is to advance the equation for the center-of-mass velocity in Eq. (A.26) and

calculate the ion fluid velocity by subtracting the ion particle flux.

2.2 Alfvén Wave Propagation and the Growth of Firehose Modes in an Anisotropic

Plasma

Particle energy gain in reconnection is dominantly parallel to the ambient magnetic field,

with the implication that the parallel pressure exceeds the perpendicular pressure in the heated

42



plasma as reconnection develops. As a consequence, the cores of magnetic islands approach the

marginal firehose condition at late time in simulations. The tension in a bent magnetic field line,

which is the fundamental driver of reconnection, goes to zero at the marginal firehose condition.

Thus, the pressure anisotropy that develops as reconnection proceeds regulates reconnection and

particle energy gain [53, 106, 103, 88]. Therefore, to ensure that the updated kglobal model

properly reproduces the influence of pressure anisotropy on the magnetic field dynamics, we

benchmark the model with a circularly polarized Alfvén wave mode propagating along a uniform

magnetic field in a system with pressure anisotropy (P∥ ̸= P⊥), a problem with a well-established

solution [31]. With increased pressure anisotropy, the reduction in magnetic tension of a propa-

gating Alfvén wave reduces the wave phase speed. In a second benchmarking test on the dynam-

ics of the pressure anisotropy of both particle ions and electrons, we explore the growth of the

firehose instability.

The upgraded kglobal model builds upon the original kglobal framework, the details of

which can be found in the references [103, 88]. Implementation of the model uses a code written

in normalized units. A reference magnetic field strength B0 and the initial total ion density ni0 =

nip + nif define the Alfvén speed, CA =
√

B2
0/(4πmini0). Lengths and times are normalized

to an arbitrary macroscale length, L, and the Alfvén crossing time, τA = L/CA, respectively.

Electric fields and temperatures are normalized to CAB0/c and miC
2
A, respectively. The kglobal

model employs second-order diffusivities and fourth-order hyperviscosities (∝ ν∇4) in each

fluid equation to mitigate noise at the grid scale [103].

The simulations were performed within a two-dimensional (2D) spatial domain, coupled

with a three-dimensional velocity space, featuring an initially uniform magnetic field of magni-

tude B0 oriented along the x-axis. The domain measures 2πL × πL with 512 × 256 cells with

43



periodic boundary conditions on all sides. Each grid cell is initialized with 200 self-consistent

(i.e., not test) particles, 100 ions and 100 electrons. The particle-to-fluid density ratios are ini-

tially set to nip/nif = nep/nef = 1/3 everywhere, but evolve in time at each location. The

ion-to-electron mass ratio is set to 25, although the results of the simulations are insensitive to

this value. The scalar temperature of both fluid species is set to 0.1. For the particle ions and

electrons, the perpendicular temperature is initially set to 0.1 while the parallel temperature is

adjusted to control the magnitude of the anisotropy.

We introduce perturbations to the magnetic field and velocity in the form of a circularly

polarized Alfvén wave, with a wavelength equal to the length of the simulation box. During a

propagation interval of 3τA, we measured the wave’s velocity in order to compare it with the

theoretical phase speed Cp of an Alfvén wave, which is given by

Cp = CA

√
1−

4π(P∥ − P⊥)

B2
≡ αCA (2.31)

where α =
√
1− 4π(P∥ − P⊥)/B2, aligns with the Chew-Goldberger-Low (CGL) framework in

the linear regime where pressure perturbations are negligible [107]. Here, P∥ and P⊥ denote the

combined ion and electron pressures along and perpendicular to the magnetic field, respectively.

44



The results are shown in Fig 2.1, where Cp, normalized to CA is plotted against α. There is

excellent agreement between the code results and linear wave theory, highlighting the precision

of the kglobal model and suggesting that the impact of ion and electron pressure anisotropy on

reconnection dynamics will be properly included in the new model.

45



0.0 0.2 0.4 0.6 0.8 1.0
= 1 4 (P P )/B2

0.0

0.2

0.4

0.6

0.8

1.0

C p
/C

A

Regression Line: y = 0.98x + 0.02

Figure 2.1: The observed phase speed of the Alfvén wave, Cp, plotted against the anisotropy
parameter α. The blue dots are the results of simulations. The intersection of the vertical and
horizontal dotted lines marks the isotropic Alfvén wave.

46



In the second benchmark, the rate of growth of the firehose instability was explored with

a pre-set pressure anisotropy above the instability threshold (α2 = −0.2). The system was ini-

tialized with small amplitude periodic disturbances with seventeen distinct wavelengths, defined

by k = m/2πL, where m is the mode number. The simulations incorporated a kinematic vis-

cosity of ν = 5.0 × 10−5. The growth rate, which includes the pressure anisotropy drive and

viscous damping is given by γ = kCA|α| − νk4. The fourth-order hyperviscosity sets the up-

per bound of the instability at small spatial scales. Figure 2.2 presents the expected theoretical

growth rate (represented by a solid blue line) and the growth rates found in the simulations for

a range of unstable modes. The red crosses indicate systems in which the electrons and ions

contribute equally to the total anisotropy. To explore the effect of unequal pressure anisotropies,

we also performed several simulations where the ions contribute 3/4 and the electrons contribute

1/4 of the total anisotropy (green squares). The growth rates measured in the unequal anisotropy

simulations closely follow those from equal anisotropy simulations and the results from all the

simulations are consistent with the predictions of linear theory. Note that there are no spurious

unstable modes observed at short scales such as those that arise in some numerical models when

the mode frequency exceeds the cyclotron frequency [108]. All waves in the kglobal model are

below the cyclotron frequency.

47



2 4 6 8 10 12 14 16
m

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

A

Theoretical Growth Rate Line
Numerical Growth Rates with Equal Electron-Ion Pressure Anisotropy
Numerical Growth Rates with Unequal Electron-Ion Pressure Anisotropy

Figure 2.2: Normalized growth rate of the firehose mode, γτA, versus mode number, m = 2πkL,
over a span of unstable m values. Red crosses show values from simulations with equal electron
and ion pressure anisotropies; values for unequal anisotropies are shown as green squares. The
blue line is the theoretical growth rate.

48



Chapter 3: Simultaneous Proton and Electron Energization during Macroscale

Magnetic Reconnection

The material in this chapter has been adapted from [109] with permission from the authors.

3.1 Simulation Model Setup

In this chapter we present results from simulations of magentic reconnection in a macro-

scopic domain (i.e., one not constrained by kinetic scales). The simulations include four distinct

plasma species: fluid protons and fluid electrons (which collectively form the MHD backbone),

and particle protons and electrons which are are represented by macro-particles that move through

the MHD grid. The particles are treated in the guiding-center limit, thus eliminating the need to

resolve their respective Larmor radii. Although the computational domain does not allow spatial

variations in one direction, such a configuration is known to support reconnection and has been

used in many previous simulations. Particles move across the magnetic field with their E × B

drift and along the local magnetic field at their parallel velocity.

As with the previous version of kglobal [85], the reconnecting component of the upstream

magnetic field B0 (along the x direction) and the total proton density (the combined number den-

sity of particle and fluid protons) ni0 serve as normalization parameters by defining the Alfvén

speed CA0 = B0/
√
4πmini0. Because kinetic scales are excluded in the model, lengths are nor-

49



malized to an arbitrary macroscale L0 and time scales are normalized to τA = L0/CA0. Both tem-

peratures and particle energies are normalized to miC
2
A0. The perpendicular electric field in the

simulations follows the MHD scaling CA0B0/c, while the parallel field scales as miC
2
A0/(eL0).

Although the parallel electric field is small compared to the perpendicular component, the energy

associated with the parallel potential drop over the scale L0 is of the order of miC
2
A0, making it

comparable to the available magnetic energy per particle.

The proton-to-electron mass ratio is set to 25, but our results are insensitive to this pa-

rameter. The speed of light in our simulations is 30 times the Alfvén speed. All species (both

fluid and particle) begin with an isotropic temperature of 0.0625miC
2
A0, with the particles hav-

ing a Maxwellian velocity distribution. The simulations are initialized with constant densities

and pressures in a force-free current sheet with periodic boundary conditions. Thus, B =

B0 tanh(y/w)x̂ +
√

B2
0 sech

2(y/w) +B2
g ẑ where Bg is the asymptotic out-of-plane magnetic

field (the guide field) and w, the width of the current sheet, is set to 0.005L0. The initial total

density of electrons (ne) and protons (ni) is normalized to unity, with particles comprising 25% of

the density (nep for particle electrons and nip for particle protons) and the remaining 75% in the

fluid component (nef for fluid electrons and nif for fluid protons). We demonstrate that results

of the simulations are insensitive to this fraction. The simulations are conducted on grids with

varying resolutions, ranging from 1024× 512 to 8192× 4096 grid points as shown in Table 3.1.

The number of particles is set to 100 particles per grid. The time step δt is adjusted as the grid

size is varied to ensure that particle stepping is accurate (e.g., cδt ≃ δ with δ the grid scale). The

values of the time-step are given in Table 3.1.

50



In our simulations, diffusion and hyperviscosity terms are included to ensure numerical

stability while reducing any high-frequency noise arising at the grid scale. The diffusion coeffi-

cients, denoted as Dn for number density diffusion and Dp for pressure diffusion, are set to the

values listed in Table 3.1. A hyperviscosity ν is included in the magnetic field evolution equation

to facilitate reconnection while minimizing dissipation at large scales. It is applied as a fourth-

order Laplacian term (∇4). The same hyperviscosity is included in the evolution equations for the

fluid proton flux, fluid proton number density and fluid proton pressure. The effective Lundquist

number Sν = CAL
3
0/ν associated with the hyperviscosity is varied to change the effective system

size (the ratio of the macro to the dissipation scale). The hyperviscosity coefficients, denoted as

νB, νnv, νn, and νp, as well as Sν , are set to the values given in Table 3.1.

3.2 Simulation Results

In our simulations the hyperviscosity triggers reconnection, leading to the formation of

small-scale flux ropes. These flux ropes coalesce over time and eventually form system-scale

structures. Figure 3.1 shows the results of a simulation with Bg/B0 = 0.25. Shown is the density

of particle protons in the x–y plane at five times, t/τA = 0, 4, 7, 10, and 13, in panels (a)–(e).

51



Table 3.1: Parameters for Simulation Domains

Grid Points 1024× 512 2048× 1024 4096× 2048 8192× 4096

Time Step (in τA) 2× 10−4 1× 10−4 5× 10−5 2.5× 10−5

Proton Number Density Diffusion (Dn) 10.5× 10−5 5.25× 10−5 2.625× 10−5 1.3125× 10−5

Proton Pressure Diffusion (Dp) 10.5× 10−4 5.25× 10−4 2.625× 10−4 1.3125× 10−4

Hyperviscosity
(νB, νnv, νn, and νp) 84× 10−9 10.5× 10−9 1.312× 10−9 0.1640× 10−9

Effective Lundquist Number (Sν) 1.2× 107 9.5× 107 7.6× 108 6.1× 109

52



At t = 0 the magnetic field reverses across a uniform current sheet as shown in panel

(a) of Figure 3.1. As the system evolves, reconnection begins at multiple sites. The resulting

reconnected magnetic field lines convect away from the x-points. The contraction of magnetic

islands along the current sheet leads to the energization of particle electrons and protons trapped

within them. As the simulation evolves, these islands grow and undergo coalescence, culminating

in the formation of a single, large magnetic island late in time. As can be seen in Figure 3.1

the particle proton density remains relatively uniform along the field lines due to the protons’

mobility parallel to the magnetic field.

53



0.5
0

-0.5

(a)

0.5
0

-0.5

(b)

0.5
0

-0.5

(c)

0.5
0

-0.5

(d)

-3 -2 -1 0 1 2 3
x/L0

0.5
0

-0.5

(e) 0.2

0.4

0.6

0.8

1.0

pa
rti

cl
e 

pr
ot

on
 d

en
si

ty

y/L0

Figure 3.1: The temporal evolution of the particle proton density in the x–y plane for Bg/B0 =
0.25. Five times are displayed: t/τA = 0, 4, 7, 10, and 13 in panels (a) to (e). In each panel,
magnetic field lines are shown in black.

54



Figure 3.2 shows the energy per particle (energy density divided by number density) of the

particle protons and particle electrons, ⟨W ⟩, and the firehose parameter, α = 1−4π(P∥−P⊥)/B
2,

at late time, t/τA = 13, from the same simulation. In the top two panels darker colors represent

higher average energies per particle. Thus, the average energy of the protons is significantly

greater than that of the electrons.

The firehose parameter in the bottom panel of Figure 3.2 is a measure of the impact of

pressure anisotropy on the local magnetic tension, given by αB · ∇B/4π. Specifically, when α

approaches zero the magnetic tension, which is the driver of magnetic reconnection also goes to

zero. Thus, the dominant feedback mechanism of energetic particles on reconnection dynamics is

through the reduction of magnetic tension as measured by the firehose parameter. In the context

of the simulation from kglobal, both P∥ and P⊥ are derived from the particle species since both

the fluid proton and electrons are taken to be isotropic and do not impact the magnetic tension

force. Thus, in Figure 3.2 extensive areas within the magnetic islands approach marginal stability

(α = 0), punctuated by small unstable regions (red shades). Consequently, the local magnetic

tension, a crucial factor in driving reconnection and particle energy gain, is strongly reduced

within the large flux rope and somewhat reduced in the smaller flux ropes. This highlights the

critical importance of particle feedback on the MHD fluid dynamics. Traditional MHD models

that include test-particle dynamics do not include this feedback and thereby risk an uncontrolled

increase of the particle energy.

55



0.5

0

-0.5 t = 13 A

0.5

0

-0.5 t = 13 A

-3 -2 -1 0 1 2 3
x/L0

0.5

0

-0.5 t = 13 A

0.0

0.2

0.4

0.6

0.8

1.0

<
W

>
/m

iC
2 A0

1.0

0.5

0.0

0.5

1.0

y/L0

Figure 3.2: The energy per particle of the particle protons (top) and particle electrons (middle),
⟨W ⟩, is shown in the top two panels, respectively, for the simulation with Bg/B0 = 0.25. The
darker the color, the higher the average energy of the particles so the protons contain significantly
more energy than the electrons. The firehose parameter α at t = 13τA is shown in the bottom
panel. The color white corresponds to marginal stability and reveals that within flux ropes where
particles have been energized, the particle feedback has significant impact on the MHD environ-
ment.

56



At late time, after a single large magnetic island dominates the current sheet, we calculate

the proton energy spectra by aggregating particle counts throughout the entire simulation domain.

Figure 3.3(a) displays the differential number densities F (W ) = dN(W )/dW plotted against

normalized energy W/miC
2
A0 at various times for Bg/B0 = 0.25. F (W ) begins as a Maxwellian

but develops a power-law tail, characterized by an index δ′, that first appears for low energies

before eventually extending to near the maximum energy observed in the domain. The powerlaw

index at late time approaches 3.4. Similar behavior was observed for electrons in the earlier

electron-only version of kglobal [85]. Note that the particle distribution at low energy in Figure

3.3(a) does not change significantly in time. This is not because the protons are not heated during

reconnection but because the distribution in the figure is from the entire computational domain,

which includes large regions of plasma upstream of the reconnecting current layer that has not

been accelerated. The heating of the plasma will be more evident when we focus on the particle

distributions within the large islands that develop at late time in the simulations.

In Figure 3.3(b) we show the late-time F (W ) spectra for protons for different Bg values,

measured at times when approximately equal amounts of magnetic flux have undergone recon-

nection. A similar figure from the electron-only version of kglobal was presented in [85]. In both

cases as Bg decreases, the power-law index δ′ also decreases. This results in a harder energy

spectrum and leads to the generation of a larger number of high-energy particles for small values

of Bg.

57



10
2

10
1

10
0

10
1

10
2

10
3

W/miC
2

A0

10
10

10
8

10
6

10
4

10
2

10
0

10
2

F(
W

) (
A

rb
. u

ni
ts

)

~W
3.4

Proton
(a)

t/ A = 0

t/ A = 1

t/ A = 2

t/ A = 3

t/ A = 4

t/ A = 5

t/ A = 6

t/ A = 7

t/ A = 8

t/ A = 9

t/ A = 10

t/ A = 11

t/ A = 12

t/ A = 13

10
2

10
1

10
0

10
1

10
2

10
3

W/miC
2

A0

10
10

10
8

10
6

10
4

10
2

10
0

10
2

F(
W

) (
A

rb
. u

ni
ts

)

Proton
(b)

Bg/B0 = 0.25

Bg/B0 = 0.5

Bg/B0 = 0.75

Figure 3.3: A log-log plot of the differential number density versus the normalized energy. Panel
a shows the temporal evolution of the proton distribution in a simulation with Bg/B0 = 0.25.
The initial Maxwellian spectrum develops an energetic tail that hardens over time and eventually
stabilizes into a power law. Panel b shows the late time spectra for different guide fields. Smaller
guide fields result in harder spectra that extend to higher energies.

58



In addition, we explore the late-time distributions of energetic protons for a guide field

Bg/B0 = 0.25 and various values of Sν (effective system size). The results in Fig. 3.4 show

that the power-law indices of the spectra are only weakly dependent on the system size for large

domains. As the size of the system increases, the spectra become slightly harder and the upper

limit in the particle energy gain increases. This is because in larger domains more flux ropes form

in the reconnecting current sheet and it is the merger of those flux ropes that efficiently drives

particle energy gain.

59



10−2 10−1 100 101 102 103 104

W/miC2
A0

10−12

10−10

10−8

10−6

10−4

10−2

100

102

F(
W

) (
A

rb
. u

ni
ts

)

Proton Sν=1.2×107
Sν=9.5×107
Sν=7.6×108
Sν=6.1×109

Figure 3.4: The late-time distributions of energetic protons in a simulation with Bg/B0 = 0.25
for various values of Sν (effective system size) reveal that the power-law indices of the spectra
vary weakly with the domain size. The rightmost point of each curve marks the one-count level,
indicating the flux when there is just one particle count per energy bin. The notable difference is
that as the size of the system increases, the maximum energy of the protons also increases.

60



A free parameter in the kglobal model is the fraction of the particle electrons and protons

compared with their respective fluid component [103, 88, 105]. An important question is whether

reconnection simulations are sensitive to these ratios. To address this question, we investigated

the spectra of protons and electrons with a guide field Bg/B0 = 0.25 and with different fractions

of the particle densities. Figure 3.5 shows the late-time energy spectra of particle protons and

electrons from three simulations with differing particle fractions. The power-law indices of the

spectra are nearly identical, indicating that the simulation results are insensitive to this ratio.

The only notable difference is that as the fraction of particles decreases, the maximum achieved

energy of the particles increases modestly. This can be attributed to the fixed total energy in the

simulation box. The particles gain more energy than their fluid counterparts with the parallel

energy gain exceeding that of the perpendicular energy. Thus, when there are fewer particles in

the simulation box, the particles can gain more energy before the firehose parameter acts to limit

reconnection and energy gain.

61



10−2 10−1 100 101 102 103

W/miC2
A0

10−8

10−6

10−4

10−2

100

102

F(
W

) (
A

rb
. u

ni
ts

)

nip/ni = 40%
nip/ni = 25%
nip/ni = 10%

10−2 10−1 100 101 102 103

W/miC2
A0

10−8

10−6

10−4

10−2

100

102

F(
W

) (
A

rb
. u

ni
ts

)

nep/ne = 40%
nep/ne = 25%
nep/ne = 10%

Figure 3.5: The simultaneous distributions of protons (top) and electrons (bottom) from simula-
tions with a guide field Bg/B0 = 0.25 and different fractions of particle densities. The power-law
indices of the spectra for the different ratios are nearly the same, indicating that the simulation is
insensitive to this ratio. The only notable difference is that as the fraction of particles decreases,
the maximum achieved energy of each species increases modestly.

62



We show in Figure 3.6 the initial and the late-time spectra of particle protons and electrons

for the case Bg/B0 = 0.25. The initial curves for both protons and electrons are Maxwellian

distributions. Late-time spectra exhibit power-law tails for both particle types. The spectral

indices of the powerlaws of the two species are nearly the same. Although the length of the

power-law tail for the protons is shorter than that of the electrons, it extends nearly a full decade

higher in energy. A key feature is the significantly greater energy gain of the protons compared

to the electrons, which is consistent with the increased energy content of protons compared to

electrons shown in Figure 3.2.

63



10
2

10
1

10
0

10
1

10
2

10
3

W/miC2
A0

10
11

10
9

10
7

10
5

10
3

10
1

10
1

F(
W

) (
A

rb
. u

ni
ts

)

Electron
Proton

Figure 3.6: The simultaneous distributions of protons (blue curves) and electrons (orange curves)
for a guide field Bg/B0 = 0.25. The curves on the left indicate the initial Maxwellian distribu-
tions for both particle types. The hard spectra are the late time results, showing extended power
law distributions for both protons and electrons with powerlaw indices that are comparable.

64



We next explore the details of particle heating and energization by examining the rela-

tive numbers and energy content of the non-thermal and hot thermal particles. The spectra are

modeled with the sum of a kappa [110] and a Maxwellian distribution (See Appendix B.1 for

a detailed explanation of kappa distribution. The kappa distribution turns into a powerlaw at

high energy which can represent the extended powerlaw tails seen in the simulation data [8].

This combination is particularly useful in modeling the distributions of particles in space and

astrophysical plasmas where non-thermal behaviors are observed [111, 7, 8].

In Figure 3.7 we show F (W ) (black) for protons and electrons at late time after the dis-

tribution functions have reached a steady state. The data is taken from within the large late-time

island shown in Figure 3.1(e). Thus, the spectra shown include only particles that have undergone

heating and energization during reconnection. Shown in Figure 3.8 are the intial and final elec-

tron and proton distributions from the same large island but using a linear scale. This figure more

clearly shows the changes in the distribution of low energy particles as they are heated during

reconnection. The spectra in Figure 3.7 are fit (solid red) by the sum of a Maxwellian distribution

and a kappa distribution [85]. This procedure is designed to capture the characteristics of both

hot thermal and nonthermal particle populations. From this fitting approach, we extract several

parameters:

1. The spectral index δ′ of the powerlaw of the non-thermal particles.

2. The relative densities and total energy of the hot thermal and non-thermal electrons and

protons.

3. The break point energy (vertical line in Figure 3.7) , expressed in units of miC
2
A0, that

marks the transition from thermal to non-thermal particle dominance in the spectra. At the

65



10 2 10 1 100 101 102 103

W/miC2
A0

10 10

10 8

10 6

10 4

10 2

100

102

F(
W

) (
Ar

b.
 u

ni
te

s)

~W 3.4

x = 1.34

Proton
Fitting
Kappa Thermal
Kappa Nonthermal
Maxwellian

10 2 10 1 100 101 102 103

W/miC2
A0

10 10

10 8

10 6

10 4

10 2

100

102

F(
W

) (
Ar

b.
 u

ni
ts

)

~W 3.3

x = 0.28

Electron
Fitting
Kappa Thermal
Kappa Nonthermal
Maxwellian

Figure 3.7: Fitting and decomposition of the late time spectra of electrons and protons. The solid
black curves are the proton (top) and electron (bottom) distribution functions for a guide field of
0.25. The solid red curves show the fit to a sum of a Maxwellian and a kappa distribution. The
green dashed curve is the thermal part of the kappa distribution, while the blue dashed curve is the
non-thermal part of the kappa distribution. The yellow curves are the Maxwellian distributions
of the two species. The pink vertical line indicates the energy where the nonthermal and thermal
components of the kappa distribution are equal. The solid blue line is the asymptote for the
distribution function in the power law region.

66



break point the number densities of the hot thermal and non-thermal components of the

kappa distribution are equal (See Appendix B.1 for a detailed explanation of the procedure

used to separate the hot thermal and non-thermal components of the kappa distribution).

To the left of the break point, the thermal component dominates, while to the right, the

thermal components decrease sharply, and the non-thermal component dominates.

67



Table 3.2: Density Distributions for Particle Acceleration at Late Time

Type Guide Field Maxwellian Kappa Total Kappa Maxwellian Kappa Nonthermal
(in B0) (%) (%) (%) (%)

Electron 0.25 38.2 61.8 45.1 16.7

Electron 0.50 68.0 32.0 26.8 5.2

Electron 0.75 84.3 13.7 12.4 1.3

Proton 0.25 51.7 48.3 35.8 12.5

Proton 0.50 43.1 56.9 46.7 10.2

Proton 0.75 54.5 45.5 40.3 5.2

68



0.0 0.2 0.4 0.6 0.8 1.0
W/miC2

A0

0

2

4

6

8

10

F(
W

) (
A

rb
. u

ni
ts

)

Initial Electron
Initial Proton
Late-Time Electron
Late-Time Proton

Figure 3.8: The comparison between the initial spectra and late-time spectra of protons and
electrons is presented on a linear-linear scale to more clearly reveal the changes in the particle
distributions at low energy.

69



Table 3.2 documents how the strength of the guide field affects the proportion of particles

that comprise the Maxwellian and kappa distributions, with the latter further subdivided into

thermal and non-thermal components. As the guide field increases, the proportion of non-thermal

particles, indicated by the kappa non-thermal percentage, consistently decreases. For electrons,

the non-thermal component reduces significantly from a higher percentage at a lower guide field

to a much smaller fraction at the higher guide field. A similar trend is observed for protons, where

the non-thermal proportion diminishes as the guide field strength increases. This trend establishes

that a strong guide fields suppresses the generation of non-thermal electrons and protons, which

is consistent with the findings for energetic electrons in [85].

Table 3.3 shows the relative energy of thermal and non-thermal electrons and protons at

late time, again as a function of guide field strength. In this data the total number density of

particles is normalized to unity and the numbers in the table reflect the energy content of the

various components in units of miC
2
A. As the guide field increases, the energy in the non-thermal

component of the kappa distribution decreases for both electrons and protons. For electrons,

the kappa non-thermal energy decreases from 0.095 at a guide field of 0.25 to 0.006 at a guide

field of 0.75. For protons, the kappa non-thermal energy decreases from 0.340 at a guide field

of 0.25 to 0.045 at a guide field of 0.75. The downstream thermal temperatures are also shown.

For electrons the thermal temperature, represented as a weighted sum of the Maxwellian and

kappa-Maxwellian distributions, remains comparable for different guide fields. The temperature

increase is modest compared with the initial value of 0.0625miC
2
A0. The heating is driven by the

macroscale parallel electric field as explored by the reference [85]. In contrast, protons exhibit

a much larger increase in temperature as has been documented in the Earth space environment

[112, 113, 114]. For the lowest guide field in Table 3.3 the ratio of the temperature increment of

70



the protons to electrons is around five, which is slightly smaller than in measurements in the Earth

space environment [113]. Also noticeable is the decrease in thermal temperature of protons with

increasing guide field. We emphasize that because the data in Table 3.3 is normalized to miC
2
A0,

the temperature increments of both species scales with this parameter as in the observations

[115, 112, 113, 114].

Also shown in Table 3.3 is the breakpoint energy, which marks the transition from t