# Automated title and abstract screening for scoping reviews using the GPT-4
Large Language Model

David Wilkins  Discipline of General Practice
The University of Adelaide

###### Abstract

Scoping reviews, a type of systematic literature review, require intensive
human effort to screen large numbers of scholarly sources for their relevance
to the review objectives. This manuscript introduces GPTscreenR, a package for
the R statistical programming language that uses the GPT-4 Large Language
Model (LLM) to automatically screen sources. The package makes use of the
chain-of-thought technique with the goal of maximising performance on complex
screening tasks. In validation against consensus human reviewer decisions,
GPTscreenR performed similarly to an alternative zero-shot technique, with a
sensitivity of 71%, specificity of 89%, and overall accuracy of 84%. Neither
method achieved perfect accuracy nor human levels of intraobserver agreement.
GPTscreenR demonstrates the potential for LLMs to support scholarly work and
provides a user-friendly software framework that can be integrated into
existing review pipelines.

##  1 Introduction

A scoping review is a relatively novel type of literature review that aims to
map the key concepts and existing activity within an area of research [1].
Like systematic reviews, scoping reviews typically use rigorous, transparent
[13], and sometimes pre-registered methods for gathering and synthesising
evidence, and increasingly use formal frameworks for both performing and
reporting reviews [11]. Scoping reviews can inform future systematic reviews
or primary research in the same area [18]. However, they differ from
systematic reviews in aiming to describe the breadth of coverage of the
available literature rather than research findings in depth [1].

Frameworks for performing a scoping review typically involve defining a
research area or question, searching bibliographic databases for potentially
relevant published material (‘sources’), screening these sources to identify
those relevant to the area or question, and systematically extracting and
reporting data from the relevant sources [20, 1]. The screening stage will
usually involve initial screening of source titles and abstracts against pre-
determined inclusion and exclusion criteria, followed by full text screening,
with both steps performed in replicate by at least two human reviewers [10,
12, 13]. Because database searches can return many hundreds or thousands of
potentially relevant sources, these screening steps can require intensive
human time and effort. Many software methods have been proposed or used to
support or partially automate source screening for scoping reviews, including
text mining to prioritise relevant sources for human screening [16, 5, 2],
automated clustering and labelling of sources to support human decision-making
[17], and ‘crowdsourcing’ screening to untrained workers via online platforms
[8]. A similar but more extensive set of methods have been developed and
employed for systematic reviews [6, 3] for which the process of source
screening is broadly comparable.

Since the release of the first Generative Pre-trained Transformer (GPT) Large
Language Model (LLM) by OpenAI (San Francisco, California, United States of
America) in 2018 [15], transformer-based LLMs and the GPT lineage in
particular have seen rapid and widespread adoption for a range of automation
tasks. Broadly, these models generate a probabilistically weighted list of
‘tokens’ (parts of text such as letter combinations and punctuation) to
continue or complete some input text (a ‘prompt’), having been trained to do
so by practising such predictions on large human-written corpora. When this
generative process is iterated, it allows for a range of applications
involving analysis and production of text, such as summarising documents,
generating fiction in a particular genre or style, or conversing with humans
[9].

While LLMs are not yet widely used to screen sources for literature reviews,
early work suggests they may perform well in this role. Guo et al. [4]
reported the use of a GPT-lineage model (they do not specify which, though
their published code suggests OpenAI’s ‘gpt-3.5-turbo’ model) to screen 24,307
titles and abstracts from five systematic reviews and one scoping review,
achieving weighted average sensitivity of 76% and specificity of 91% when
compared to human reviewers. Their approach involves giving the model a brief
prompt instructing it to take on the persona of a researcher screening titles
and abstracts, followed by a source’s title and abstract as well as the
inclusion and exclusion criteria. The model is instructed to respond with a
decision to include or exclude the source, and the process is iterated across
the full set of sources to be screened. Syriani et al. [19] similarly reported
the use of ‘gpt-3.5-turbo’ to screen titles and abstracts for a systematic
review and achieved sensitivities of above 70%. They also systematically
evaluated prompts given to the LLM to identify a prompt that performed best at
the screening task; their chosen prompt, like that of Guo et al., placed the
LLM in the role of an academic reviewer.

Both of these approaches made use of a single, fixed text prompt template,
which the LLM then completes with additional text representing its response
(the decision to include or exclude a source), a method sometimes called
‘zero-shot prompting’. Recent work has identified a number of methods which
can be superior to zero-shot prompting when using LLMs for tasks that require
complex or multi-step reasoning. These methods include ‘chain-of-thought
prompting’ [21], in which a complex task is broken down into a series of
intermediate steps, and the ‘tree of thoughts’ strategy [23], in which
multiple parallel chains of thought are generated, compared, and integrated.
The LLM is induced to follow these complex reasoning strategies either by
being given examples of multi-step reasoning on similar tasks, or by being
lead through the process with a series of intermediate prompts.

In this paper, I introduce a package for the R programming language [14]
called GPTscreenR that implements a chain-of-thought approach to using GPT-4
for scoping review title and abstract screening, and evaluate its performance
by comparison to human reviewers. The purpose of this package is to assist and
augment rather than replace human reviewers in performing scoping reviews.
This package represents the first LLM-based screening tool designed
specifically for scoping reviews. Further, this study provides the first
report on the accuracy of LLM-based screening using the most recent iteration
of the GPT model lineage, GPT-4, and using the recently developed chain-of-
thought approach.

##  2 Methods

###  2.1 The GPTscreenR package

GPTscreenR is an R [14] package released under the MIT open source licence.
The source code is available for download from GitHub at
https://github.com/wilkox/GPTscreenR. At the time of writing the most recent
package version was 0.0.3, and the results presented in this paper were
obtained using this version.

GPTscreenR consists of two main components. The first is a set of internal
functions for interfacing with the OpenAI Application Programming Interface
(API), which allows for ‘conversations’ with the GPT-4 LLM, as well as
functions for representing and manipulating those conversations. These
internal functions are designed to be model-agnostic, so that future versions
of the package or users with particular needs can use different GPT models.
The OpenAI API requires an OpenAI account, and OpenAI charges fees for use of
the API. In order to access the API, GPTscreenR requires a secret key to be
registered prior to source screening, and instructions for doing so are
provided in the package documentation and on loading of the package in R if
the key has not been correctly registered.

The second component is a set of user-facing functions to perform source
screening with GPT-4. The review_description() function assists in generating
a text description of the review’s objectives and inclusion and exclusion
criteria, using the Population, Concept, and Context (PCC) framework [10, 12]
for defining the review’s inclusion and exclusion criteria. The use of this
function is optional, and users may instead choose to provide a description of
the review and criteria for source selection using any framework or format
they see fit.

The screen_source() function performs the main task of the package. This
function mediates an conversation with GPT-4 in which chain-of-thought
prompting [21] is used to guide GPT-4 through screening a source title and
abstract against the study inclusion criteria. The template for this
conversation is given in Figure 1. The OpenAI API defines a conversation as a
series of messages, each of which originates from one of three roles: system,
representing an authoritative voice that can instruct GPT-4 on its task and
behaviour; user, representing a human user that can interact with GPT-4; and
assistant, representing the responses generated by GPT-4. In screen_source(),
the system role gives GPT-4 general instructions, while the user role provides
the user-written review description and the source title and abstract.

Figure 1: Template for the conversation with GPT-4 mediated by the
screen_source() function. Messages with variable content, including user-
provided data as well as GPT-4’s responses, are given in italics.

The specific phrase ‘work step by step’ is derived from the ‘let’s think step
by step’ prompt phrase which significantly improves LLM performance on multi-
step reasoning tasks with a zero-shot prompt [7], adapted to this chain-of-
thought approach.

GPT-4 is instructed to summarise the inclusion criteria for the scoping
review, compare the title and abstract against these summarised criteria, and
make a final recommendation on whether to include or exclude the source
(Figure 1). This approach was chosen after noting that a major source of error
when attempting to screen sources with zero-shot, one-shot, or few-shot
prompts (i.e. with a single prompt and no, one, or a few examples) was that
GPT-4 would fail to consider important inclusion criteria. As an example,
Figure 2a presents a conversation with GPT-4 using a zero-shot prompt. The
screening task in this example is intentionally adversarial and designed to
lead the model towards making an error. In order to correctly recommend
exclusion of the source, GPT-4 must notice that the review is looking for
research on therapy alpaca interventions, but that the source reports on a
therapy camel intervention. The presence of multiple other inclusion criteria
which are met by the source, as well as the mention of alpacas in the source
abstract, serve as distractors. In this example, GPT-4 incorrectly recommends
inclusion. If the conversation is then continued to draw GPT-4’s attention to
the error, it is able to identify and correct it (Figure 2b), suggesting that
the error arises from a failure of GPT-4 to properly consider the relevant
inclusion criterion rather than an inability to do so. Using the chain-of-
thought approach overcomes this problem without the need for human
intervention and correction (Figure 2c). GPT-4 identifies ‘The context of the
study must involve the use of therapy alpacas as a part of a programme of care
in residential aged care facilities.’ as an inclusion criterion, and correctly
assesses that this is the only inclusion criterion not met. GPT-4 then
correctly recommends exclusion of the source.

Figure 2a: A conversation with GPT-4, using a zero-shot prompt to instruct the
model to screen a single source for a hypothetical scoping review. This
adversarial task has been intentionally designed to influence GPT-4 towards
making an error. Figure 2b: A continuation of the conversation in Figure 2a,
in which GPT-4 is lead to recognise and correct the error. Figure 2c: A
conversation with GPT-4, instructing it to perform the same adversarial task
presented in Figure 2a, but using the chain-of-thought approach. This is not a
continuation of the conversation presented in Figures 2a and 2b but a new
conversation. In this case, GPT-4 correctly recommends excluding the source.
The model also correctly identifies that the source meets all of the relevant
inclusion criteria except for the requirement that the source report on an
alpaca intervention.

The screen_source() function returns a list comprising the complete transcript
of the conversation with GPT-4 and GPT-4’s final recommendation as a logical
(Boolean) value. The conversation transcript can be used to interrogate cases
where GPT-4 may have returned an incorrect or unexpected result. The package
also provides a function screen_sources(), which applies screen_source()
iteratively across a data frame of sources. screen_sources() caches screening
outcomes to a file as they accumulate, so that screening can be split across
multiple sessions and recover from interruptions.

###  2.2 Validation

To validate GPTscreenR’s approach, six scoping reviews were identified from
the Open Science Framework (OSF; https://osf.io) where the review inclusion
criteria and the results of title and abstract screening were publicly
available. A summary of the review characteristics is provided in Table 1.
Random subsets were used where the large number of sources made it prohibitive
in time and cost to screen the full set. The total number of sources available
for screening and the number used for validation from each review are given in
Table 1.

Table 1: Scoping reviews used to validate the performance of GPTscreenR. ‘total available’ is the number of sources available for screening following attempted retrieval of missing abstracts, correction of malformed data, and removal of sources with missing, irreparable, or duplicate data; it may differ from the number reported by the review authors. Short name |  OSF project title | OSF URL |  File(s) |  Sources screened (total available)
---|---|---|---|---
‘COVID’ |  Complementary, Alternative, and Integrative Medicine-Specific COVID-19 Misinformation on Social Media: A Scoping Review | https://osf.io/ytz5e |  ‘TITLE ABSTRACT FULL TEXT SCREENING DATA_CAIM COVID-19 SM Misinfo_Mar2223.xlsx‘ |  324 (458)
‘melanoma’ |  Scoping Review – Melanoma | https://osf.io/knje4 |  ‘Melanoma Data.xlsx‘ |  72 (128)
‘smartphones’ |  Smartphone Addiction Scoping Review | https://osf.io/f9huw |  ‘Scop_rev_full.csv‘ |  256 (5,376)
‘solastalgia’ |  Australian Solastalgia Scoping Review | https://osf.io/qxe3n |  ‘Solastalgia re-screening and extraction spreadsheet Dec 2022.xlsx‘ |  150 (150)
‘suicide’ |  Social Norms and Suicidality - Scoping Review | https://osf.io/btpzc |  ‘June 2022 Screenings Extraction TOP UP.xlsx‘, ‘Nov20 Stage 2 Removal of duplicates and screenings.xlsx‘ |  100 (2,094)
‘teachers’ |  Teachers’ soft skills: a Scoping Review | https://osf.io/n9rkd |  ‘Teachers´ Soft Skills a SR - Screening NM & OJ(1).xlsx‘ |  245 (355)

Some of the reviews did not include source full abstracts in the publicly
available files, and when these abstracts could not be obtained from external
databases these sources were excluded from validation. There were also many
cases where missing, malformed, or duplicate data required either manual
correction or exclusion of sources. The scoping review data, code used to
prepare this data for validation, and code for calculating summary statistics
are available in a reproducible form in the package repository on GitHub
(https://github.com/wilkox/GPTscreenR/tree/master/validation).

The consensus human reviewer decision at the title and abstract screening
level was used as the gold standard outcome. Accuracy, sensitivity and
specificity were calculated by comparing GPT-4’s recommendation against this
gold standard. Three of the scoping reviews (COVID, solastalgia, and teachers)
included individual human reviewer decisions in addition to the final decision
in their publicly available datasets, and these were used to calculate human
intraobserver agreement (Cohen’s kappa) using the R function cohen.kappa()
from the psych package [22]. This was compared to human/GPT-4 agreement across
all screened sources, calculated with the same method.

###  2.3 Comparison to zero-shot method

In order to directly compare the chain-of-thoughts approach to a zero-shot
approach (i.e. a conversation consisting of a single prompt with no examples,
followed by GPT-4’s response), the validation screening task was repeated
using the prompt designed by Guo et al. [4], substituting the permitted
responses INCLUDE and EXCLUDE for included and excluded respectively in order
to maintain compatibility with GPTscreenR’s parsing of the response (Figure
3). In cases where a random subset of sources had been used for validation of
GPTscreenR, the same subset was used. The code used to prepare this data and
calculate summary statistics was otherwise identical to that used for the
chain-of-thoughts method validation and is available in the ‘zeroshot’ branch
of the package repository on GitHub
(https://github.com/wilkox/GPTscreenR/tree/zeroshot).

Figure 3: Zero-shot prompt used as a comparator for the chain-of-thoughts
approach, derived from the approach of Guo et al. [4]. Variable content,
including user-provided data as well as GPT-4’s response, is given in italics.

##  3 Results

1,147 sources were screened from the six scoping reviews. GPTscreenR achieved
an overall accuracy of 84% compared to the gold standard of the consensus
human reviewer decision, with weighted average sensitivity of 71% and weighted
average specificity of 89% (Figure 4). For the three reviews that provided
individual reviewer decisions, the weighted average Cohen’s kappa was 0.67,
while the weighted average Cohen’s kappa between final human and GPT-4
decisions was 0.52.

Figure 4: Accuracy, sensitivity, and specificity of GPTscreenR compared to the
gold standard of the consensus human reviewer decision. The box-and-whisker
plot represents for each statistic the median (central bar), 25th and 75th
percentiles (upper and lower hinges of the box), and largest value $\pm$ 1.5
$\times$ the interquartile range (whiskers). Points given the values for the
individual scoping reviews used for validation.

The zero-shot method achieved an overall accuracy of 83%, with weighted
average sensitivity of 72% and weighted average specificity of 87%. The
weighted average Cohen’s kappa between human and GPT-4 decisions was 0.52.

##  4 Discussion

###  4.1 Performance of GPTscreenR

GPTscreenR and its chain-of-thought approach performed comparably to the zero-
shot method reported by Guo et al. [4]. Also similarly to their report,
Cohen’s kappa was higher for inter-human than human-GPT-4 agreement,
suggesting that despite this generally good performance GPT-based methods
still do not match the performance of a human reviewer when the consensus
human decision is used as the gold standard.

###  4.2 Comparison to the zero-shot method

The replication of Guo et al.’s zero-shot prompt on the reviews used for
validation in this study resulted in similar weighted average sensitivity
(78%) compared to their report (76%), though somewhat poorer specificity (84%
compared to 91%). There are a number of possible factors contributing to this
difference. Firstly, this study used the GPT-4 model while Guo et al. likely
used the ‘gpt-3.5-turbo’ model, although it might be expected that the more
advanced GPT-4 would generally perform better on the same task. Secondly, this
study included only scoping reviews, while Guo et al. examined five systematic
reviews and only one scoping review, although the reported sensitivity (100%)
and specificity (94%) for that scoping review were higher than the weighted
average for this study. Thirdly, validation with both the chain-of-thoughts
and zero-shot methods was dependent on the availability of source and
screening data from the included reviews, and it was noted that most of the
reviews had issues with missing, malformed, missing, and duplicate data,
including lacking the text of some or all source abstracts and sufficient
information to retrieve them from public databases. While effort was made to
rectify these issues, many sources from the included reviews could not be used
for validation, and the differing quality of the datasets used in the two
studies may have affected the performance of the zero-shot method. Finally,
the different reviews used for validation in the two studies likely posed
different levels of difficulty the zero-shot approach. If LLM-based approaches
to source screening become more widely adopted, the performance of different
methods can be more accurately determined across a growing sample of reviews.

Compared against the zero-shot method, the chain-of-thought method achieved
poorer sensitivity but higher specificity. This result is consistent with
rationale for selecting the chain-of-thought method, which was to reduce type
I errors (false positives) of the type demonstrated in Figure 2a. However,
this came at the cost of higher type II error. It is not possible to
retrospectively examine the reasoning process that lead GPT-4 to make an
incorrect recommendation such as a false negative, as asking it to explain its
reasoning post-hoc will result in a confabulated statistically likely
response; as they have no direct access to their internal processes, LLMs are
unable to meaningfully introspect. However, unlike the zero-shot method, the
chain-of-thought method as implemented by GPTscreenR does allow reviewers to
view and assess the reasons for GPT-4’s final decision, by reviewing the
transcript of the conversation with GPT-4 and in particular its assessment of
the source against its summarised list of criteria. This may be useful for
resolving disagreements between human reviewers and GPT-4, or disagreements
among human reviewers.

###  4.3 Limitations of this approach

Compared to the prompts used by Guo et al. [4] and Syriani et al. [19], both
the prompt used for this approach and the response generated by GPT-4 have
substantially higher token counts. This results in both a longer time to
screen a single source (typically 20–30 seconds) and a higher dollar cost, as
the OpenAI API currently charges on a per-token basis. This may make
GPTscreenR less appealing to some reviewers, particularly when screening a
large number of sources.

Validation of GPTscreenR was limited by a small number of scoping reviews with
publicly available, high-quality data on human title and abstract screening
decisions. Because a subset of these reviews were used for both testing and
validation, there is a risk of over-fitting of the prompt to these particular
reviews. The public release of GPTscreenR encourages users to contribute data
from their own scoping reviews to support more accurately measuring
GPTscreenR’s real-world performance, and to guide further refinement of the
approach.

##  5 Conclusions

The use of Large Language Models to screen sources for scoping reviews is a
promising technology to reduce the human time and effort required to perform
these reviews. The GPTscreenR package performs comparably to a zero-shot based
approach to source screening. However, as with the zero-shot approach,
agreement between its decisions and the human reviewer consensus fell short of
inter-human agreement, suggesting that LLM-based screening is not yet as
reliable as a human reviewer. Future work may help to further quantify and
refine GPTscreenR’s performance, and to expand this approach to other types of
scholarly review and to other stages of the review process.

##  6 References

## References

  * [1] Hilary Arksey and Lisa O’Malley  “Scoping studies: towards a methodological framework”  In _International Journal of Social Research Methodology_ 8.1, 2005, pp. 19–32  DOI: 10.1080/1364557032000119616
  * [2] Kevin E. K. Chai, Robin L. J. Lines, Daniel F. Gucciardi and Leo Ng  “Research Screener: a machine learning tool to semi-automate abstract screening for systematic reviews”  In _Systematic Reviews_ 10.1, 2021, pp. 93  DOI: 10.1186/s13643-021-01635-3
  * [3] Allison Gates et al.  “Performance and usability of machine learning for screening in systematic reviews: a comparative evaluation of three tools”  In _Systematic Reviews_ 8.1, 2019, pp. 278  DOI: 10.1186/s13643-019-1222-2
  * [4] Eddie Guo et al.  “Automated Paper Screening for Clinical Reviews Using Large Language Models”  In _arXiv_ , 2023  DOI: 10.48550/arxiv.2305.00844
  * [5] Brian E. Howard et al.  “SWIFT-Review: a text-mining workbench for systematic review”  In _Systematic Reviews_ 5.1, 2016, pp. 87  DOI: 10.1186/s13643-016-0263-z
  * [6] Hanan Khalil, Daniel Ameen and Armita Zarnegar  “Tools to support the automation of systematic reviews: a scoping review”  In _Journal of Clinical Epidemiology_ 144, 2022, pp. 22–42  DOI: 10.1016/j.jclinepi.2021.12.005
  * [7] Takeshi Kojima et al.  “Large Language Models are Zero-Shot Reasoners”  In _arXiv_ , 2022  DOI: 10.48550/arxiv.2205.11916
  * [8] Michael L. Mortensen et al.  “An exploration of crowdsourcing citation screening for systematic reviews”  In _Research Synthesis Methods_ 8.3, 2017, pp. 366–386  DOI: 10.1002/jrsm.1252
  * [9] OpenAI  “Introducing ChatGPT”, ”https://openai.com/blog/chatgpt”, 2023
  * [10] Micah D. J. Peters et al.  “Updated methodological guidance for the conduct of scoping reviews”  In _JBI Evidence Synthesis_ 18.10, 2020, pp. 2119–2126  DOI: 10.11124/jbies-20-00167
  * [11] Micah D. J. Peters et al.  “Scoping reviews: reinforcing and advancing the methodology and application”  In _Systematic Reviews_ 10.1, 2021, pp. 263  DOI: 10.1186/s13643-021-01821-3
  * [12] Micah DJ Peters et al.  “Chapter 11: Scoping Reviews (2020 version)”  In _JBI Manual for Evidence Synthesis_ JBI, 2020  DOI: 10.46658/jbimes-20-12
  * [13] Mai T. Pham et al.  “A scoping review of scoping reviews: advancing the approach and enhancing the consistency”  In _Research Synthesis Methods_ 5.4, 2014, pp. 371–385  DOI: 10.1002/jrsm.1123
  * [14] R Core Team  “R: A Language and Environment for Statistical Computing”, 2023  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [15] Alec Radford, Karthik Narasimhan, Tim Salimans and Ilya Sutskever  “Improving language understanding by generative pre-training”, ”https://openai.com/research/language-unsupervised”, 2018
  * [16] Ian Shemilt et al.  “Pinpointing needles in giant haystacks: use of text mining to reduce impractical screening workload in extremely large scoping reviews”  In _Research Synthesis Methods_ 5.1, 2014, pp. 31–49  DOI: 10.1002/jrsm.1093
  * [17] Claire Stansfield, James Thomas and Josephine Kavanagh  “‘Clustering’ documents automatically to support scoping reviews of research: a case study”  In _Research Synthesis Methods_ 4.3, 2013, pp. 230–241  DOI: 10.1002/jrsm.1082
  * [18] Anthea Sutton, Mark Clowes, Louise Preston and Andrew Booth  “Meeting the review family: exploring review types and associated information retrieval requirements”  In _Health Information & Libraries Journal_ 36.3, 2019, pp. 202–222  DOI: 10.1111/hir.12276
  * [19] Eugene Syriani, Istvan David and Gauransh Kumar  “Assessing the Ability of ChatGPT to Screen Articles for Systematic Reviews”  In _arXiv_ , 2023
  * [20] The Joanna Briggs Institute  “The Joanna Briggs Institute Reviewers’ Manual 2015: Methodology for JBI Scoping Reviews”  The Joanna Briggs Institute, 2015
  * [21] Jason Wei et al.  “Chain of Thought Prompting Elicits Reasoning in Large Language Models”  In _arXiv_ , 2022  DOI: 10.48550/arxiv.2201.11903
  * [22] William Revelle  “psych: Procedures for Psychological, Psychometric, and Personality Research” R package version 2.3.6, 2023  Northwestern University  URL: https://CRAN.R-project.org/package=psych
  * [23] Shunyu Yao et al.  “Tree of Thoughts: Deliberate Problem Solving with Large Language Models”  In _arXiv_ , 2023  URL: https://arxiv.org/abs/2305.10601


# Fully Kinetic Simulations of Proton-Beam-Driven Instabilities from Parker
Solar Probe Observations

L. Pezzini Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven, Leuven, Belgium  Solar-Terrestrial Centre of
Excellence–SIDC, Royal Observatory of Belgium, Brussels, Belgium  A. N. Zhukov
Solar-Terrestrial Centre of Excellence–SIDC, Royal Observatory of Belgium,
Brussels, Belgium  Skobeltsyn Institute of Nuclear Physics, Moscow State
University, Moscow, Russia  F. Bacchini Centre for mathematical Plasma
Astrophysics, Department of Mathematics, KU Leuven, Leuven, Belgium  Royal
Belgian Institute for Space Aeronomy, Solar-Terrestrial Centre of Excellence,
Brussels, Belgium  G. Arrò Los Alamos National Laboratory, Los Alamos, USA  R.
A. López Research Center in the intersection of Plasma Physics, Matter, and
Complexity ($P^{2}mc$),
Comisión Chilena de Energía Nuclear, Casilla 188-D, Santiago, Chile  A. Micera
Institut für Theoretische Physik, Ruhr-Universität Bochum, Bochum, Germany  M.
E. Innocenti Institut für Theoretische Physik, Ruhr-Universität Bochum,
Bochum, Germany  G. Lapenta Centre for mathematical Plasma Astrophysics,
Department of Mathematics, KU Leuven, Leuven, Belgium

###### Abstract

The expanding solar wind plasma ubiquitously exhibits anisotropic non-thermal
particle velocity distributions. Typically, proton Velocity Distribution
Functions (VDFs) show the presence of a core and a field-aligned beam. Novel
observations made by Parker Solar Probe (PSP) in the innermost heliosphere
have revealed new complex features in the proton VDFs, namely anisotropic
beams that sometimes experience perpendicular diffusion. This phenomenon gives
rise to VDFs that resemble a “hammerhead”. In this study, we use a 2.5D fully
kinetic simulation to investigate the stability of proton VDFs with
anisotropic beams observed by PSP. Our setup consists of a core and an
anisotropic beam populations that drift with respect to each other. This
configuration triggers a proton-beam instability from which nearly parallel
fast magnetosonic modes develop. Our results demonstrate that before this
instability reaches saturation, the waves resonantly interact with the beam
protons, causing significant perpendicular heating at the expense of the
parallel temperature. Furthermore, the proton perpendicular heating induces a
hammerhead-like shape in the resulting VDF. Our results suggest that this
mechanism probably contribute to producing the observed hammerhead
distributions.

Solar wind (1534) — Plasma astrophysics (1261) — Space plasmas (1544)

##  1 Introduction

For decades, spacecraft missions like Helios, Ulysses and Wind have conducted
in situ observations, measuring velocity distribution functions (VDFs) of both
ions and electrons in the solar wind (Marsch et al., 1982; Marsch, 2006).
These observations revealed several features regarding the kinetic properties
of the different particle species in the solar wind plasma. In this
environment, not all plasma species have equal temperatures. Electrons are
generally colder than protons in the fast solar wind but hotter than protons
in the slow solar wind (Montgomery et al., 1968; Newbury et al., 1998). Proton
populations can also present temperature anisotropy with respect to the
direction of the local magnetic field (Bale et al., 2013). According to the
double-adiabatic expansion model, in the absence of heat flux and collisions,
protons in the inner heliosphere would show
$T_{\perp,\mathrm{p}}<T_{\parallel,\mathrm{p}}$ (Chew et al., 1956). Here,
$T_{\perp,\mathrm{p}}$ and $T_{\parallel,\mathrm{p}}$ denote the perpendicular
and parallel temperatures of protons, respectively. However, observations
reveal a predominance of cases where
$T_{\perp,\mathrm{p}}>T_{\parallel,\mathrm{p}}$ in the fast solar wind within
the inner heliosphere, underscoring the need for a mechanism driving
substantial perpendicular heating (Richardson et al., 1995; Matteini et al.,
2007). Proton VDFs are predominantly gyrotropic, often exhibiting skewness,
resulting in a heat flux aligned along the magnetic field, directed away from
the Sun (Marsch et al., 1982). They also often exhibit a magnetic field-
aligned beam in the anti-sunward direction, a secondary proton population
streaming at higher velocities compared to the core one (Asbridge et al.,
1974; Alterman et al., 2018). It has been suggested that the process at the
origin of the proton beam is the reflection of turbulent fluctuations at the
MHD/kinetic spectral break: the greater the wavenumber of the spectral break,
the faster the generated proton beam (Voitenko & Pierrard, 2015). These non-
Maxwellian features are more pronounced in the fast solar wind compared to the
slow solar wind, and their magnitude increases as the distance to the Sun
decreases (Verscharen et al., 2019b). They represent symmetry-breaking
properties capable of inducing plasma instabilities, leading to subsequent
energy transfer processes restoring the equilibrium within the system
(Verniero et al., 2020).

The PSP mission (Fox et al., 2016) aims to study fundamental processes
involved in energy transfer within the solar wind, such as acceleration and
heating processes, to the solar surface. Observations from the Solar Wind
Electron Alpha and Proton (SWEAP) instrument onboard PSP (Kasper et al., 2016)
detected anisotropic beams in the proton VDFs, as reported by Klein et al.
(2021). Moreover, extreme cases of anisotropic beam structures were also
observed, specifically termed “hammerhead” (Verniero et al., 2022) due to
their resemblance to the shape of a hammerhead. Consistently with the
measurements of SWEAP, the FIELDS instrument (Bale et al., 2016) detected a
strong wave activity from fluid MHD scales down to ion kinetic scales (Bowen
et al., 2020; Verniero et al., 2020). The solar wind is a weakly collisional
plasma (Alterman et al., 2018; Verscharen et al., 2019b), thus it is
wave/particle interaction processes, rather than collisions, which act in the
direction of reducing the non-Maxwellian character of VDFs (Kasper et al.,
2002; Klein et al., 2018). The linear Alfvén dispersion relation at MHD scales
connects with the fast-magnetosonic/whistler (FM/W) waves at kinetic scales
(Belcher & Davis, 1971; Verscharen et al., 2019b). These modes are right-
handed (RH) polarized and can be driven by temperature anisotropy with
$T_{\perp}\ll T_{\parallel}$ in a high plasma beta, denoted as $\beta\doteq
8\pi nk_{\mathrm{B}}T/B^{2}\gg 1$ (Gary, 1993; Kennel & Engelmann, 1966).
Here, $n$ represents the number density, $T$ signifies the temperature, $B$
denotes the magnetic field, and $k_{\mathrm{B}}$ stands for the Boltzmann
constant. The observed relative drift speeds between core protons, secondary
proton beams, and $\alpha$-particles (Marsch et al., 1982; Neugebauer et al.,
1996; Steinberg et al., 1996; Podesta & Gary, 2011a, b) constitute another
reservoir of free energy that can excite waves through wave-particle resonant
interaction or cyclotron resonance (Bourouaine et al., 2013; Verscharen et
al., 2013; Ofman et al., 2022, 2023). Damping mechanisms such as the cyclotron
resonant damping (Goldstein et al., 1994; Leamon et al., 1998) dissipate
energy at ion kinetic scales, where particles interact efficiently with
electromagnetic waves (Leamon et al., 1998; Chandran et al., 2011). Through
quasi-linear theory it is possible to detect signatures of cyclotron resonant
damping in proton density phase-space (Hollweg & Isenberg, 2002; Cranmer,
2014; Bowen et al., 2022).

The previously mentioned hammerhead structure (Verniero et al., 2022) is
believed to originate from resonant interaction between parallel-propagating
fluctuations and proton at kinetic scales, similar to those observed in (Klein
et al., 2021). However, its precise origin has yet to be firmly established.
Hybrid numerical simulation can efficiently study non-linear processes
associated with unstable proton-beam VDFs (Ofman et al., 2022, 2023). However,
in these codes, electrons are typically treated as a fluid, which leads to a
loss of their fundamental kinetic nature. To gain insight into the mechanism
behind the origin of the “hammerhead” (Verniero et al., 2022) from a first-
principles approach and thus address the issue of perpendicular heating in
solar wind plasma. To achieve this, we extend and complement the stability
analysis of the proton-beam unstable VDF investigated by Klein et al. (2021)
by exploring its nonlinear regime. In this study, we employ ECsim (Lapenta,
2017; Gonzalez-Herrero et al., 2019; Bacchini, 2023; Croonen, 2024), a semi-
implicit, exactly energy-conserving PIC simulation tool, to investigate how
the dynamics of both electrons and protons impact the development of the
proton hammerhead. This paper is organized as follows: in Section 2, we
present the setup of the PIC simulation in detail. Section 3 showcases the
numerical modeling results alongside linear-theory analyses. Finally, Section
4 contains the discussion and conclusions.

##  2 Numerical setup

The ECsim code (Lapenta, 2017; Gonzalez-Herrero et al., 2019; Bacchini, 2023;
Croonen, 2024) was employed to solve the Vlasov-Maxwell system of equations
for a proton-electron plasma. Our initial parameters are inspired by PSP
measurements made on the 30 of January 2020, as reported by Klein et al.
(2021). We conducted a simulation in a periodic two-dimensional $(x,y)$
Cartesian geometry. At initialization, the electric field is set to zero and
the background magnetic field is homogeneous
$\boldsymbol{B}_{0}=B_{0}\hat{e}_{x}$. Introducing a field-aligned reference
system, the $x$-axis denotes the parallel direction, while the $y$ and $z$
axes represent the perpendicular directions.

The non-Maxwellian protons, denoted by the subscript “p”, have constant number
density $n_{\mathrm{p}}$, are divided into two sub-populations. The core, with
subscript “c”, has number density $n_{\mathrm{c}}/n_{\mathrm{p}}=0.864$ and
the beam, subscript “b”, has number density
$n_{\mathrm{b}}/n_{\mathrm{p}}=0.136$. The core and the beam are both
initialized with a bi-Maxwellian distribution, each characterized by a non-
zero drift velocity along the $x$-axis (parallel direction):

$\displaystyle f_{\mathrm{p}}\left(v_{x},v_{y},v_{\mathrm{z}}\right)=$
$\displaystyle\sum_{i=\mathrm{c},\mathrm{b}}\frac{n_{i}}{\pi^{3/2}w_{x,i}w_{y,i}w_{z,i}}$
(1)
$\displaystyle\times\exp\left(-\frac{(v_{x}-V_{x,i})^{2}}{w_{x,i}^{2}}-\frac{v_{y}^{2}}{w_{y,i}^{2}}-\frac{v_{z}^{2}}{w_{z,i}^{2}}\right).$

In equation (1) the subscript $i$ denotes the proton species, i.e. the core or
the beam. The thermal velocities of the $i$-th proton species are expressed as
$w_{i}\doteq\sqrt{k_{\mathrm{B}}T_{i}/m_{\mathrm{p}}}$, where $m_{\mathrm{p}}$
stands for the proton mass. Initially, these velocities are set to
$w_{x,\mathrm{c}}/c_{\mathrm{A,p}}\approx 0.64$ and
$w_{y,\mathrm{c}}/c_{\mathrm{A,p}}\approx 0.56$ for the core, while for the
beam they are assigned as $w_{x,\mathrm{b}}/c_{\mathrm{A,p}}\approx 1.00$ and
$w_{y,\mathrm{b}}/c_{\mathrm{A,p}}\approx 0.79$. The Alfvén speed for the
protons defined as $c_{\mathrm{A,p}}\doteq B/\sqrt{4\pi
n_{\mathrm{p}}m_{\mathrm{p}}}$ at the initialisation with a magnitude of
$c_{\mathrm{A,p}}/c=0.00291$, where $c$ represents the speed of light in
vacuum. Given that the plasma is described by a bi-Maxwellian distribution,
there exists a symmetry in the perpendicular velocities, expressed as
$w_{y,i}=w_{z,i}$ for all species $i$. The drift velocities of the core and
beam are approximately $V_{x,\mathrm{c}}/c_{\mathrm{A,p}}\approx 0.18$ and
$V_{x,\mathrm{b}}/c_{\mathrm{A,p}}\approx-1.34$, respectively. As a result,
these two populations move towards each other. Electrons are initialized as a
Maxwellian with thermal velocity $w_{\mathrm{e}}/c_{\mathrm{A,p}}\approx 7.12$
with unitary number density $n_{\mathrm{e}}$. The described electron-proton
plasma satisfies the zero-net-current condition
$n_{\mathrm{c}}V_{x,\mathrm{c}}+n_{\mathrm{b}}V_{x,\mathrm{b}}=0$ and the
quasi-neutrality condition $n_{\mathrm{e}}=n_{\mathrm{p}}$.

The proton VDF previously discussed is unstable and expected to produce fast
magnetosonic modes, as demonstrated by the linear analysis conducted in (Klein
et al., 2021). Based on linear theory calculations (for detailed analysis, see
Section 3), the instability would manifest within a specific range of spatio-
temporal scales, namely around those associated with the fastest growing mode.
Consequently, we adjust the spatial and temporal parameters of our simulation
accordingly. The fastest growing mode has a parallel wave number with a value
of $k_{\parallel}/d_{\mathrm{p}}^{-1}\approx 0.5$ and a perpendicular one of
$k_{\perp}/d_{\mathrm{p}}^{-1}\approx 0.03$. Here, $d_{\mathrm{p}}$ is the
proton skin depth, defined as $d_{\mathrm{p}}\doteq c/\omega_{\mathrm{p}}$,
and $\omega_{\mathrm{p}}\doteq\sqrt{4\pi e^{2}n_{\mathrm{p}}/m_{\mathrm{p}}}$
is the proton plasma frequency. This leads to a wavelength in the parallel
direction of
$\lambda_{x}/d_{\mathrm{p}}=2\pi/(k_{\parallel}/d_{\mathrm{p}}^{-1})\approx
13$ and in the perpendicular direction of
$\lambda_{y}/d_{\mathrm{p}}=2\pi/(k_{\perp}/d_{\mathrm{p}}^{-1})\approx 209$.
As a result, the domain length in each direction is selected to host several
wavelength of this mode. This choice yields a rectangular box with dimensions
$L_{x}/d_{\mathrm{p}}=64$ and $L_{y}/d_{\mathrm{p}}=256$. The spatial grid
consists of $256\times 1024$ cells with a spatial increment of $\Delta
x/d_{\mathrm{p}}=\Delta y/d_{\mathrm{p}}=0.25$. Similarly to the spatial
domain, the estimation of the temporal scale is made by employing the growth
rate of the most unstable mode $\gamma_{\mathrm{m}}/\omega_{\mathrm{p}}\approx
0.000121$ or $\gamma_{\mathrm{m}}/\Omega_{\mathrm{p}}\approx 0.042$, where
$\Omega_{\mathrm{p}}\doteq eB_{0}/(m_{\mathrm{p}}c)$ is the proton cyclotron
frequency. The ratio between the proton plasma and the cyclotron frequency is
$\omega_{\mathrm{p}}/\Omega_{\mathrm{p}}\approx 344.04$. The instability
evolves over a time scale
$t/\omega_{\mathrm{p}}^{-1}=2\pi/(\gamma_{\mathrm{m}}/\omega_{\mathrm{p}})\approx
52000$. Therefore, the simulation is run for a total time
$t_{\mathrm{tot}}/\omega_{\mathrm{p}}^{-1}\approx 150000$, to fully capture
the development of the instability until the non-linear phase. The time step
is set to $\Delta t/\omega_{\mathrm{p}}^{-1}=0.5$ to ensure accurate
resolution of the proton gyro-period. We employed $4096$ particles per cell
per species, initially distributed uniformly on the grid. We utilized a
reduced proton-to-electron mass ratio, $m_{\mathrm{p}}/m_{\mathrm{e}}=183.6$,
to narrow the gap between electron and proton dynamical scales, thereby saving
computational time.

##  3 Results

Figure 1: Results for the linear Vlasov theory, corresponding to our initial
setup, in the left column, compared with the simulation results in the right
column. Panel (a) to (f) show the following: the real frequency
$\omega_{\mathrm{r}}/\omega_{\mathrm{p}}$ (a); the imaginary frequency
$\gamma$ (c); the polarization of the wave (e). The right column shows the
following: in panel (b), the power spectrum measured for a cut at $y=L_{y}/2$
is plotted against the frequency at different positions along the $x$-axis;
panel (d) shows the power spectrum at $t/\Omega_{\mathrm{p}}^{-1}\approx 218$,
with some of the instability growth rates contours of panel (c) overlaid, with
the black dot denoting the fastest-growing mode; panel (f) is the hodogram at
$t/\Omega_{\mathrm{p}}^{-1}\approx 218$, where the red triangle shows the
value of the magnetic fluctuation corresponding to the point $x=L_{x}$, see
the text for details on how this is computed.

The results from our numerical simulation are compared with the Vlasov linear
theory in Figure 1. We solved the linear dispersion relation using the DIS-K
solver111The code is publicly available and can be found at
https://github.com/ralopezh/dis-k. for the setup described above (López et
al., 2021; López, 2023). Specifically, the spectrum and polarisation of the
waves developed by the system are analyzed. The real part
$\omega_{\mathrm{r}}/\omega_{\mathrm{p}}$ of the complex frequency , as
obtained from linear theory, is illustrated in Figure 1(a) against
$k_{\parallel}$ and $k_{\perp}$. In Figure 1(b) we plot the power spectrum of
the magnetic fluctuations
$\delta\boldsymbol{B}/B_{0}=(\boldsymbol{B}-\boldsymbol{B}_{0})/B_{0}$
measured along the parallel direction fixing $y=L_{y}/2$, denoted as
$\|\mathcal{F}\left\\{\delta\boldsymbol{B}/B_{0}(y=L_{y}/2)\right\\}\|^{2}$.
Here, $\mathcal{F}$ represent the Fast Fourier Transform (FFT) in time, and
the colorbar indicates the variation in the spatial dimension $x/L_{x}$. The
peak of maximum intensity occurs at
$\omega_{\mathrm{r}}/\omega_{\mathrm{p}}=0.00151$ at all $x/L_{x}$ as pointed
out in panel (b) by the black dashed line, corresponding to
$k_{\parallel}/d_{\mathrm{p}}^{-1}\approx 0.501$ in panel (a). According to
Vlasov linear theory, the imaginary part of the complex frequency, $\gamma$,
represents the growth rate of the unstable modes. Figure 1(c) plots the power
spectrum of $\gamma/\omega_{\mathrm{p}}$ in the $k_{\parallel}$–$k_{\perp}$
plane. In Figure 1(d), we present the analog of panel (c) for the simulated
data. We display the power spectrum of $\delta\boldsymbol{B}/B_{0}$ at
$t/\Omega_{\mathrm{p}}^{-1}\approx 218$. Similarly to panel (b),
$\mathcal{F}\left\\{\delta\boldsymbol{B}/B_{0}\right\\}$ represents the 2D FFT
in the spatial domain of $\delta\boldsymbol{B}/B_{0}$. Black isocontours of
$\gamma/\omega_{\mathrm{p}}$ values from linear theory in panel (c), have been
superimposed to highlight the good agreement between analytic and simulation
results. This specific time has been chosen because it corresponds to the
moment when the magnetic energy reaches its maximum, see Figure 4(a). The
dominant mode (i.e. fastest-growing mode), indicated by the black dot,
propagates mostly along the parallel direction ($x$-axis) with
$k_{\parallel}/d_{\mathrm{p}}^{-1}\approx 0.501$ and
$k_{\perp}/d_{\mathrm{p}}^{-1}\approx 0.0245$, having a magnitude of
$\gamma_{\mathrm{m}}/\Omega_{\mathrm{p}}\approx 0.042$. In Figure 1(e), we
plot the polarisation $P$ of the electric field of the wave as calculated from
linear theory. If $P>0$, then the wave is right-handed (RH) elliptically
polarised, whereas if $P<0$, it is left-handed (LH). In this case, the
polarisation is RH, but it is very close to unity, making the expected
hodogram almost circular. The hodogram in Figure 1(f) provides information
about the polarisation of the wave for our simulation. By the reported time
frame of $t/\Omega_{\mathrm{p}^{-1}}\approx 218$ all unstable modes have
completed their development. A wave polarisation state can be split into two
linearly polarized orthogonal components with respect to the guiding magnetic
field (Narita, 2017; Born & Wolf, 2019). In our case, the two components are
measured as the values, along the $x$-axis at $y=L_{y}/2$, of $\delta
B_{y}/B_{0}$ and $\delta B_{z}/B_{0}$, representing the magnetic fluctuation
in the transverse directions at $t/\Omega_{\mathrm{p}}^{-1}\approx 218$. The
red triangle highlights the value of the magnetic fluctuation at
$x/{d_{\mathrm{p}}}=L_{x}$ and $y/{d_{\mathrm{p}}}=L_{y}/2$. From this point,
a counterclockwise rotation corresponds to an RH circular polarization state,
while a clockwise rotation refers to an LH circular polarization state. In
Figure 1(f) we have an RH circularly polarized wave, which once again aligns
with the linear theory shown in Figure 1(e) and discussed in Klein et al.
(2021).

Figure 2: Magnetic field fluctuations resulting from the PIC simulation. Panel
(a) displays the spectral power of magnetic fluctuations for different
wavenumbers (i.e., different modes), corresponding to those within the region
delimited by the inner contour in Figure 1(d) with
$k_{\parallel}/d_{\mathrm{p}}^{-1}\in[0.3,0.6]$. The calculated fastest-
growing mode, represented by the violet line and highlighted with its
corresponding wavenumber, is compared with the expected maximum growth rate
from the linear theory, indicated by the dashed black line. In Panel (b), a
spatio-temporal analysis is presented, focusing on the cut along $y=L_{y}/2$
of the magnetic fluctuation stacked in time. The black line in panel (b) has a
slope equal to the phase velocity of the wave calculated from linear theory.

In Figure 2(a), the spectral power of magnetic fluctuations has been plotted
for for a number of unstable modes within the range of parallel wavenumbers
$k_{\parallel}/d_{\mathrm{p}}^{-1}\in[0.3,0.6]$ Figure 1(d). Our main focus is
on the mode with the steepest slope, representing the fastest-growing one.
This mode possesses wavenumber values
$k_{\parallel}/d_{\mathrm{p}}^{-1}=0.501$ and
$k_{\perp}/d_{\mathrm{p}}^{-1}=0.0245$ and along with the growth rate
$\gamma_{\mathrm{m}}/\Omega_{\mathrm{p}}=0.042$, as determined by linear
theory. The maximum theoretical growth rate $\gamma_{\mathrm{m}}$ from the
linear analysis of Figure 1(c) is drawn in dashed line as comparison. The
fastest growing mode is exponentially growing during the linear phase of the
instability, occurring from time $t/\Omega_{\mathrm{p}}^{-1}\approx 116$ to
$t/\Omega_{\mathrm{p}}^{-1}\approx 189$, closely matching the predicted growth
rate according to Vlasov linear theory. From the analysis in Figure 1(b)-(c),
we can calculate the phase velocity of the wave corresponding to the most-
unstable mode, denoted as
$v_{\mathrm{ph}}=\omega_{\mathrm{r}}/k_{\parallel}\approx 0.00306$. In Figure
2(b), a spatial and temporal analysis of the magnetic fluctuations
$\delta\boldsymbol{B}/B_{0}(y=L_{y})$ is presented. This plot was obtained
stacking in time the value of $\delta\boldsymbol{B}/B_{0}(y=L_{y})$. The line
we overlaid, with equation $x=v_{\mathrm{ph}}t+b$, where
$b/d_{\mathrm{p}}\approx-130.6$ is the offset, demonstrates that the
perturbation in the magnetic field is propagating at the speed predicted by
the linear analysis.

Figure 3: Fluctuations in the proton density $\rho_{\mathrm{p}}$ in panel (a),
with $\delta\rho_{\mathrm{p}}=\rho_{\mathrm{p}}(t)-\rho_{0}$, and fluctuations
in the magnetic field components $B_{x},B_{y},B_{z}$, in panel (b)-(d), at the
moment of the instability saturation $t/\Omega_{\mathrm{p}}^{-1}\approx 218$.
The amplitude of $\rho_{\mathrm{p}}$ and $B_{x}$ fluctuations, panels (a) and
(b), have been multiplied by a factor 10 for visualisation purposes.

The spatial distribution of fluctuations in proton mass density
$\rho_{\mathrm{p}}$ and magnetic field components is displayed in Figure
3(a)–(d). They are presented at time $t/\Omega_{\mathrm{p}}^{-1}\approx 218$,
corresponding to the moment of instability saturation after the end of the
linear phase, to examine the waves when they are fully developed. The
perturbation in proton density, as visible in panel (a), indicates that the
instability develops as a compressive wave. Additionally, there is a positive
correlation between the fluctuation in proton density and the magnetic field
fluctuation in the $x$-direction, as denoted by both the pattern and intensity
shown in Figure 3(a)-(b). The magnetic fluctuation in the $y$ and $z$
directions, representing the transverse directions, is more intense,
indicating that it carries more energy than that in the parallel direction. We
observe that these wavepackets are predominantly parallel-propagating, and
thus, $\theta\doteq k_{\perp}/k_{\parallel}\approx 0.05$, with $\theta$
representing the angle between the wavevector components $k_{\parallel}$ and
$k_{\perp}$ as calculated from linear theory. Therefore, we can infer that the
developed wave is a fast-magnetosonic/whistler (FM/W) wave.

Figure 4: Evolution of different types of energy and proton temperature ratio.
Panels (a)-(c) illustrate the time evolution of different energy components,
categorised into three distinct phases named I, II, and III, corresponding to
the stationary, linear, and nonlinear phases, respectively. The gray region
identifies phase II, with the dashed-dotted line indicating the inflection
point of the curves, where the conversion of kinetic energy into magnetic
energy peaks at $t_{\mathrm{peak}}/\Omega_{\mathrm{p}}^{-1}\approx 170$. Panel
(a) displays the changes in magnetic and total kinetic energy of electrons and
protons, normalised by the initial total energy of the system. Panel (b)
presents the thermal energies (parallel and perpendicular with respect to the
guide field direction) of both the core and beam populations, along with the
core and beam drift (i.e. bulk) energy. In panel (c), the temporal evolution
of the temperature ratio of the total proton population, encompassing both the
core and beam, is shown. The formal definition of these quantities is provided
in the text.

In Figure 4, we show the energy exchanges within the considered system. In
Figure 4 (a)-(c) we identify three phases of evolution the system is passing
through. Phase I or quasi-stationary phase ends at approximately
$t_{\mathrm{beg}}/\omega_{\mathrm{p}}^{-1}\approx 126$, which coincides with
the moment the instability starts growing with the fastest growing mode.
Subsequently, Phase II or linear phase, in light gray, begins at the end of
phase I and evolves until the instability saturates and then the conversion
from magnetic to kinetic energy reaches its maximum around
$t_{\mathrm{end}}/\omega_{\mathrm{p}}^{-1}\approx 251$ as illustrated in
Figure 2(a). At the end of this phase, the total magnetic and kinetic energies
respectively reach their maximum and minimum values. Moreover, the dash-dotted
line in Figure 4 highlights the inflection point of the energies, which is the
point where the conversion of kinetic energy into magnetic energy peaks at
$t_{\mathrm{peak}}/\omega_{\mathrm{p}}^{-1}\approx 170$. It is followed by
Phase III, named non-linear phase, which starts at
$t_{\mathrm{end}}/\omega_{\mathrm{p}}^{-1}\approx 251$ and continues until the
end of the simulation.

All the energies in Figure 4 we plot the difference in each energy $\delta
E=E(t)-E(t=0)$ and normalized by $U_{0}$, which is the total energy of the
system at initialization, using the formula $\delta
E/U_{0}\doteq\left(E(t)-E(t_{0})\right)/U_{0}$. The quasi-stationary phase is
characterized by a transient of quiescence in which the total magnetic energy
$E^{\mathrm{m}}\doteq\int B^{2}/(8\pi)~{}\mathrm{d}^{3}\boldsymbol{x}$ and the
total kinetic energy of electrons and protons
$E^{\mathrm{k}}=\sum_{s}\left(E^{\mathrm{th}}_{s}+E^{\mathrm{d}}_{s}\right)$
remain nearly constant over time, as shown in Figure 4(a). The kinetic thermal
energy for the species $s$ is represented as $E^{\mathrm{th}}_{s}\doteq
1/2\int\mathrm{Tr}(\boldsymbol{P}_{s})~{}\mathrm{d}^{3}\boldsymbol{x}$, where
$\mathrm{Tr}(\boldsymbol{P}_{s})$ is the trace of the pressure tensor. The
drift (i.e. bulk) energy of the generic species $s$ is defined as
$E^{\mathrm{d}}_{s}\doteq
1/2\int\rho_{s}V_{s}^{2}~{}\mathrm{d}^{3}\boldsymbol{x}$, where
$\boldsymbol{V}_{s}\doteq
1/n_{s}\int\boldsymbol{v}f_{s}~{}\mathrm{d}^{3}\boldsymbol{v}$ represent the
drift velocity of the $s$-th species and $f_{s}$ denotes the distribution
function of the $s$-th species. The pressure tensor is denoted as
$\boldsymbol{P}_{s}\doteq
m_{s}\int(\boldsymbol{v}-\boldsymbol{V}_{s})(\boldsymbol{v}-\boldsymbol{V}_{s})f_{s}~{}\mathrm{d}^{3}\boldsymbol{v}$,
where $\boldsymbol{v}$ denotes the velocity. In a reference system aligned
with the direction
$\hat{\boldsymbol{b}}\doteq\boldsymbol{B}/\lVert\boldsymbol{B}\rVert$ of the
magnetic field, the perpendicular pressure is thus defined as
$p_{\perp,s}\doteq\boldsymbol{P}_{s}:\left(\mathbb{1}-\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\right)/2$,
and the parallel pressure is
$p_{\parallel,s}\doteq\boldsymbol{P}_{s}:\left(\hat{\boldsymbol{b}}\hat{\boldsymbol{b}}\right)$.
The thermal energy for the general species $s$, in the magnetic field-aligned
reference frame, is calculated in the perpendicular direction as
$E^{\mathrm{th}}_{\perp,s}\doteq 1/(2k_{\mathrm{B}})\int
p_{\perp,s}\,\mathrm{d}^{3}\boldsymbol{x}$ and in the parallel direction as
$E^{\mathrm{th}}_{\parallel,s}\doteq 1/(2k_{\mathrm{B}})\int
p_{\parallel,s}\,\mathrm{d}^{3}\boldsymbol{x}$. The components of kinetic
energy in Figure 4(b) are: the beam parallel and perpendicular thermal energy,
respectively $\delta E^{\mathrm{th}}_{\parallel,\mathrm{b}}/U_{0}$ and $\delta
E^{\mathrm{th}}_{\perp,\mathrm{b}}/U_{0}$, core parallel and perpendicular
thermal energy, respectively $\delta
E^{\mathrm{th}}_{\parallel,\mathrm{c}}/U_{0}$ and $\delta
E^{\mathrm{th}}_{\perp,\mathrm{c}}/U_{0}$, and the drift energy of both the
core and the beam $\delta E^{\mathrm{d}}_{\mathrm{c}}/U_{0}$ and $\delta
E^{\mathrm{d}}_{\mathrm{b}}/U_{0}$.

During the non-linear phase, the magnetic energy is converted back into
kinetic energy, as illustrated in Figure 4(a). In Figure 4(b), it is evident
that $\delta E^{\mathrm{th}}_{\perp,\mathrm{c}}/U_{0}$ and $\delta
E^{\mathrm{th}}_{\parallel,\mathrm{c}}/U_{0}$ do not change significantly.
Meanwhile, the $\delta E^{\mathrm{th}}_{\perp,\mathrm{b}}/U_{0}$ is steadily
increasing, while $\delta E^{\mathrm{d}}_{\mathrm{b}}/U_{0}$ is consistently
decreasing. During this phase, the bulk energy stored in the beam’s drift
motion is converted into the thermal perpendicular energy of the beam. In
Figure 4(c), the time evolution of the proton temperature ratio, defined as
$T_{\perp,\mathrm{p}}/T_{\parallel,\mathrm{p}}=p_{\perp,\mathrm{p}}/p_{\parallel,\mathrm{p}}$
and integrated across the spatial domain, is presented. The perpendicular and
parallel pressures for the total population of protons are calculated from
from the pressure tensor $\boldsymbol{P}_{\mathrm{p}}\doteq
m_{\mathrm{p}}\int(\boldsymbol{v}-\boldsymbol{V}_{\mathrm{cm},\mathrm{p}})(\boldsymbol{v}-\boldsymbol{V}_{\mathrm{cm},\mathrm{p}})f_{\mathrm{p}}\,\mathrm{d}^{3}\boldsymbol{v}$,
where the velocity of the center of mass of the proton population is
$\boldsymbol{V}_{\mathrm{cm},\mathrm{p}}\doteq\left(n_{\mathrm{c}}\boldsymbol{V}_{\mathrm{c}}+n_{\mathrm{b}}\boldsymbol{V}_{\mathrm{b}}\right)/\left(n_{\mathrm{c}}+n_{\mathrm{b}}\right)$.
From these results, we observe that the pronounced outcome of the considered
instability is the increase in proton temperature ratio, redistributing the
energy from the parallel to the perpendicular direction. In conclusion,
studying the energy exchanges of the system in Figure 4(a)-(c) reveals how
energy is redistributed inside the system. The two primary sources of free
energy are the beam temperature ratio and the relative drifting speed of the
core and beam. Therefore, based on this analysis, we can conclude that the
system exhibits an instability which we term as “beam-firehose instability”
because it is powered not only by thermal anisotropy but also by the drifting
population, resembling a beam-plasma type of instability.

Figure 5: Evolution of the proton VDF $f_{\mathrm{p}}(v_{x},v_{y},t)$ (left
column) and its parallel cut $f_{\mathrm{p}}(v_{x},v_{y}=0,t)$ (right column).
Panels (a)–(b) represent the quasi-stationary phase of the proton-beam
instability instability. Panel (a) shows the initial distribution of protons
$f_{\mathrm{p}}(v_{x},v_{y},t)$, whereas panel (b) shows the horizontal cut
$f_{\mathrm{p}}(v_{x},v_{y}=0,t)$ (with the core and beam sub-populations
indicated by different lines), at $t/\Omega_{\mathrm{p}}^{-1}=0$. Panels
(c)–(d) represent the phase space of protons during the linear phase at
$t/\Omega_{\mathrm{p}}^{-1}\approx 218$. The proton VDF during the non-linear
phase at $t/\Omega_{\mathrm{p}}^{-1}\approx 436$ is shown in panels (e) and
(f). The resonant velocity of the system,
$v_{\mathrm{res}}=(\omega_{\mathrm{r}}+\Omega_{\mathrm{p}})/k_{\parallel}$, is
shown by the blue line in panel (f).

Figure 5 shows the time evolution of the core, beam, and total proton VDF,
denoted as $f_{\mathrm{p}}(v_{x},v_{y},t)$, in the left column, as well as a
one-dimensional cut $f_{\mathrm{p}}(v_{x},v_{y}=0,t)$, in the right column, to
provide further evidence of the wave-particle resonant interaction. Panels (a)
and (b) in Figure 5 represent the initial conditions of the simulation,
corresponding to the beginning of the quasi-stationary phase of the
instability. Subsequently, panels (c) and (d) in Figure 5 at
$t/\Omega_{\mathrm{p}}^{-1}\approx 218$ illustrate the linear phase, during
which the beam population undergoes resonant interaction with the excited fast
magnetosonic waves. Figure 5(e)-(f) represents the proton VDF during the non-
linear phase at $t/\Omega_{\mathrm{p}}^{-1}\approx 436$. A hint of the
“hammerhead-like” structure (extensions of the beam towards the values of
higher $v_{y}$), initially observed forming during the linear phase in Figure
5(c) and (d), becomes more pronounced during the non-linear phase in Figure
5(e)–(f). In Figure 5(d) and (f), we have shown a vertical line indicating the
resonant velocity
$v_{\mathrm{res}}\doteq(\omega_{\mathrm{r}}-n\Omega_{\mathrm{p}})/k_{\parallel}$,
where $n$ is any integer, for the fast-magnetosonic wave quantities:
$k_{\parallel}\approx 0.5$, $\Omega_{\mathrm{p}}\approx 0.0029$,
$\omega_{\mathrm{r}}\approx 0.0015$. Cyclotron resonance with ${n}=-1$, which
corresponds to the backward-propagating Alfvén wave and the forward-
propagating fast-magnetosonic mode (Klein et al., 2021), has been assumed.
Comparing this to the distribution in Figure 5(b), we observe a decrease in
the number of particles for velocities $v>v_{\mathrm{res}}$ which identifies a
phase space depletion region. Conversely, we observe an increase in the number
of particles for velocities $v<v_{\mathrm{res}}$ creating a phase space pile-
up region.

Figure 6: Proton VDF $f_{\mathrm{p}}(v_{x},v_{y},t)$ at time
$t/\Omega_{\mathrm{p}}^{-1}\approx 436$ normalised by its value at the initial
time. The blue line shows the parallel resonant velocity $v_{\mathrm{r}}$. The
black dashed circle centered at $(v_{x},v_{y})=(v_{\mathrm{ph}},0)$, with
$v_{\mathrm{ph}}=\omega_{\mathrm{r}}/k_{\parallel}$ the parallel phase
velocity of the waves, is the diffusion path due to the resonant interaction.

Here, we will investigate the nature of proton diffusion through wave-particle
resonant interaction. In Figure 6, the proton VDF is shown at time
$t/\Omega_{\mathrm{p}}^{-1}=436$, normalised by its value at the initial time.
The normalization allows us to see small changes in the VDF more easily. In
the white regions, no significant particle diffusion occurs, while in the blue
regions, there are fewer particles than in the initial state, and the red
regions are where particles have accumulated. The blue line here indicates the
parallel resonant velocity, while the black circle,
$(v_{x}-v_{\mathrm{ph}})^{2}+v_{y}^{2}=1$, represents an iso-line defining the
path in velocity phase space where particle energy remains constant. Here,
$v_{\mathrm{ph}}=\omega_{\mathrm{r}}/k_{\parallel}$ denotes the parallel phase
velocity of the waves. Particles diffuse in phase space along trajectories
tangent to the iso-line of constant energy. Resonant particles experience a
gain in kinetic energy, manifested as an increase in $v_{x}^{2}+v_{y}^{2}$.
This acquired energy is extracted from the waves at $\omega_{\mathrm{res}}$,
resulting in their damping. As illustrated in Figure 6, a VDF depletion region
(blue) is observed for velocities $v>v_{\mathrm{res}}$, while a phase space
pile-up region occurs for $v<v_{\mathrm{res}}$.

##  4 Discussion and Conclusions

We presented the first fully kinetic study of unstable proton VDFs from the
recent PSP observation of 30 January 2020. We employed the ECsim code, a semi-
implicit exactly energy conserving PIC (Lapenta, 2017; Gonzalez-Herrero et
al., 2019; Bacchini, 2023; Croonen, 2024) in order to confirm and extend to
the non-linear regime the previous analysis by Klein et al. (2021), studying
in detail the energy exchanges between the different particle species. For
simplicity’s sake, we model electrons as a Maxwellian distribution function,
notwithstanding the very non-thermal nature of the electron VDF (Micera et
al., 2020, 2021; Verscharen et al., 2022). The proton population was modeled
using a double-component bi-Maxwellian distribution, as described in (Klein et
al., 2021). In this representation the anisotropic proton core and proton beam
species drift along a background magnetic field. Proton beams are known to
drive electromagnetic proton-beam instabilities, leading to the generation of
waves, as described by Gary et al. (1984). In this case, the system develops a
proton-beam instability with maximum growth rate
$\gamma_{\mathrm{m}}/\omega_{\mathrm{p}}=0.000121$ and wavenumbers
$k_{\parallel}/d_{\mathrm{p}}^{-1}=0.5$ and
$k_{\perp}/d_{\mathrm{p}}^{-1}=0.03$. This corresponds to a fast-magnetosonic
wave which is RH circularly polarised and compressible, in good agreement with
the results found by Klein et al. (2021).

Based on our analysis of the temporal evolution of the energy, we can
summarise a few key points. The proton core sub-population does not contribute
significantly to the dynamics of the system, which is primarily dominated by
the beam sub-population. The thermal energy of the beam protons is transferred
from the parallel to the perpendicular direction. The parallel thermal energy
of the proton beam increases at the expense of the proton beam drift energy,
which constitutes the main source of free energy. During the nonlinear phase,
the beam drift energy is contributing to the increase in perpendicular thermal
energy of the beam. The notable result of the wave-particle interaction is an
increase in proton temperature ratio, redistributing energy from the parallel
to the perpendicular direction.

In our simulated proton VDFs, starting from the linear phase, we note a
reduction in the number of particles for velocities exceeding the resonant
velocity $v_{\mathrm{res}}$, indicating a phase space depletion region.
Conversely, there is an increase in the particle count for velocities below
$v_{\mathrm{res}}$, manifesting as a phase space pile-up region. Our
simulation indicates that the fast magnetosonic waves are inducing resonance
in the proton beam population. Consequently, these protons experience
diffusion in phase space, leading to perpendicular heating effects as they
move into less densely populated regions, resulting in larger perpendicular
velocity (Verscharen et al., 2019a). The diffusion process initiates in the
linear phase and continues into the nonlinear phase, where it saturates,
gradually diminishing in intensity.

Cyclotron resonance plays a significant role in shaping magnetic field spectra
and distribution functions observed in solar wind (Goldstein et al., 1994;
Leamon et al., 1998). The nonlinear evolution of the fast magnetosonic mode is
significantly influenced by cyclotron resonance, which plays a crucial role in
inducing perpendicular proton heating and redistributing energy from the
parallel to the perpendicular direction. This resonance process effectively
mitigates the free energy associated with the parallel anisotropy,
constituting a new channel of regulation of the solar wind temperature
anisotropy. We demonstrated that wave-particle resonant interaction can modify
the proton beam VDF producing a hint of a “hammerhead” distribution. This
supports the hypothesis regarding the origin of the proton hammerhead (i.e.
strongly asymmetric beam) VDF proposed by Verniero et al. (2022). The width of
the hammerhead-like feature we see in the proton VDFs, induced by the resonant
interaction between fast magnetosonic waves and protons, is smaller with
respect to observations Verniero et al. (2022). We speculate that the presence
of turbulence, which is ubiquitous in the solar wind, could provide an
additional source of energy to be dissipated via the interaction with the
proton beam, eventually producing a more pronounced hammerhead. It is known
that turbulence drives ion heating (Servidio et al., 2012; Matthaeus et al.,
2020; Arrò et al., 2022), inducing anisotropic deformations on Maxwellian
particle distributions (Servidio et al., 2015; Arrò et al., 2023). However,
the interplay between turbulence and non-Maxwellian distributions, like those
studied in the present work, has been poorly investigated and could
significantly impact dissipation. This point will be addressed in future
works.

The paper by Ofman et al. (2022) presents a comparison between 2D and 3D
hybrid simulations of ion beam instability in solar wind plasma, based on
observations from the PSP near perihelia. The study demonstrates that the 3D
simulation does not exhibit additional effects compared to the 2D simulation,
indicating that the third dimension is not necessary to study proton-beam-
driven instabilities. However, solar wind protons may be more accurately
described by a kappa distribution function instead of a Maxwellian
distribution (Livadiotis & McComas, 2013; Nicolaou et al., 2020). A natural
progression of this work could involve extending the simulation with a study
of this system modeling the proton beam using a kappa distribution. However,
these aspects are beyond the scope of this paper and are left for future
studies.

## Acknowledgements

The computational resources and services used in this work were partially
provided by the VSC (Flemish Supercomputer Center), funded by the Research
Foundation Flanders (FWO) and the Flemish Government – department EWI. L.P.
acknowledges support from a PhD grant awarded by the Royal Observatory of
Belgium and the Fonds voor Wetenschappelijk Onderzoek (FWO), with id number
11PCB24N. A.N.Z. thanks the Belgian Federal Science Policy Office (BELSPO) for
the provision of financial support in the framework of the PRODEX Programme of
the European Space Agency (ESA) under contract number 4000136424. F.B.
acknowledges support from the FED-tWIN programme (profile Prf-2020-004,
project “ENERGY”) issued by BELSPO, and from the FWO Junior Research Project
G020224N granted by the Research Foundation – Flanders (FWO). A.M. is
supported by the Deutsche Forschungsgemeinschaft (German Science Foundation;
DFG) project 497938371. M.E.I. acknowledges support from DFG within the
Collaborative Research Center SFB1491.

## References

  * Alterman et al. (2018) Alterman, B. L., Kasper, J. C., Stevens, M. L., & Koval, A. 2018, The Astrophysical Journal, 864, 112, doi: 10.3847/1538-4357/aad23f
  * Arrò et al. (2022) Arrò, G., Califano, F., & Lapenta, G. 2022, Astronomy & Astrophysics, 668, A33
  * Arrò et al. (2023) Arrò, G., Pucci, F., Califano, F., Innocenti, M. E., & Lapenta, G. 2023, The Astrophysical Journal, 958, 11
  * Asbridge et al. (1974) Asbridge, J. R., Bame, S. J., & Feldman, W. C. 1974, Solar Physics, 37, 451, doi: 10.1007/BF00152503
  * Bacchini (2023) Bacchini, F. 2023, The Astrophysical Journal Supplement Series, 268, 60, doi: 10.3847/1538-4365/acefba
  * Bale et al. (2013) Bale, S. D., Pulupa, M., Salem, C., Chen, C. H. K., & Quataert, E. 2013, The Astrophysical Journal, 769, L22, doi: 10.1088/2041-8205/769/2/L22
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Science Reviews, 204, 49, doi: 10.1007/s11214-016-0244-5
  * Belcher & Davis (1971) Belcher, J. W., & Davis, L. 1971, Journal of Geophysical Research, 76, 3534, doi: 10.1029/JA076i016p03534
  * Born & Wolf (2019) Born, M., & Wolf, E. 2019, Principles of Optics, seventh anniversary edition, 60th anniversary of first edition, 20th anniversary of seventh edition edn. (Cambridge, United Kingdom: Cambridge University Press)
  * Bourouaine et al. (2013) Bourouaine, S., Verscharen, D., Chandran, B. D. G., Maruca, B. A., & Kasper, J. C. 2013, The Astrophysical Journal, 777, L3, doi: 10.1088/2041-8205/777/1/L3
  * Bowen et al. (2020) Bowen, T. A., Mallet, A., Huang, J., et al. 2020, The Astrophysical Journal Supplement Series, 246, 66, doi: 10.3847/1538-4365/ab6c65
  * Bowen et al. (2022) Bowen, T. A., Chandran, B. D. G., Squire, J., et al. 2022, Physical Review Letters, 129, 165101, doi: 10.1103/PhysRevLett.129.165101
  * Chandran et al. (2011) Chandran, B. D. G., Dennis, T. J., Quataert, E., & Bale, S. D. 2011, The Astrophysical Journal, 743, 197, doi: 10.1088/0004-637X/743/2/197
  * Chew et al. (1956) Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 236, 112, doi: 10.1098/rspa.1956.0116
  * Cranmer (2014) Cranmer, S. R. 2014, The Astrophysical Journal Supplement Series, 213, 16, doi: 10.1088/0067-0049/213/1/16
  * Croonen (2024) Croonen, J. 2024, The Astrophysical Journal Supplement Series
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Science Reviews, 204, 7, doi: 10.1007/s11214-015-0211-6
  * Gary (1993) Gary, S. P. 1993, Theory of Space Plasma Microinstabilities
  * Gary et al. (1984) Gary, S. P., Smith, C. W., Lee, M. A., Goldstein, M. L., & Forslund, D. W. 1984, The Physics of Fluids, 27, 1852, doi: 10.1063/1.864797
  * Goldstein et al. (1994) Goldstein, M. L., Roberts, D. A., & Fitch, C. A. 1994, Journal of Geophysical Research: Space Physics, 99, 11519, doi: 10.1029/94JA00789
  * Gonzalez-Herrero et al. (2019) Gonzalez-Herrero, D., Micera, A., Boella, E., Park, J., & Lapenta, G. 2019, Computer Physics Communications, 236, 153, doi: 10.1016/j.cpc.2018.10.026
  * Hollweg & Isenberg (2002) Hollweg, J. V., & Isenberg, P. A. 2002, Journal of Geophysical Research: Space Physics, 107, doi: 10.1029/2001JA000270
  * Kasper et al. (2002) Kasper, J. C., Lazarus, A. J., & Gary, S. P. 2002, Geophysical Research Letters, 29, doi: 10.1029/2002GL015128
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Science Reviews, 204, 131, doi: 10.1007/s11214-015-0206-3
  * Kennel & Engelmann (1966) Kennel, C. F., & Engelmann, F. 1966, The Physics of Fluids, 9, 2377, doi: 10.1063/1.1761629
  * Klein et al. (2018) Klein, K. G., Alterman, B. L., Stevens, M. L., Vech, D., & Kasper, J. C. 2018, Physical Review Letters, 120, 205102, doi: 10.1103/PhysRevLett.120.205102
  * Klein et al. (2021) Klein, K. G., Verniero, J. L., Alterman, B., et al. 2021, The Astrophysical Journal, 909, 7, doi: 10.3847/1538-4357/abd7a0
  * Lapenta (2017) Lapenta, G. 2017, Journal of Computational Physics, 334, 349, doi: 10.1016/j.jcp.2017.01.002
  * Leamon et al. (1998) Leamon, R. J., Matthaeus, W. H., Smith, C. W., & Wong, H. K. 1998, The Astrophysical Journal, 507, L181, doi: 10.1086/311698
  * Livadiotis & McComas (2013) Livadiotis, G., & McComas, D. J. 2013, Space Science Reviews, 175, 183, doi: 10.1007/s11214-013-9982-9
  * López et al. (2021) López, R., Shaaban, S., & Lazar, M. 2021, Journal of Plasma Physics, 87, 905870310, doi: 10.1017/S0022377821000593
  * López (2023) López, R. A. 2023, Ralopezh/Dis-k: First Public Version, [object Object], doi: 10.5281/ZENODO.8184896
  * Marsch (2006) Marsch, E. 2006, Living Reviews in Solar Physics, 3, 1, doi: 10.12942/lrsp-2006-1
  * Marsch et al. (1982) Marsch, E., Mühlhäuser, K.-H., Schwenn, R., et al. 1982, Journal of Geophysical Research: Space Physics, 87, 52, doi: 10.1029/JA087iA01p00052
  * Matteini et al. (2007) Matteini, L., Landi, S., Hellinger, P., et al. 2007, Geophysical Research Letters, 34, 2007GL030920, doi: 10.1029/2007GL030920
  * Matthaeus et al. (2020) Matthaeus, W. H., Yang, Y., Wan, M., et al. 2020, The Astrophysical Journal, 891, 101, doi: 10.3847/1538-4357/ab6d6a
  * Micera et al. (2021) Micera, A., Zhukov, A. N., López, R. A., et al. 2021, The Astrophysical Journal, 919, 42, doi: 10.3847/1538-4357/ac1067
  * Micera et al. (2020) —. 2020, The Astrophysical Journal Letters, 903, L23, doi: 10.3847/2041-8213/abc0e8
  * Montgomery et al. (1968) Montgomery, M. D., Bame, S. J., & Hundhausen, A. J. 1968, Journal of Geophysical Research, 73, 4999, doi: 10.1029/JA073i015p04999
  * Narita (2017) Narita, Y. 2017, Nonlinear Processes in Geophysics, 24, 203, doi: 10.5194/npg-24-203-2017
  * Neugebauer et al. (1996) Neugebauer, M., Goldstein, B. E., Smith, E. J., & Feldman, W. C. 1996, Journal of Geophysical Research: Space Physics, 101, 17047, doi: 10.1029/96JA01406
  * Newbury et al. (1998) Newbury, J. A., Russell, C. T., Phillips, J. L., & Gary, S. P. 1998, Journal of Geophysical Research: Space Physics, 103, 9553, doi: 10.1029/98JA00067
  * Nicolaou et al. (2020) Nicolaou, G., Livadiotis, G., & Wicks, R. T. 2020, Entropy, 22, 212, doi: 10.3390/e22020212
  * Ofman et al. (2022) Ofman, L., Boardsen, S. A., Jian, L. K., Verniero, J. L., & Larson, D. 2022, The Astrophysical Journal, 926, 185, doi: 10.3847/1538-4357/ac402c
  * Ofman et al. (2023) Ofman, L., Boardsen, S. A., Jian, L. K., et al. 2023, The Astrophysical Journal, 954, 109, doi: 10.3847/1538-4357/acea7e
  * Podesta & Gary (2011a) Podesta, J. J., & Gary, S. P. 2011a, The Astrophysical Journal, 742, 41, doi: 10.1088/0004-637X/742/1/41
  * Podesta & Gary (2011b) —. 2011b, The Astrophysical Journal, 734, 15, doi: 10.1088/0004-637X/734/1/15
  * Richardson et al. (1995) Richardson, J. D., Paularena, K. I., Lazarus, A. J., & Belcher, J. W. 1995, Geophysical Research Letters, 22, 325, doi: 10.1029/94GL03273
  * Servidio et al. (2012) Servidio, S., Valentini, F., Califano, F., & Veltri, P. 2012, Physical Review Letters, 108, 045001, doi: 10.1103/PhysRevLett.108.045001
  * Servidio et al. (2015) Servidio, S., Valentini, F., Perrone, D., et al. 2015, Journal of Plasma Physics, 81, 325810107, doi: 10.1017/S0022377814000841
  * Steinberg et al. (1996) Steinberg, J. T., Lazarus, A. J., Ogilvie, K. W., Lepping, R., & Byrnes, J. 1996, Geophysical Research Letters, 23, 1183, doi: 10.1029/96GL00628
  * Verniero et al. (2020) Verniero, J. L., Larson, D. E., Livi, R., et al. 2020, The Astrophysical Journal Supplement Series, 248, 5, doi: 10.3847/1538-4365/ab86af
  * Verniero et al. (2022) Verniero, J. L., Chandran, B. D. G., Larson, D. E., et al. 2022, The Astrophysical Journal, 924, 112, doi: 10.3847/1538-4357/ac36d5
  * Verscharen et al. (2013) Verscharen, D., Bourouaine, S., & Chandran, B. D. G. 2013, The Astrophysical Journal, 773, 163, doi: 10.1088/0004-637X/773/2/163
  * Verscharen et al. (2019a) Verscharen, D., Chandran, B. D. G., Jeong, S.-Y., et al. 2019a, The Astrophysical Journal, 886, 136, doi: 10.3847/1538-4357/ab4c30
  * Verscharen et al. (2019b) Verscharen, D., Klein, K. G., & Maruca, B. A. 2019b, Living Reviews in Solar Physics, 16, 5, doi: 10.1007/s41116-019-0021-0
  * Verscharen et al. (2022) Verscharen, D., Chandran, B. D. G., Boella, E., et al. 2022, Frontiers in Astronomy and Space Sciences, 9, 951628, doi: 10.3389/fspas.2022.951628
  * Voitenko & Pierrard (2015) Voitenko, Y., & Pierrard, V. 2015, Solar Physics, 290, 1231, doi: 10.1007/s11207-015-0661-8